Глупые предсказания полиномиальной регрессии

Предположим, я хочу подобрать модель линейной регрессии с полиномом второй степени (ортогональной), а затем спрогнозировать ответ. Вот коды для первой модели (m1)

x=1:100
y=-2+3*x-5*x^2+rnorm(100)
m1=lm(y~poly(x,2))
prd.1=predict(m1,newdata=data.frame(x=105:110))

Теперь давайте попробуем ту же модель, но вместо использования $ poly (x, 2) $ я буду использовать ее столбцы, например:

m2=lm(y~poly(x,2)[,1]+poly(x,2)[,2])
prd.2=predict(m2,newdata=data.frame(x=105:110))

Давайте посмотрим на краткое изложение m1 и m2.

> summary(m1)

Call:
lm(formula = y ~ poly(x, 2))

Residuals:
     Min       1Q   Median       3Q      Max 
-2.50347 -0.48752 -0.07085  0.53624  2.96516 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.677e+04  9.912e-02 -169168   <2e-16 ***
poly(x, 2)1 -1.449e+05  9.912e-01 -146195   <2e-16 ***
poly(x, 2)2 -3.726e+04  9.912e-01  -37588   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9912 on 97 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:     1 
F-statistic: 1.139e+10 on 2 and 97 DF,  p-value: < 2.2e-16 

> summary(m2)

Call:
lm(formula = y ~ poly(x, 2)[, 1] + poly(x, 2)[, 2])

Residuals:
     Min       1Q   Median       3Q      Max 
-2.50347 -0.48752 -0.07085  0.53624  2.96516 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.677e+04  9.912e-02 -169168   <2e-16 ***
poly(x, 2)[, 1] -1.449e+05  9.912e-01 -146195   <2e-16 ***
poly(x, 2)[, 2] -3.726e+04  9.912e-01  -37588   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 0.9912 on 97 degrees of freedom
Multiple R-squared:     1,      Adjusted R-squared:     1 
F-statistic: 1.139e+10 on 2 and 97 DF,  p-value: < 2.2e-16

Таким образом, m1 и m2 в основном одинаковы. Теперь посмотрим на прогнозы prd.1 и prd.2.

> prd.1
        1         2         3         4         5         6 
-54811.60 -55863.58 -56925.56 -57997.54 -59079.52 -60171.50 

> prd.2
         1          2          3          4          5          6 
  49505.92   39256.72   16812.28  -17827.42  -64662.35 -123692.53

Q1: Почему prd.2 значительно отличается от prd.1?

Q2: Как я могу получить prd.1, используя модель m2?


r regression prediction
person Stat    schedule 15.12.2012    source источник
comment
Не ответ, но достаточно высокие значения R-квадрата (0,99 с чем-то) меня всегда пугают ...   -  person Deer Hunter    schedule 16.12.2012
comment
Это вообще не проблема. Мы можем изменить $ y $ на что-то вроде $ y = -2 + 3 * x-5 * x ^ 2 + x ^ 5 + rnorm (100,15) $, и R-квадрат уменьшится до 95%, но проблема остается для предсказания.   -  person    schedule 16.12.2012
comment
Результаты первой модели где-то выглядят как плохо обусловленная матрица. Прогнозы просто следуют из бессмысленных коэффициентов, оцененных по первой модели.   -  person Deer Hunter    schedule 16.12.2012
comment
Коллинеарность огромна, как и у полиномиальных членов. Это может привести к плохой подготовке, о которой упоминал @DeerHunter. После этого вы использовали результаты для экстраполяции, усугубляя опасность. Различия в prd могут быть связаны с разным округлением в двух моделях.   -  person Peter Flom    schedule 16.12.2012
comment
@ Peter Flom, poly создает ортогональные многочлены (по умолчанию), поэтому он значительно уменьшит мультиколлинеарность (на самом деле, это и есть вся причина использования poly!). Например, в модели m2 коэффициент инфляции дисперсии (VIF) составляет около 1, поэтому мультиколлинеарность равна не проблема. (Вы можете проверить это с помощью пакета car, а затем vif (m2))   -  person    schedule 16.12.2012
comment
Различия между m1 и m2 любого значения (то есть, помимо того, как названы IV) возникают только для $assign и $terms членов. Разница в поведении зависит от внутреннего устройства класса lm и, в частности, от того, как ведет себя predict.lm. Это делает этот вопрос подходящим для SO, где сообщество R экспертов должно суметь с ним справиться.   -  person whuber    schedule 16.12.2012

Ответы (1)


arrow_upward
8
arrow_downward

m1 - правильный способ сделать это. m2 входит в целый мир боли ...

Чтобы делать прогнозы на основе m2, модель должна знать, что она приспособлена к ортогональному набору базисных функций, чтобы использовать те же базисные функции для экстраполированных новых значений данных. Сравните: poly(1:10,2)[,2] с poly(1:12,2)[,2] - первые десять значений не совпадают. Если вы точно соответствуете модели с poly(x,2), тогда predict все это понимает и поступает правильно.

Что вам нужно сделать, так это убедиться, что ваши предсказанные местоположения преобразованы с использованием того же набора базовых функций, который использовался для создания модели в первую очередь. Вы можете использовать для этого predict.poly (обратите внимание, что я называю свои объясняющие переменные x1 и x2, чтобы было легко сопоставить имена):

px = poly(x,2)
x1 = px[,1]
x2 = px[,2]

m3 = lm(y~x1+x2)

newx = 90:110
pnew = predict(px,newx) # px is the previous poly object, so this calls predict.poly

prd.3 = predict(m3, newdata=data.frame(x1=pnew[,1],x2=pnew[,2]))
person Spacedman    schedule 15.12.2012
comment
Большое спасибо за ваш ответ. Это полностью отвечает на мои вопросы. Причина, по которой я задал свой второй вопрос, заключается в следующем: предположим, мы подбираем многочлен 5-й степени, например m4 = lm (y ~ poly (x, 5)). Затем после подгонки мы хотим подобрать новую модель (m5), отбросив два члена: то есть многочлены степени 2 и 4. И, наконец, сделаем некоторые прогнозы с этой окончательной моделью (m5). Я не мог этого сделать, используя только m4. Однако это можно сделать, используя ваш аргумент и модель m3, как вы упомянули. - person Stat; 16.12.2012