R函数`poly`真的有用吗？

Question

我已阅读了手册页?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真正做的是什么？

Answer 1

要获得普通多项式,请使用问题raw = TRUE.不幸的是,回归中的普通多项式存在不希望的方面.如果我们符合二次,再说了,然后立方立方的较低阶系数均高于二次不同,即56.900099702,-0.466189630,0.001230536的二次主场迎战6.068478e + 01,-5.688501e-01, 2.079011e-03用以下立方体改装后.

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

Answer 2

当在统计模型中引入多项式项时,通常的动机是确定响应是否是"弯曲的"以及当加入该项时曲率是否"显着".投入一个+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为中心,线性项已向下移动,因此在相同的x点处为0((Intercept)模型中的隐含项依赖于在请求预测时将所有内容移回.)

Answer 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