Сплайновая (кусочная) регрессия

* Тема: Сплайновая регрессия в SPSS (кусочные полиномы).
* Ключевые слова: сплайн, кусочно-линейная регрессия, CNLR, NLR, узел, knot, spline, piecewise, разрыв, кубический, квадратичный, полином, нелинейный.
* Опубликован: в SPSSX-L 03.10.2001, перевод: 05.11.2008.
* Автор: корп. SPSS / David Matheson.
* Перевод: А. Балабанов.
* Размещение: http://www.spsstools.ru/Syntax/RegressionRepeatedMeasure/PiecewiseRegression.txt (.sps).
* Проверено: SPSS 15.0.0.


* (Вопрос) Меня интересует, работал ли кто-нибудь в SPSS с моделями кусочной регрессии. 
  У меня есть данные, где наряду с линейными взаимосвязями, присутствуют также и разрывы
  (скачки) функции, которые я хотел бы учесть в модели.


* (Ответ) Размещён в SPSSX-L 03.10.2001. Автор: David Matheson (техническая поддержка SPSS)

Я помещаю пару решений из SPSS AnswerNet. Первое решение связано с моделями кусочной
регрессии, когда положение узлов (известно). Второе решение - с ситуацией, когда положение
узловых точек оценивается по имеющимся данным. Хотя я так понимаю, что в вашей ситуации
есть трудность и с оценкой нахождения узла, и с оценкой скачка функции в узле. См. дополнительную
информацию о подобных моделях в Seber (ссылка дана во втором решении).


Решение 100000148: Сплайновая регрессия (a.k.a. кусочный полином).

Вопрос.
Как можно в SPSS устроить подгонку регрессионной модели, для
которой я ожидаю изменения коэффициента наклона в большую или
меньшую сторону (либо даже изменения знака коэффициента) при
некотором конкретном значении предиктора. Либо если я ожидаю, что
все значения отклика изменятся на некоторую константу после прохождения
предиктором некоторого значения?


Ответ.
Регрессионные модели, в которых функция меняется в одной или нескольких
точках на шкале значений предикторов, называются сплайнами, либо кусочными
полиномами, а сами точки переменения характера функции называются узловыми точками (узлами).
Если аналитику известны эти узлы, подгонка сплайнов без проблем осуществляется
в процедуре REGRESSION. Свойства этих сплайнов описаны ниже, а также
дано описание их подгонки в SPSS.
Существование сплайновой модели предполагается в случае, если аналитик
ожидает, что связь между предиктором и откликом меняется в некоторых точках
из ряда значений предиктора. Такое изменение в узловой точке может выражаться
в изменении формы связи (например, от линейной к квадратичной), в добавлении или
вычитании некоторой константы для прогноза отклика справа от узловой точки предиктора, 
либо просто в изменении наклона, ускорения и т.д. регрессионной функции. 
Допустим, имеется предиктор XA и отклик Y, и мы полагаем, что функция регрессии
меняется при достижении XA некоторого известного значения K1. Допустим также,
что мы вычислили переменные XA2 и XA3 (соответственно, квадраты и кубы значений XA)
для того, чтобы иметь возможность проверить кубическую модель.
Подгонка сплайна в этом случае осуществляется посредством создания дополнительного
набора предикторных переменных, от XB0 до XB3, которые все равны нулю если 
XA <= K1. Эти переменные предназначены для оценки изменения, соответственно, константы, 
линейного, квадратичного и кубического параметров модели для значений XA > K1.
model.
Если XA > K1, то:   XB0=1 ;
                    XB1=(XA-K1) ;
                    XB2=(XA-K1)**2 ;
                    XB3=(XA-K1)**3 .
Включение в модель XB0 позволяет учесть разрыв (скачок) функции регрессии
в узловой точке. Отсутствие XB0 будет означать, что регрессия может поменять
своё направление в узловой точке, но два куска линии (до и после K1) все
же будут соединены в узле. Заметим, что некоторые авторы определяют фиктивную
переменную XB0 как 1 для значений XA, больших, ИЛИ РАВНЫХ K1, другие же - 
для значений XA > K1, как это сделано выше. Правильный вариант нужно выбирать, исходя
из смысловой наполненности модели: сути того процесса, который моделируется. Именно
от этого зависит, будет ли ваша функция в узле непрерывной справа (XB0=1, если XA=K1), 
или непрерывной слева (XB=0,если XA=K1). Для всех элементов XB высших порядков значение
соответствующей XB-переменной в любом случае будет равно нулю в точке K1.

Если мы желаем проверить модель, в которой регрессия непрерывна в K1 и линейна
на обоих участках справа и слева от K1, а изменяется только коэффициент
наклона, мы включим  в качестве предикторов переменные XA и XB1. Значимая t-статистика
для коэффициента XB1 при наличии XA будет означать, что линия регрессии
значимо меняет наклон в узловой точке. Знак и величина коэффициента перед XB1, в таком
случае, укажет на направление и величину данного изменения.

При выборе сплайна существует компромисс между гладкостью функции в узловых
точках и достижением наилучшей подгонки к данным. Для регрессионной функции
порядка R максимальная гладкость достигается при уравнивании всех производных
порядков ниже R с обеих сторон узловой точки. Это ограничение достигается добавлением
лишь одного параметра из набора XB: параметра R-ой степени. Например, отклик Y
может прогнозироваться с помощью переменных XA, XA2, XA3 и XB3. То есть, при прохождении
узловой точки корректируется только кубическая составляющая сплайна. Если подобная модель
плохо подходит под данные, либо не соответствует теории, аналитик может снять одно или
несколько ограничений на непрерывность, введя в модель параметры (XA-K1), т.е. XB-переменные
низших порядков. 
Прекрасное введение в параметризацию и проверку сплайновых моделей можно найти в 
статье Патриции Смит: Patricia Smith, "Splines as Useful and Convenient
Statistical Tools" в журнале American Statistician, 1979 г. Смит приводит также обсуждение
ограничений, необходимых, когда функция после узла является полиномом более низкого порядка, 
чем до него.

Следует иметь в виду, что с увеличением количества узлов и параметров регрессионной модели
исследователь быстро сталкивается с проблемой мультиколлинеарности.


Ниже приведены 4 примера определения и проверки сплайновых моделей.
В каждом примере полагается, что переменные Y и XA известны и присутствуют в 
активном файле данных. Синтаксис команды REGRESSION в примерах достаточно
скудный, т.е. без построения графиков остатков или другой диагностической информации,
что сделано ради фокуса на аспектах, специфичных для подгонки сплайнов.
(Но, разумеется, использование  графиков, диагностики коллинеарности и выбросов, 
должно быть постоянно в арсенале практических методов аналитика при построении
регрессий с помощью процедуры REGRESSION для оценки корректности модели, 
выполнения предположений и контроля эффекта нехарактерных значений.)

ПРИМЕР 1
В первом примере линейная модель имеет 2 узла, в XA=15 и XA=25, и 
является непрерывной в обеих узловых точках.

COMPUTE xb1 = xa - 15.
COMPUTE xc1 = xa - 25.
RECODE xb1 xc1 (lo thru 0 = 0).
REGRESSION VARIABLES = y xa xb1 xc1
      /STATISTICS = DEFAULT
      /DEP = y
      /METHOD ENTER xa
      /METHOD ENTER xb1
      /METHOD ENTER xc1.

XB1 и XC1 вводятся разновременно, чтобы проверить изменение наклона
в первом узле вне зависимости от второго узла, рассмотрев гипотезу
о том, что XB1 = 0 при наличии XA. Затем подобная же гипотеза проверяется
для коэффициента перед XC1, уже при наличии XA и XB1.

ПРИМЕР 2
Во втором примере есть только один узел в XA=15. Сплайн непрерывен в
узле, но порядок полинома меняется. До узла взаимосвязь линейная, после
узла - квадратичная. Для моделирования этого эффекта созданы переменные
XB1 и XB2. Заметим, что квадратичный эффект для куска сплайна после K1
проверяется при наличии линейной составляющей для этого куска. Проверка
полиномов более высокого порядка при наличии всех полиномов более низкого 
порядка является стандартной практикой в подгонке полиномов.
От такой практики, однако, отходят, например, в ситуации, описанной в Примере 4,
где все составляющие высоких порядков для куска XB введены в модель в первую очередь.
В последнем случае мы проверяли выполнение ограничений на непрерывность, и добавление
каждого нового XB-параметра в модель улучшало качество подгонки регрессии, но ухудшало
её непрерывность.

COMPUTE xb1 = xa - 15.
RECODE xb1 (LO THRU 0 = 0).
COMPUTE xb2 = xb1*xb1.
REGRESSION VARIABLES = y xa xb1 xb2
     /STATISTICS = DEFAULT CHA
     /DEP = y
     /METHOD ENTER xa
     /METHOD ENTER xb1
     /METHOD ENTER xb2.

ПРИМЕР 3
В третьем примере сплайн имеет разрыв в узле XA=15.
Модель предполагает при этом, что функция квадратична до и после узла.
Сплайн непрерывен справа в узле, т.е., если XA=15, то XB0=1
(функция делает "скачок" в точке 15, а не после 15).

COMPUTE XA2 = XA*XA.
COMPUTE XB0 = (XA GE 15).
COMPUTE XB1 = XB0 * (XA - 15).
COMPUTE XB2 = XB1 * XB1.
REGRESSION VARIABLES = y xa xa2 xb0 xb1 xb2
      /STATISTICS = DEFAULT CHA
      /DEP = y
      /METHOD ENTER xa xa2
      /METHOD ENTER xb0 xb1 xb2.

Здесь в начале проверяется единая квадратичная функция. Затем, со введением
XB-переменных во второй блок модели, проверяется потребность в построении сплайна.
Если изменение R-квадрат не является значимым, это означает, что единая квадратичная
функция одинаково описывает взаимосвязь X и Y на всем ряде значений переменной X.


ПРИМЕР 4
В четвертом примере существует один узел в точке XA=15. Подгоняется кубический
сплайн с ограничением, что все производные второго и низших порядков равны с обеих 
сторон узла, что эквивалентно утверждению, что в узле меняется только кубическая
составляющая регрессионной функции. Ограничения на непрерывность постепенно ослабляются
серией подкоманд METHOD, которые вводят в модель XB-параметры низших порядков. Заметьте, 
что параметры низших порядков для переменной X (XA) присутствуют во всех моделях, даже если
XB-параметры низших порядков в модели отсутствуют.

COMPUTE XA2 = XA*XA.
COMPUTE XA3 = XA2*XA.
COMPUTE XB0 = (XA GE 15).
COMPUTE XB1 = XB0 * (XA - 15).
COMPUTE XB2 = XB1 * XB1.
COMPUTE XB3 = XB2 * XB1.
REGRESSION VARIABLES = y xa xa2 xa3 xb0 xb1 xb2 xb3
      /STATISTICS = DEFAULT CHA
      /DEP = y
      /METHOD ENTER XA XA2 XA3
      /METHOD ENTER XB3
      /METHOD ENTER XB1
      /METHOD ENTER XB0.

******************************************.

Решение 100009298 : Сплайновая регрессия с оценкой узлов в SPSS

Вопрос.
Как можно в SPSS устроить подгонку регрессионной модели, для
которой я ожидаю изменения коэффициента наклона в большую или
меньшую сторону (либо даже изменения знака коэффициента) при
некотором значении предиктора. При этом я хотел бы, чтобы точка, 
с которой начинается это изменение, также оценивалась бы в качестве параметра модели.
Можно ли в SPSS построить подобную модель?


Ответ.
Регрессионные модели, в которых функция меняется в одной или нескольких
точках на шкале значений предикторов, называются сплайнами, либо кусочными
полиномами, а сами точки переменения характера функции называются узловыми точками (узлами).
Если аналитику известны эти узлы, подгонка сплайнов без проблем осуществляется
в процедуре REGRESSION. Решение 100000148 описывает свойства таких сплайнов и 
процедуру их оценивания в SPSS. Если же эти узлы должны быть оценены по данным, т.е.
если узел сам представлен переменной, то множество подобных моделей можно построить в
процедурах SPSS NLR и CNLR, как это описано ниже. Эти процедуры, соответственно, 
оценивают нелинейные регрессионные модели, а также нелинейные регрессионные модели с
ограничениями.

Одно из свойств NLR, которое позволяет нам осуществлять подгонку моделей с 
переменными-узлами, это возможность вводить мультипликативные взаимосвязи между
параметрами. Команда MODEL PROGRAM для первого примера, иллюстрирует эту возможность
с 3-фазной линейной сплайновой моделью.

MODEL PROGRAM ba0=3 ba1=3 bb1=-4 bc1=2 knot1=15 knot2=25.
COMPUTE predex1 = ba0 + ba1*xa + bb1*(xa-knot1)*(xa GE knot1)
                + bc1*(xa-knot2)*(xa GE knot2).
NLR y1 with xa / pred = predex1 /save pred resid (residex1).

Параметры ba0 и ba1, соответственно, задают константу и линейный коэффициент для
предиктора XA. Параметр bb1 - это вторая фаза коррекции наклона (линейного коэффициента).
Он применяется к разности (xa - knot1) когда knot1 является первым узлом. Умножая этот 
параметр на результат логического выражения (xa GE knot1), мы ограничиваем
действие bb1 только областью значений XA, равных, или превосходящих первую узловую точку.

Параметр bc1 представляет дальнейшую коррекцию угла наклона (третью фазу). Его
взаимосвязь с точкой второго узла (knot2) аналогично взаимосвязи bb1/knot1.
Поименовав knot1 и knot2 в команде MODEL PROGRAM, мы объявили их параметрами, 
подлежащими оценке. Значения, которые мы приписали каждому параметру - только лишь
стартовые значения для итеративной процедуры их оценки методом максимального правдоподобия, 
который заложен в процедуры NLR и CNLR.

Однако, на переменные-узлы следует наложить одно ограничение, не обязательное
для моделей с фиксированными узлами. Функция должна быть непрерывна в узле. Это значит,
что в ней не может быть "скачков", или ступеней в узле. Этот ограничение, фактически,
запрещает введение константы на фазе 2 (параметра шага) при оценке процедурами NLR или CNLR.
Кроме того, примеры ниже ограничены ситуациями, когда само число узлов фиксируется
аналитиком, хотя их положение и оценивается по данным автоматически. Наконец, если в 
вашей модели предполагается существование более чем одного узла, вы можете захотеть 
наложить условие, чтобы "позднейшие" узлы имели большие значения, чем "ранние". Обычно
для этого требуется просто приписать "позднейшим" узлам более высокие стартовые значения, 
чем у "ранних", но в слабо идентифицирующихся моделях может потребоваться введение 
дополнительного ограничения. Для этого нужно использовать CNLR-процедуру с ограничениями,
указываемыми в подкоманде /BOUNDS. Такое ограничение иллюстрируется в первом примере, 
где имеется 2 узла.

Дополнительную информацию о регрессионных моделях с фиксированными и переменными узлами см.
в 9-й главе Seber, G.A.F., & Wild, C.J. (1989). Nonlinear regression. New York: Wiley.

Далее показано три примера анализа со сплайновой регрессией и переменными-узлами.
Вначале приводится INPUT-программа для создания случайного массива данных для примера.
Выполните эти команды в SPSS перед началом экспериментов с подгонкой модели.
В каждом примере предиктором является переменная XA. Каждое её значение в данных 
повторяется 3 раза. Зависимые переменные для примеров 1, 2 и 3, соответственно - Y1, Y2 и Y3.
Они генерируются командами SPSS прямо перед соответствующим анализом в примере. 
Данные для примеров достаточно рафинированные и, кроме того, мы немного
"жульничаем", устанавливая начальные значения для подбора параметров равными
тем значениям, с которыми данные были сгенерированы. Тем не менее, примеры достаточно
хорошо иллюстрируют процесс использования процедур NLR и CNLR для подгонки сплайновых моделей.

******************************************.

SET SEED = 2000000.
title 'Сплайн с оценкой узлов, хорошие стартовые значения'.
INPUT PROGRAM.
LOOP xa = 1 TO 35.
LOOP rep = 1 TO 3.
LEAVE xa.
END case.
END LOOP.
END LOOP.
END file.
END INPUT PROGRAM.
EXECUTE.

*   ПРИМЕР 1
    В первом примере линейная модель имеет 2 узла,
    стартовые значения для которых XA=15 и XA=25,
    и является непрерывной в обеих узловых точках.
* y1 - независимая переменная.
COMPUTE y1=3 + 3*xa + normal(2).
IF (xa gt 15) y1=y1 - 4*(xa-15).
IF (xa gt 25) y1=y1 + 2*(xa-25).

GRAPH /SCATTER=xa with y1.

* Ограничения на значения узлов не установлены.
MODEL PROGRAM ba0=3 ba1=3 bb1=-4 bc1=2 knot1=15 knot2=25.
COMPUTE predex1 = ba0 + ba1*xa + bb1*(xa-knot1)*(xa GE knot1)
                + bc1*(xa-knot2)*(xa GE knot2).
NLR y1 / PRED = predex1 /SAVE pred resid (residex1).
* Диаграммы для диагностики остатков и предсказанных значений.
GRAPH /SCATTER=predex1 WITH residex1.
GRAPH /SCATTER= predex1 WITH y1.
GRAPH /SCATTER=xa WITH predex1.

*  Для пользователей SPSS версий ранее 5.0, предиктор(ы) в командах NLR и CNLR
*  вам следует указывать так:.
* NLR y1 WITH xa / PRED = predex1 /SAVE pred resid (residex1).
* В позднейших версиях SPSS такая форма команды будет вызывать предупреждение, но 
* команда всё же будет выполняться.

*  Пример 1 с дополнительным ограничением: KNOT2 должен быть больше, чем KNOT1.
MODEL PROGRAM  ba0=3 ba1=3 bb1=-4 bc1=2 knot1=15 knot2=25.
COMPUTE predex1 = ba0 + ba1*xa + bb1*(xa-knot1)*(xa ge knot1)
                + bc1*(xa-knot2)*(xa ge knot2).
CNLR y1 / PRED = predex1 /SAVE pred resid (residex1)
   / BOUNDS knot2 - knot1 > 0 .
GRAPH /SCATTER=predex1 WITH residex1.
GRAPH /SCATTER= predex1 WITH y1.
GRAPH /SCATTER=xa WITH predex1.



*   ПРИМЕР 2.
* Во втором примере есть один узел в значении XA=15. Сплайн непрерывен в 
* узле, но порядок полинома в нём меняется.
* До узла функция линейна, после - квадратична.
* Y2 - независимая переменная.
COMPUTE y2 = 4 + 2*xa + normal(.5).
IF (xa gt 15) y2=y2 + (xa - 15) - .5*(xa-15)*(xa-15).
GRAPH /SCATTER=xa WITH y2.

* В этой модели мы вычислим промежуточную переменную, pc2, 
которая просто содержит результат логического выражения (xa GE knot1).
* Это сделано просто во избежание многократного повторения этого выражения.
* Параметр knot1 подлежит оценке.

MODEL PROGRAM ba0=4 ba1=2 bb1=1 bb2=-0.5 knot1=15  .
COMPUTE pc2 = (xa GE knot1).
COMPUTE predex2 = ba0 + ba1*xa + bb1*(xa-knot1)*pc2
                + bb2*((xa-knot1)**2)*pc2.
NLR y2  / PRED = predex2 /SAVE pred resid (residex2).

GRAPH /SCATTER=predex2 WITH residex2.
GRAPH /SCATTER= predex2 WITH y2.
GRAPH /SCATTER=xa WITH predex2.


*   ПРИМЕР 3.
* В третьем примере есть один узел в значении XA=15.
* Осуществляется подгонка кубического сплайна с ограничением, что
  все производные второго и низших порядков с обеих сторон от узла равны, 
  что эквивалентно утверждению, что после прохождения узловой точки меняется
  только кубическая составляющая. Затем ограничение на непрерывность постепенно
  ослабляется вводом в модель NLR параметров XB низших порядков. Если вы посмотрите
  на корреляции между оценками параметров, вы увидите, что последние 2 модели плохо
  идентифицируются, т.к. многие из этих корреляций очень высоки. В последней модели 
  некоторые из оценок параметров идеально коррелированы. Эти модели приводятся здесь
  для иллюстрации последствий снятия ограничений.
COMPUTE y3 = 6 + 3.5*xa - 1.5*xa*xa + .5*xa*xa*xa + normal(.5).
IF (xa gt 15) y3 = y3 - 2*(xa-15)*(xa-15)*(xa-15).
GRAPH /SCATTER=xa WITH y3.

MODEL PROGRAM ba0=6 ba1=3.5 ba2=-1.5 ba3=0.5 bb3=-2 knot1=15  .
COMPUTE pc2 = (xa ge knot1).
COMPUTE predex3 = ba0 + ba1*xa + ba2*xa*xa + ba3*xa*xa*xa
                + bb3*((xa-knot1)**3)*pc2.
NLR y3 / PRED = predex3 /SAVE pred resid (residex3).

GRAPH /SCATTER=predex3 WITH residex3.
GRAPH /SCATTER= predex3 WITH y3.
GRAPH /SCATTER=xa WITH predex3.

* Пример 3 с дополнительным квадратичным параметром после узла.
MODEL PROGRAM ba0=6 ba1=3.5 ba2=-1.5 ba3=0.5 bb2=0 bb3=-2 knot1=15  .
COMPUTE pc2 = (xa ge knot1).
COMPUTE pred = ba0 + ba1*xa + ba2*xa*xa + ba3*xa*xa*xa
                + bb2*((xa-knot1)**2)*pc2 + bb3*((xa-knot1)**3)*pc2.
NLR y3.

* Пример 3 с дополнительными квадратичным и линейным параметрами после узла.
MODEL PROGRAM ba0=6 ba1=3.5 ba2=-1.5 ba3=0.5
              bb1=0 bb2=0 bb3=-2 knot1=15 .
COMPUTE pc2 = (xa ge knot1).
COMPUTE pred = ba0 + ba1*xa + ba2*xa*xa + ba3*xa*xa*xa
                + bb1*(xa-knot1)*pc2 + bb2*((xa-knot1)**2)*pc2
                + bb3*((xa-knot1)**3)*pc2.
NLR y3.