lm() 中的 poly()：原始与正交之间的区别答案

【问题标题】：poly() in lm(): difference between raw vs. orthogonallm() 中的 poly()：原始与正交之间的区别
【发布时间】：2017-07-21 09:50:53
【问题描述】：

我有

library(ISLR)
attach(Wage)

# Polynomial Regression and Step Functions

fit=lm(wage~poly(age,4),data=Wage)
coef(summary(fit))

fit2=lm(wage~poly(age,4,raw=T),data=Wage)
coef(summary(fit2))

plot(age, wage)
lines(20:350, predict(fit, newdata = data.frame(age=20:350)), lwd=3, col="darkred")
lines(20:350, predict(fit2, newdata = data.frame(age=20:350)), lwd=3, col="darkred")

预测线似乎相同，但为什么系数如此不同？你如何在raw=T 和raw=F 中解释它们。

我看到poly(...,raw=T) 产生的系数与~age+I(age^2)+I(age^3)+I(age^4) 产生的系数相匹配。

如果我想使用系数“手动”获得预测（不使用predict() 函数），有什么需要注意的吗？我应该如何解释poly()中正交多项式的系数。

【问题讨论】：

你读过这个：inside-r.org/r-doc/stats/poly 吗？ raw 表示是否使用正交多项式。
是的，我做到了。可能需要一个关于模型构建中正交多项式的小解释
我认为这更多的是关于 R 如何在内部存储模型，并且在预测时几乎没有什么区别，只关心数值精度，但我不是统计学家，这样的问题可能会更好适合数学堆栈交换。
它会显着影响系数，所以我认为它的内部影响不大
en.wikipedia.org/wiki/Orthogonal_polynomials Wikipedia 有一个关于它的页面。这有点像您正在更改轴系并且您的坐标发生了巨大变化，但矢量将保持不变。（可能不是最好的类比，但我又不是数学家。）

标签： r

【解决方案1】：

默认情况下，raw = FALSE、poly() 计算正交多项式。它首先在内部使用原始编码 x、x^2、x^3、... 设置模型矩阵，然后缩放列，使每一列与前一列正交。这不会改变拟合值，但好处是您可以看到多项式中的某个阶是否显着改善了低阶回归。

考虑简单的cars 数据，响应停止distance 并驱动speed。从物理上讲，这应该具有二次关系，但在这个（旧）数据集中，二次项并不重要：

m1 <- lm(dist ~ poly(speed, 2), data = cars)
m2 <- lm(dist ~ poly(speed, 2, raw = TRUE), data = cars)

在正交编码中，您在summary(m1) 中得到以下系数：

                Estimate Std. Error t value Pr(>|t|)    
(Intercept)       42.980      2.146  20.026  < 2e-16 ***
poly(speed, 2)1  145.552     15.176   9.591 1.21e-12 ***
poly(speed, 2)2   22.996     15.176   1.515    0.136

这表明存在高度显着的线性效应，而二阶不显着。后一个 p 值（即多项式中的最高阶）与原始编码中的相同：

                            Estimate Std. Error t value Pr(>|t|)
(Intercept)                  2.47014   14.81716   0.167    0.868
poly(speed, 2, raw = TRUE)1  0.91329    2.03422   0.449    0.656
poly(speed, 2, raw = TRUE)2  0.09996    0.06597   1.515    0.136

但低阶 p 值会发生显着变化。原因是在模型m1 中，回归量是正交的，而在m2 中它们是高度相关的：

cor(model.matrix(m1)[, 2], model.matrix(m1)[, 3])
## [1] 4.686464e-17
cor(model.matrix(m2)[, 2], model.matrix(m2)[, 3])
## [1] 0.9794765

因此，在原始编码中，如果 speed^2 保留在模型中，您只能解释 speed 的 p 值。并且由于两个回归变量高度相关，因此可以删除其中一个。但是，在正交编码中speed^2 只捕获了未被线性项捕获的二次部分。然后很明显，线性部分很重要，而二次部分没有额外的意义。

【讨论】：

很棒的答案。非常感谢。只是一个小问题：应该如何解释正交多项式的系数？
我认为确切的单位很难解释，但我认为这在任一多项式中都是正确的。但在正交情况下，二次项只给出与线性多项式的偏差；三次项与二次多项式等的偏差。
@ECII @Achim 在对stackoverflow.com/questions/31457230/… 的回答中，我给出了一个将正交多项式转换为“常规”幂系数的函数。尽管它似乎有效，但我不喜欢该功能，但它确实很丑：-0，我会很感激有关替代方法的建议。
@AchimZeileis 因为我会运行 dist~speed+I(speed^2)，这相当于你的示例中的 poly(speed,2,raw=T)，所以我会得出结论根据 0.656 或 0.136 的 p 值，速度和速度^2 似乎都具有统计学意义。另一方面，模型 m1 显示速度具有统计显着性，但不是 speed^2。我很难理解这...两条相互矛盾的信息。你能帮我更好地理解这个概念吗？
@Rahul 这就是正交化的重点。在原始编码中，您只能解释 speed^2 的速度 p 值保留在模型中。并且由于两个回归变量高度相关，其中一个可以下降。但是，在正交编码中，speed^2只捕捉到了线性项没有捕捉到的二次部分。然后很明显，线性部分很重要，而二次部分没有额外的意义。

【解决方案2】：

我相信基于raw=T 运行多项式回归的方式是查看最高幂项并根据pvalue 对该系数的重要性进行评估。

如果发现不显着（大pvalue），那么回归将在没有特定的非显着幂（即下一个较低的程度）的情况下重新运行，如果不是，这将一次执行一步减少重大。

如果在任何时候更高的程度很重要，那么该过程将停止并断言该程度是适当的。

【讨论】：