Как интерпретировать коэффициенты из подгонки полиномиальной модели?

Я пытаюсь создать полином второго порядка, соответствующий некоторым имеющимся у меня данным. Допустим, я заговорю это подходит с ggplot():

ggplot(data, aes(foo, bar)) + geom_point() + 
       geom_smooth(method="lm", formula=y~poly(x, 2))

Я получил:

Таким образом, подгонка второго порядка работает довольно хорошо. Я рассчитываю это с R:

summary(lm(data$bar ~ poly(data$foo, 2)))

И я получаю:

lm(formula = data$bar ~ poly(data$foo, 2))
# ...
# Coefficients:
#                     Estimate Std. Error t value Pr(>|t|)    
# (Intercept)         3.268162   0.008282 394.623   <2e-16 ***
# poly(data$foo, 2)1 -0.122391   0.096225  -1.272    0.206
# poly(data$foo, 2)2  1.575391   0.096225  16.372   <2e-16 ***
# ....

Теперь я бы предположил, что формула для моей подгонки:

Но это просто дает мне неправильные значения. Например, если равно 3, я ожидаю, что станет примерно в 3,15. Однако, вставляя в приведенную выше формулу, я получаю: $\text{foo}$$\text{bar}$

Что дает? Я неправильно интерпретирую коэффициенты модели?

Мой подробный ответ ниже, но общий (то есть реальный) ответ на этот вопрос: 1) экспериментируйте, разбирайтесь, смотрите на данные, вы не можете сломать компьютер, что бы вы ни делали, так. , , эксперимент; или 2) RTFM .

Вот некоторый Rкод, который более или менее повторяет проблему, определенную в этом вопросе:

# This program written in response to a Cross Validated question
# http://stats.stackexchange.com/questions/95939/
# 
# It is an exploration of why the result from lm(y_x+I(x^2))
# looks so different from the result from lm(y~poly(x,2))

library(ggplot2)


epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
       geom_smooth(method = "lm", formula = y ~ poly(x, 2))

summary(lm(y~x+I(x^2)))       # Looks right
summary(lm(y ~ poly(x, 2)))   # Looks like garbage

# What happened?
# What do x and x^2 look like:
head(cbind(x,x^2))

#What does poly(x,2) look like:
head(poly(x,2))

Первый lmвозвращает ожидаемый ответ:

Call:
lm(formula = y ~ x + I(x^2))

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.92734    0.15376  25.542  < 2e-16 ***
x           -0.53929    0.11221  -4.806 5.62e-06 ***
I(x^2)       0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Вторая lmвозвращает что-то странное:

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.24489    0.02241 144.765  < 2e-16 ***
poly(x, 2)1  0.02853    0.22415   0.127    0.899    
poly(x, 2)2  1.09835    0.22415   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

Поскольку lmв двух вызовах одинаковое значение, аргументы lmдолжны быть разными. Итак, давайте посмотрим на аргументы. Очевидно, yто же самое. Это другие части. Давайте посмотрим на первые несколько наблюдений за правыми переменными в первом вызове lm. Возвращение head(cbind(x,x^2))выглядит так:

            x         
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Это как и ожидалось. Первый столбец xи второй столбец x^2. Как насчет второго вызова lm, тот, с поли? Возвращение head(poly(x,2))выглядит так:

              1         2
[1,] -0.1714816 0.2169976
[2,] -0.1680173 0.2038462
[3,] -0.1645531 0.1909632
[4,] -0.1610888 0.1783486
[5,] -0.1576245 0.1660025
[6,] -0.1541602 0.1539247

ОК, это действительно другое. Первый столбец нет x, а второй столбец нет x^2. Так что, что бы ни poly(x,2)делали, это не возвращает xи x^2. Если мы хотим знать, что polyделает, мы могли бы начать с чтения его файла справки. Так и говорим help(poly). В описании сказано:

Возвращает или оценивает ортогональные полиномы степени 1 в степени по указанному набору точек x. Все они ортогональны постоянному многочлену степени 0. В качестве альтернативы оценивают необработанные многочлены.

Теперь либо вы знаете, что такое "ортогональные полиномы", либо нет. Если вы этого не сделаете, то используйте Википедию или Bing (конечно, не Google, потому что Google злой - естественно, не такой плохой, как Apple, но все же плохой). Или вы можете решить, что вам все равно, что такое ортогональные многочлены. Вы можете заметить фразу «необработанные полиномы» и заметить немного дальше в файле справки, polyу rawкоторого есть опция, которая по умолчанию равна FALSE. Эти два соображения могут вдохновить вас попробовать, head(poly(x, 2, raw=TRUE))какой возврат:

            1        2
[1,] 1.000000 1.000000
[2,] 1.040404 1.082441
[3,] 1.080808 1.168146
[4,] 1.121212 1.257117
[5,] 1.161616 1.349352
[6,] 1.202020 1.444853

Взволнованный этим открытием (оно выглядит правильно, сейчас, да?), Вы можете попробовать summary(lm(y ~ poly(x, 2, raw=TRUE))) это. Возвращает:

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.53815 -0.13465 -0.01262  0.15369  0.61645 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)              3.92734    0.15376  25.542  < 2e-16 ***
poly(x, 2, raw = TRUE)1 -0.53929    0.11221  -4.806 5.62e-06 ***
poly(x, 2, raw = TRUE)2  0.09029    0.01843   4.900 3.84e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2241 on 97 degrees of freedom
Multiple R-squared:  0.1985,    Adjusted R-squared:  0.182 
F-statistic: 12.01 on 2 and 97 DF,  p-value: 2.181e-05

На ответ выше есть как минимум два уровня. Сначала я ответил на ваш вопрос. Во-вторых, и, что гораздо важнее, я проиллюстрировал, как вы должны сами отвечать на подобные вопросы. Каждый человек, который «знает, как программировать», прошел последовательность, подобную приведенной выше, шестьдесят миллионов раз. Даже люди, столь же удручающе плохие в программировании, как и я, постоянно проходят через эту последовательность. Это нормально, когда код не работает. Нормально неправильно понимать, что делают функции. Способ справиться с этим - обойтись, поэкспериментировать, посмотреть на данные и RTFM. Выйдите из режима «бездумного следования рецепту» и «детективного» режима.

Существует интересный подход к интерпретации полиномиальной регрессии Stimson et al. (1978) . Это включает в себя переписывание

$Y = \beta_{0} + \beta_{1} X + \beta_{2} X^{2} + u$

как

$Y = m + \beta_{2} \left( f - X \right)^{2} + u$

$m = \beta_{0} - \left. \beta_{1}^{2} \right/ 4 \beta_{2}$$\beta_{2}$$f = \left. -\beta_{1} \right/ 2 \beta_{2}$

Если вы просто хотите подтолкнуть в правильном направлении без особого суждения: poly()создайте ортогональные (не коррелированные) полиномы, в отличие от I(), что полностью игнорирует корреляцию между результирующими полиномами. Корреляция между переменными предиктора может быть проблемой в линейных моделях (см. Здесь для получения дополнительной информации о том, почему корреляция может быть проблематичной), поэтому, вероятно, лучше (в общем) использовать poly()вместо I(). Теперь, почему результаты выглядят такими разными? Ну, и poly()и I()возьмите x и преобразовать его в новый x (в случае I(), новый x просто x ^ 1 или x ^ 2, в случае poly(), новые x намного сложнее (если вы хотите знать, откуда они (а вы, вероятно, нет), вы можете начатьздесь или вышеупомянутая страница Википедии или учебник). Дело в том, что когда вы вычисляете (прогнозируете) y на основе определенного набора значений x, вам необходимо использовать преобразованные значения x, полученные либо poly()или, либо I()(в зависимости от того, какое из них было в вашей линейной модели). Так:

library(ggplot2)    

set.seed(3)
epsilon <- 0.25*rnorm(100)
x       <- seq(from=1, to=5, length.out=100)
y       <- 4 - 0.6*x + 0.1*x^2 + epsilon

# Minimum is at x=3, the expected y value there is
4 - 0.6*3 + 0.1*3^2

ggplot(data=NULL,aes(x, y)) + geom_point() + 
   geom_smooth(method = "lm", formula = y ~ poly(x, 2))

modI <- lm(y~x+I(x^2)) 
summary(modI) # Looks right
modp <- lm(y ~ poly(x, 2))
summary(modp)  # Looks like garbage

# predict y using modI
coef(modI)[1] + coef(modI)[2] * 3^1 + coef(modI)[3] * 3^2

# predict y using modp
# calculate the new x values using predict.poly()
x_poly <- stats:::predict.poly(object = poly(x,2), newdata = 3)
coef(modp)[1] + coef(modp)[2] * x_poly[1] + coef(modp)[3] * x_poly[2]

В этом случае обе модели возвращают один и тот же ответ, что предполагает, что корреляция между переменными предиктора не влияет на ваши результаты. Если бы корреляция была проблемой, оба метода предсказывали бы разные значения.

'poly' выполняет орто-нормализацию Грэма-Шмидта для полиномов 1, x, x ^ 2, ..., x ^ deg. Например, эта функция делает то же самое, что и 'poly', не возвращая атрибуты 'coef', конечно.

MyPoly <- 
function(x, deg)
{
    n <- length(x)
    ans <- NULL
    for(k in 1:deg)
    {
        v <- x^k
        cmps <- rep(0, n)
        if(k>0) for(j in 0:(k-1)) cmps <- cmps + c(v%*%ans[,j+1])*ans[,j+1]
        p <- v - cmps
        p <- p/sum(p^2)^0.5
        ans <- cbind(ans, p)
    }
    ans[,-1]
}

Я попал в эту ветку, потому что меня интересовала функциональная форма. Итак, как мы можем выразить результат 'poly' как выражение? Просто переверните процедуру Грэма-Шмидта. Вы получите беспорядок!