Нелинейная регрессия Опции

[Pinus]

Подумал вот о чем (если вернуться к первой задаче). Если в уравнение Y = b0*(1-exp(-A/b1))^b2 + b3*Z^2 добавлять мультипликативный эффект, то наверно надо делать так: Y = b00*(1-exp(-A/b10))^b20 + b01*Z^2*(1-exp(-A/b11))^b21 + b3*Z^2 (также, как с dummy)?

Если нормальное распределение - можно и так (если уверены, что Y(Z,A)=f(A)+F(Z), а не, например f(A)*F(Z)). Но вообще-то, при любой нелинейной регрессии значимость влияния переменной Z уже весьма сложна для интерпретации. Даже в Вашем примере, когда-то (Z <0) рост Z приводит к уменьшению Y, когда-то (Z>0) рост Z приводит к росту Y. Иными словами интерпретация затруднена. В случае f(A)*f(Z) все будет еще сложнее.

Сложнее с интерпретацией или еще какие-то другие проблемы возникают?

[плав]

Подумал вот о чем (если вернуться к первой задаче). Если в уравнение Y = b0*(1-exp(-A/b1))^b2 + b3*Z^2 добавлять мультипликативный эффект, то наверно надо делать так: Y = b00*(1-exp(-A/b10))^b20 + b01*Z^2*(1-exp(-A/b11))^b21 + b3*Z^2 (также, как с dummy)?


Сложнее с интерпретацией или еще какие-то другие проблемы возникают?

Так вообщем-то думми-недумми роли не играет. Если у нас две независимых переменных А и В, то надо записывать А + В + А*В, затем вместо А и В подставляете любые формулы (причем в SAS, например, можно делать именно так, как я написал, вначале кодировать формулы простыми буквами, а затем использовать упрощенную запись).
Проблемы две - получить устойчивое решение и проинтерпретировать его. У меня при сложных нелинейных зависимостях часто появлялись проблемы с устойчивостью (меняем немного стартовые параметры и получаем другое решение), поэтому я нелинейные модели и недолюбливаю (теоретически то на модельных данных все прекрасно работает, а вот на реальных...). Может, Вам повезет больше.

[Pinus]

Если у нас две независимых переменных А и В, то надо записывать А + В + А*В, затем вместо А и В подставляете любые формулы

Просто мне казалось, что для отражения мультипликативного эффекта достаточно будет записать: Y = b0*(1-exp(-A/b1))^b2 + b3*Z^2 + b4*A*Z, а пример с dummy навел на мысль, что это неправильно, т.к. Z при взаимодействии с A, влияет на все параметры функции F(A).

Может, Вам повезет больше.

Видимо тоже не везет. С подбором параметров действительно не получается. Даже если брать только две группы качественного фактора, все равно процесс не сходится. А при большем количестве групп Q, вообще говорить не приходится. Если к одному сложному [F(A)] слагаемому добавлять две-три простых функции, то нормально, устойчивые оценки.
Что же делать, Плав? Вы писали про процедуры множественной импутации, что это такое? Может этим можно попробовать? Что еще можно здесь придумать?

Сообщение отредактировал Pinus - 23.07.2010 - 05:03

[плав]

Видимо тоже не везет. С подбором параметров действительно не получается. Даже если брать только две группы качественного фактора, все равно процесс не сходится. А при большем количестве групп Q, вообще говорить не приходится. Если к одному сложному [F(A)] слагаемому добавлять две-три простых функции, то нормально, устойчивые оценки.
Что же делать, Плав? Вы писали про процедуры множественной импутации, что это такое? Может этим можно попробовать? Что еще можно здесь придумать?

Не, множественная импутация - это методика, которая позволяет анализировать пропущенные значения. Я ее упоминал только в связи с тем, что есть подходы к суммированию оценок нескольких регрессий.
А действительно так уж нужна нелинейная функция? Это был мой вопрос с самого начала, но, по-моему DrgLena указала, что для Вас это критично. Если это действительно так, то надо попробовать поиграть со стартовыми параметрами. Иногда таким образом удается добиться схождения. Оптимально порисовать графики при разных стартовых параметрах и посмотреть какой параметр наилучшим образом описывает. Еще вариант - свалить все данные в кучу и сделать простейшую модель, потом использовать полученные параметры в качестве стартовых для всех (когда уже будете вводить групповые значения).
Вариант номер два - посмотреть на форму зависимости. У Вас она то, что называется intractable - т.е. форма, которую невозможно линеаризовать. Может быть, посмотреть функцию exp(A/B)? Тогда log(Y)=b1*(А/В)?
Еще вариант - упростить выражение, предположить, что ряд параметров одинаковы и считать исходя из этого.

[Pinus]

Вариант номер два - посмотреть на форму зависимости.

Уже тоже думал об этом. Действительно, есть линеаризуемые функции роста, попробую покрутить какую-нибудь из них (типа Гомперца). Посмотрим, что получится.
А как Вы смотрите на такую процедуру "множественной ампутации": объединить все выборки, сгруппировать значения непрерывных объясняющих переменных по интервалам (группам), сделать дисперсионный анализ (3-4 факторный) и посмотреть, что там от чего зависит, а потом уже для однородных групп строить модели со сложной функцией роста? Вдруг по ANOVA не будет никаких эффектов взаимодействия, тогда все значительно упростится.

Сообщение отредактировал Pinus - 23.07.2010 - 15:36

[плав]

Уже тоже думал об этом. Действительно, есть линеаризуемые функции роста, попробую покрутить какую-нибудь из них (типа Гомперца). Посмотрим, что получится.
А как Вы смотрите на такую процедуру "множественной ампутации": объединить все выборки, сгруппировать значения непрерывных объясняющих переменных по интервалам (группам), сделать дисперсионный анализ (3-4 факторный) и посмотреть, что там от чего зависит, а потом уже для однородных групп строить модели со сложной функцией роста? Вдруг по ANOVA не будет никаких эффектов взаимодействия, тогда все значительно упростится.

Мне еще вот какая мысль пришла в голову (может и не новая). Вам ведь надо ответить на вопрос, различаются ли показатели в разных группах?
Так вот действуем так:
1) Объединяем все наблюдения в одну группу и считает модель по Вашей формуле (общую)
2) Рассчитываем предсказанные значения и считаем квадраты разности между предсказанным значением и истинным (SSt)
3) Разбиваем данные на группы и проводим регресионный анализ внутри каждой группы
4) Для каждой группы считаем квадраты разности, затем суммируем полученные значения для всех групп (SSe)
Получили модель дисперсионного анализа и используем омнибусный тест для проверки гипотезы, что разбиение на группы не приводит к снижению дисперсии.
Единственно что надо внимательно считать число степеней свободы - каждый оцененный параметр - одна степень свободы.

Что же касается "ампутации" я не очень люблю потерю данных при превращении количественной переменной в качественную, данных не так много, чтобы информацией разбрасываться

[Pinus]

Так вот действуем так:
1) Объединяем все наблюдения в одну группу и считает модель по Вашей формуле (общую)
2) Рассчитываем предсказанные значения и считаем квадраты разности между предсказанным значением и истинным (SSt)
3) Разбиваем данные на группы и проводим регресионный анализ внутри каждой группы
4) Для каждой группы считаем квадраты разности, затем суммируем полученные значения для всех групп (SSe)
Получили модель дисперсионного анализа и используем омнибусный тест для проверки гипотезы, что разбиение на группы не приводит к снижению дисперсии.
Единственно что надо внимательно считать число степеней свободы - каждый оцененный параметр - одна степень свободы.

Как посчитать SSt и SSe понял. Что такое омнибусный тест? Посоветуйте, пожалуйста, где можно об этом почитать, чтобы разобраться и потом сослаться.

Сообщение отредактировал Pinus - 23.07.2010 - 23:31

[плав]

Как посчитать SSt и SSe понял. Что такое омнибусный тест? Посоветуйте, пожалуйста, где можно об этом почитать, чтобы разобраться и потом сослаться.

Омнибусный тест - это тот самый F-тест для сравнения дисперсий. Почитать можно в любой книге по дисперсионному анализу, только не в тех, что ориентированы на статпакеты. Тут где-то бродила ссылка на Montgomery, можно его посмотреть. Мне нравится Design an Analysis авторов Keppel и Wickens, но у меня бумажная копия, поэтому ничего не могу сказать насчет легкости приобретения - подозреваю, что не очень.

[Pinus]

Большое спасибо, Плав!
Возможные варианты понятны. Поработаю, посмотрим, что получится.

[Pinus]

Что-то не получается Nonlinear Estimation с dummy переменными.
Взять дежурный пример с прямыми линиями (исходники в файле 1).
Y - зависимая переменная, X - независимая. Качественный фактор Q с двумя уровнями. Уравнение: Y=b0+b1*X+b2*Q+b3*Q*X
Получается следующее (см. изображение):

В чем причина?

Сообщение отредактировал Pinus - 30.07.2010 - 16:52

Эскизы прикрепленных изображений

 

Прикрепленные файлы

 1.rar ( 6,46 килобайт ) Кол-во скачиваний: 443

 

[плав]

Что-то не получается Nonlinear Estimation с dummy переменными.
Взять дежурный пример с прямыми линиями (исходники в файле 1).
Y - зависимая переменная, X - независимая. Качественный фактор Q с двумя уровнями. Уравнение: Y=b0+b1*X+b2*Q+b3*Q*X
Получается следующее (см. изображение):

В чем причина?

а что там у Вас v4, v3, v2, v7?
Кроме того, не стоит брать две прямые линии, наложите хоть немного вариабельности, иначе может быть все, что угодно.

У меня, на самом деле (в SAS) - все получается вот результат GLM:
proc glm;
model y= x q x*q;
run;
=================
стандартная
Параметр оценка ошибка t-значение Pr > |t|

Свободный 12.22299828 1.27733679 9.57 <.0001
x 0.21098012 0.01789595 11.79 <.0001
q 5.16139852 1.87225244 2.76 0.0147
x*q 0.36251510 0.02760308 13.13 <.0001

А вот NLIN
proc nlin;
parms a=0.5 b=0.5 c=0.5 z=0.5;
model y=a*x+b*q+c*q*x+z;
run;
=================
Приближение
стандартной
Параметр Оценка ошибки Приближенные 95% доверительные интервалы

a 0.2110 0.0179 0.1728 0.2491
b 5.1614 1.8723 1.1708 9.1520
c 0.3625 0.0276 0.3037 0.4213
z 12.2230 1.2773 9.5004 14.9456

Видна крайняя "похожесть" результатов
Единственно, я чуть поменял данные, чтобы не было прямой линии:
data new;
input y x q;
cards;
15 20 0
20 30 0
20 41 0
25 50 0
25 60 0
25 70 0
30 80 0
30 92 0
35 100 0
35 110 0
30 20 1
35 29 1
40 40 1
45 50 1
50 60 1
55 70 1
65 79 1
70 90 1
75 100 1
;
run;

[DrgLena]

Нелинейное оценивание в Statistica будет абсолютно то же, что у плава, dummy можете создавать или не создавать, как новую перменную, а задать только вид функции.

Эскизы прикрепленных изображений

 

[Pinus]

а что там у Вас v4, v3, v2, v7?

Это обозначения переменных в таблице исх.данных (Statistica).

У меня, на самом деле (в SAS) - все получается...

Нелинейное оценивание в Statistica будет абсолютно то же, что у плава...

Да, действительно, все получается. С параболами (по другим исх.данным) тоже.

[Pinus]

Вопрос к пользователям профессиональных стат.программ (кроме Statistica). Какие еще итерационные методы получения МНК-оценок параметров нелинейной регрессии имеются в софте? В Statistica 6 есть только методы Гаусса-Ньютона и Левенберга-Маркуардта. Немного посложнее степень, в которую возводится функция, и все - оба не переваривают. Есть вообще что-нибудь еще стоящее из методов, которые можно было бы в этой ситуации попробовать?

Сообщение отредактировал Pinus - 2.08.2010 - 05:41

[плав]

Вопрос к пользователям профессиональных стат.программ (кроме Statistica). Какие еще итерационные методы получения МНК-оценок параметров нелинейной регрессии имеются в софте? В Statistica 6 есть только методы Гаусса-Ньютона и Левенберга-Маркуардта. Немного посложнее степень, в которую возводится функция, и все - оба не переваривают. Есть вообще что-нибудь еще стоящее из методов, которые можно было бы в этой ситуации попробовать?

Гаусса и наиболее крутого градиента (ЛМ - промежуточный между ними) в SAS, Mathematica, в стандартной функции nls системы R реализованы Гаусс-Ньютон, Голуб-Перейра, и функция библиотеки port. Кроме того, функция nls2 обеспечивает нахождение решения путем перебора (brute force), ML-оценки в пакете nlreg (и в PROC NLMIXED в SAS).
Но, на самом деле, я не думаю, что попытки применить другой метод что-то сильно изменят. С нелинейной регрессией надо долго играться с точки зрения стартовых параметров, проблема, скорее всего, там. Можно, конечно, попробовать brute force. Оптимально, конечно, взять хорошую книжку по регрессионному анализу (типа Дрейпера, которого я упоминал ранее) и попытаться решить задачу (упрощение) без использования пакетов - тогда станет понятно, в чем может быть проблема. Когда начнете брать производные, может станет понятно, почему методы не "берут" сложные степени (по определению нахождение максимума/минимума предполагает наличие производной - если она не определена... ).