R 函数 `poly` 的真正作用是什么？答案

【问题标题】：What does the R function `poly` really do?R 函数 `poly` 的真正作用是什么？
【发布时间】：2020-04-14 20:59:43
【问题描述】：

我已经阅读了手册页?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 到底在做什么？

【问题讨论】：

看看这个问题的答案：mathoverflow.net/questions/38864/…
我只是看了一下（不完全明白），但我还是想知道，在这种情况下，poly(horsepower,2) 生成的封闭式公式究竟是什么？
试试poly(horsepower, degree=2, raw=TRUE);您将 2 作为错误参数传递，raw 默认为 FALSE。
@baptiste，您的建议可以使poly 生成与显式公式相同的输出，但我仍然想知道poly 生成的“正交多项式”的实际形式没有那个参数。另外，根据手册，我将 2 作为度数传递：“虽然应该正式命名 '度'（因为它遵循 '...'），但长度为 1 的未命名的第二个参数将被解释为度数。”
好点 re...，在它之后命名参数仍然是个好习惯。

标签： r

【解决方案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

残差平方和

正交多项式的另一个不错的特性是，由于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

【讨论】：

您能解释一下为什么e 被计算为除以nr 的平方根吗？
所以它的长度sqrt(crossprod(e, e))是1。
您说“poly(horsepower, 2) 的列只是poly(horsepower, 2, raw = TRUE) 的列的线性组合”。是否可以准确提取线性组合是什么？这样我就可以，例如，从非正交多项式中手动推导出正交多项式，或者写一个方程？
@Gimelist, (1) 给定h <- Auto$horsepower; p <- poly(h, 2); p.raw <- poly(h, 2, raw = TRUE); co <- coef(lm(p ~ p.raw))，我们有co 是所需的转换矩阵，即cbind(1, p.raw) %*% co 等于p。 (2) 同样变换矩阵是cbind(1, p.raw)的格拉姆施密特分解的R矩阵的逆矩阵，即library(pracma); gs <- gramSchmidt(cbind(1, p.raw)); cbind(1, p.raw) %*% solve(gs$R)[, -1]等于p
您是否可以为这个问题添加答案：这些多项式与哪个内积正交？这就是我无法绕开我的大脑。我习惯于将 Hermite、Laguerre、Legendre 等多项式视为使用 Gram-Schmidt 相对于特定内积构造的多项式。谢谢！

【解决方案2】：

在统计模型中引入多项式项时，通常的动机是确定响应是否“弯曲”以及添加该项时曲率是否“显着”。添加+I(x^2) 项的结果是微小的偏差可能会根据它们的位置被拟合过程“放大”，并且当它们只是数据范围的一端或另一端的波动时，会被误解为由于曲率项。这会导致不恰当地分配“重要性”声明。

如果您只是将一个平方项与I(x^2) 相加，那么至少在x > 0 所在的域中，它也必然与x 高度相关。改为使用：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) 项依赖于将所有内容移回需要时间预测。）

【讨论】：

这是我最喜欢的话题之一，我很高兴看到@G.Grothendieck 的回答，因为我钦佩他的知识深度。他的答案优于我的答案，因为他使用提供的数据集作为问题中的“挑战”。

【解决方案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

【讨论】：

这在数学基础上可能比我接受的答案“更深”。