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

36

Я пытаюсь создать полином второго порядка, соответствующий некоторым имеющимся у меня данным. Допустим, я заговорю это подходит с 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,268-0,122Foo+1,575Foo2

Но это просто дает мне неправильные значения. Например, если равно 3, я ожидаю, что станет примерно в 3,15. Однако, вставляя в приведенную выше формулу, я получаю: барFooбар

барзнак равно3,268-0,1223+1,57532знак равно17,077

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

user13907
источник
2
Ответ на
6
@whuber Если бы я знал, что проблема была в «ортогональных многочленах», я, вероятно, нашел бы ответ. Но если вы не знаете, что искать, это немного сложно.
user13907
2
Вы также можете найти ответы, выполнив поиск по poly , что заметно в вашем коде. Я помещаю такую ​​информацию в комментариях по двум причинам: (1) ссылки могут помочь будущим читателям, а также вам, и (2) они могут помочь показать вам, как использовать нашу (несколько своеобразную) поисковую систему.
whuber
7
Вы разместили вопрос, касающийся вашего использования polyбез ввода ?polyR? Надпись « Вычислять ортогональные многочлены » вверху крупными дружескими буквами.
Glen_b
4
@Glen_b Да, хорошо, я сделал вид в ?polyпонимать синтаксис. По общему признанию, я только немного знаю о понятиях позади этого. Я не знал, что было что-то еще (или такая большая разница между «нормальными» многочленами и ортогональными многочленами), и примеры, которые я видел онлайн, все использовались poly()для подгонки, особенно с ggplot- так почему бы мне просто не использовать это и быть смущенным, если результат был "неправильным"? Имейте в виду, я не разбираюсь в математике - я просто применяю то, что, как я видел, делают другие, и пытаюсь понять это.
user13907

Ответы:

55

Мой подробный ответ ниже, но общий (то есть реальный) ответ на этот вопрос: 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. Выйдите из режима «бездумного следования рецепту» и «детективного» режима.

Билл
источник
7
Я думаю, что это заслуживает +6. Я постараюсь вспомнить через пару дней, когда это станет возможным. FTR, я думаю, что это не должно быть таким саркастичным, но он хорошо показывает, что такое ортогональные полиномы / как они работают, и показывает процесс, который вы используете, чтобы понять такие вещи.
gung - Восстановить Монику
13
Отличный ответ, спасибо. Хотя я немного обижен «RTFM» (но, может быть, это только я): проблема в том, что во всем, что я читал, по крайней мере, в отношении линейной регрессии в R, люди иногда делают это, другие делают это. Честно говоря, я не понимаю статью в Википедии об ортогональных многочленах. Мне не приходит в голову, почему можно использовать это для регрессии, если получаемые вами коэффициенты «неправильны». Я не математик - я стараюсь следовать рецептам, потому что я не ученый повар, но, тем не менее, мне нужно что-нибудь съесть.
user13907
12
@ user13907, это не только ты. Это действительно хороший ответ, который заслуживает того, чтобы за него проголосовали, но было бы полезно иметь более приятный тон.
Вальдир Леонсио
8
Вам не нужно понимать, что такое ортогональные полиномы - вам просто нужно понять, что это не то, что вам нужно. Почему кто-то может хотеть ортогональные полиномы? Передайте cov (poly (x, 2)), чтобы найти, что ковариация между двумя членами в полиноме равна нулю (с точностью до ошибки округления). Это ключевое свойство ортогональных многочленов - их члены имеют нулевую ковариацию друг с другом. Иногда для ваших переменных RHS удобно иметь нулевую корреляцию друг с другом. Их коэффициенты не являются неправильными, на самом деле, они просто должны интерпретироваться по-разному.
Билл
2
О, хорошо, это объяснение на простом английском языке теперь имеет смысл. Спасибо.
user13907
5

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

Yзнак равноβ0+β1Икс+β2Икс2+U

как

Yзнак равном+β2(е-Икс)2+U

мзнак равноβ0-β12/4β2β2езнак равно-β1/2β2

Дерден
источник
2
+1 Для соответствующих анализов, пожалуйста, смотрите stats.stackexchange.com/questions/28730 и stats.stackexchange.com/questions/157629 .
whuber
4

Если вы просто хотите подтолкнуть в правильном направлении без особого суждения: 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]

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

filups21
источник
1

'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' как выражение? Просто переверните процедуру Грэма-Шмидта. Вы получите беспорядок!

izmirlig
источник