r - R 函数 `poly` 到底有什么作用？

Question

我已经通读了手册页?poly(我承认我没有完全理解)，并且还阅读了书Introduction to Statistical Learning中该函数的描述。 。

我目前的理解是，调用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 retrofit 后的立方体如下。

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 生成的矩阵的列具有单位长度并且相互正交(并且也与截距列正交)，因此残差和的减少由于添加三次项而产生的平方数只是响应向量在模型矩阵的三次列上的投影长度的平方。

# these three give the same result

# 1. squared length of projection of y, i.e. Auto$mpg, on cubic term column
crossprod(model.matrix(fm3)[, 4], Auto$mpg)^2
##         [,1]
## [1,] 15.5934

# 2. difference in sums of squares
deviance(fm2) - deviance(fm3) 
## [1] 15.5934

# 3. difference in sums of squares from anova
anova(fm2, fm3) 
## 
## 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  <-- note Sum of Sq value