Что делает функция "pol" на самом деле?

Я прочитал страницу руководства ?poly (которую, я признаю, я не полностью понял), а также прочитал описание функции в книге Введение в статистическое обучение.

Насколько я понимаю, вызов poly(horsepower, 2) должен быть эквивалентен написанию horsepower + I(horsepower^2). Однако это, похоже, противоречит выводу следующего кода:

library(ISLR)

summary(lm(mpg~poly(horsepower,2), data=Auto))$coef

    #                       Estimate Std. Error   t value      Pr(>|t|)
    #(Intercept)            23.44592  0.2209163 106.13030 2.752212e-289
    #poly(horsepower, 2)1 -120.13774  4.3739206 -27.46683  4.169400e-93
    #poly(horsepower, 2)2   44.08953  4.3739206  10.08009  2.196340e-21

summary(lm(mpg~horsepower+I(horsepower^2), data=Auto))$coef

    #                    Estimate   Std. Error   t value      Pr(>|t|)
    #(Intercept)     56.900099702 1.8004268063  31.60367 1.740911e-109
    #horsepower      -0.466189630 0.0311246171 -14.97816  2.289429e-40
    #I(horsepower^2)  0.001230536 0.0001220759  10.08009  2.196340e-21

Мой вопрос: почему не совпадают выходы, и что на самом деле делает poly?

Ответ 1

При введении полиномиальных членов в статистической модели обычная мотивация заключается в том, чтобы определить, является ли ответ "изогнутым", и является ли кривизна "значимой", когда этот термин добавлен. Результат метания в терминах +I(x^2) что незначительные отклонения могут быть "увеличены" процессом подгонки в зависимости от их местоположения и неверно истолкованы из-за термина кривизны, когда они были просто колебаниями на одном конце или другом диапазоне данных. Это приводит к неуместному присваиванию деклараций "значимости".

Если вы просто вставляете квадрат с I(x^2), по необходимости он также будет сильно коррелирован с x, по крайней мере, в домене, где x > 0. Вместо этого: poly(x,2) создает "изогнутый" набор переменных, где линейный член не так сильно коррелирован с x, и где кривизна примерно одинакова по всему диапазону данных. (Если вы хотите прочитать статистическую теорию, найдите "ортогональные полиномы".) Просто введите poly(1:10, 2) и посмотрите на два столбца.

poly(1:10, 2)
                1           2
 [1,] -0.49543369  0.52223297
 [2,] -0.38533732  0.17407766
 [3,] -0.27524094 -0.08703883
 [4,] -0.16514456 -0.26111648
 [5,] -0.05504819 -0.34815531
 [6,]  0.05504819 -0.34815531
 [7,]  0.16514456 -0.26111648
 [8,]  0.27524094 -0.08703883
 [9,]  0.38533732  0.17407766
[10,]  0.49543369  0.52223297
attr(,"degree")
[1] 1 2
attr(,"coefs")
attr(,"coefs")$alpha
[1] 5.5 5.5

attr(,"coefs")$norm2
[1]   1.0  10.0  82.5 528.0

attr(,"class")
[1] "poly"   "matrix"

"Квадратичный" член центрируется на 5.5, а линейный член сдвигается вниз, так что он равен 0 в той же точке x (при неявном члене (Intercept) в модели зависит от смещения всего назад на запрашиваются временные предсказания.)

Ответ 2

Необработанные полиномы

Чтобы получить обычные полиномы, как в вопросе, используйте raw = TRUE. К сожалению, в регрессии есть нежелательный аспект с обычными полиномами. Если мы подгоняем квадратик, скажем, а затем кубику, то все коэффициенты младшего порядка кубики будут отличаться от квадратичной, т.е..688501e -0 1, 2.079011e -0 3 после переоборудования с кубиком ниже.

library(ISLR)
fm2raw <- lm(mpg ~ poly(horsepower, 2, raw = TRUE), Auto)
cbind(coef(fm2raw))
##                                          [,1]
## (Intercept)                      56.900099702
## poly(horsepower, 2, raw = TRUE)1 -0.466189630
## poly(horsepower, 2, raw = TRUE)2  0.001230536

fm3raw <- lm(mpg ~ poly(horsepower, 3, raw = TRUE), Auto)
cbind(coef(fm3raw))
##                                           [,1]
## (Intercept)                       6.068478e+01
## poly(horsepower, 3, raw = TRUE)1 -5.688501e-01
## poly(horsepower, 3, raw = TRUE)2  2.079011e-03
## poly(horsepower, 3, raw = TRUE)3 -2.146626e-06

Ортогональные многочлены

Что нам действительно нужно, так это добавить кубический член таким образом, чтобы коэффициенты более низкого порядка, которые были подобраны с использованием квадратичного, остались прежними после переоборудования с кубическим подгонкой. Для этого возьмите линейные комбинации столбцов в poly(horsepower, 2, raw = TRUE) и сделайте то же самое с poly(horsepower, 3, raw = TRUE), чтобы столбцы в квадратичной аппроксимации были ортогональны друг другу и аналогично для кубической аппроксимации. Этого достаточно, чтобы гарантировать, что коэффициенты более низкого порядка не изменятся, когда мы добавим коэффициенты более высокого порядка. Обратите внимание, что первые три коэффициента теперь одинаковы в двух наборах ниже (тогда как выше они отличаются). То есть в обоих случаях ниже 3 коэффициентов более низкого порядка 23,44592, -120.13774 и 44,08953.

fm2 <- lm(mpg ~ poly(horsepower, 2), Auto)
cbind(coef(fm2))
##                            [,1]
## (Intercept)            23.44592
## poly(horsepower, 2)1 -120.13774
## poly(horsepower, 2)2   44.08953

fm3 <- lm(mpg ~ poly(horsepower, 3), Auto)
cbind(coef(fm3))
##                             [,1]
## (Intercept)            23.445918
## poly(horsepower, 3)1 -120.137744
## poly(horsepower, 3)2   44.089528
## poly(horsepower, 3)3   -3.948849

Те же прогнозы

Важно отметить, что поскольку столбцы в poly(horsepwer, 2) являются просто линейными комбинациями столбцов в poly(horsepower, 2, raw = TRUE), две квадратичные модели (ортогональные и необработанные) представляют одинаковые модели (т.е. они дают одинаковые прогнозы) и отличаются только параметризацией. Например, установленные значения одинаковы:

all.equal(fitted(fm2), fitted(fm2raw))
## [1] TRUE

Это также относится к необработанным и ортогональным кубическим моделям.

Ортогональность

Мы можем проверить, что у многочленов есть ортогональные столбцы, которые также ортогональны пересечению:

nr <- nrow(Auto)
e <- rep(1, nr) / sqrt(nr) # constant vector of unit length
p <- cbind(e, poly(Auto$horsepower, 2))
zapsmall(crossprod(p))
##   e 1 2
## e 1 0 0
## 1 0 1 0
## 2 0 0 1

Остаточная сумма квадратов

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

crossprod(model.matrix(fm3)[, 4], Auto$mpg)^2
##         [,1]
## [1,] 15.5934

deviance(fm2) - deviance(fm3) # same
## [1] 15.5934

anova(fm2, fm3) # same value found on last line
## 
## Analysis of Variance Table
## 
## Model 1: mpg ~ poly(horsepower, 2)
## Model 2: mpg ~ poly(horsepower, 3)
##   Res.Df    RSS Df Sum of Sq      F Pr(>F)
## 1    389 7442.0                           
## 2    388 7426.4  1    15.593 0.8147 0.3673

Ответ 3

Быстрый ответ заключается в том, что poly вектора x по существу эквивалентен QR-разложению матрицы, столбцы которой имеют степень x (после центрирования). Например:

> x<-rnorm(50)
> x0<-sapply(1:5,function(z) x^z)
> x0<-apply(x0,2,function(z) z-mean(z))
> x0<-qr.Q(qr(x0))
> cor(x0,poly(x,5))
                 1             2             3             4             5
[1,] -1.000000e+00 -1.113975e-16 -3.666033e-17  7.605615e-17 -1.395624e-17
[2,] -3.812474e-17  1.000000e+00  1.173755e-16 -1.262333e-17 -3.988085e-17
[3,] -7.543077e-17 -7.778452e-17  1.000000e+00  3.104693e-16 -8.472204e-17
[4,]  1.722929e-17 -1.952572e-16  1.013803e-16 -1.000000e+00 -1.611815e-16
[5,] -5.973583e-17 -1.623762e-18  9.163891e-17 -3.037121e-16  1.000000e+00