原始或正交多项式回归?

机器算法验证 r 回归 多项式
2022-01-16 03:43:54

我想将变量回归到我应该使用原始多项式还是正交多项式来执行此操作?我查看了处理这些问题的网站上的问题,但我真的不明白使用它们有什么区别。yx,x2,,x5

为什么我不能只做一个“正常”回归来获得 y=\sum_{i=0}^5 \beta_i x^i 的系数\βiy=i=05βixi以及 p 值和所有其他好东西),而是必须担心是否使用原始多项式或正交多项式?在我看来,这个选择超出了我想做的范围。

在我目前正在阅读的统计书中(Tibshirani 等人的 ISLR)中没有提到这些事情。事实上,他们在某种程度上被低估了。
原因是,AFAIK,lm()在 R 中的函数中,使用y ~ poly(x, 2)数量等于使用正交多项式,使用y ~ x + I(x^2)数量等于使用原始多项式。但是在第 116 页上,作者说我们使用第一个选项,因为后者是“麻烦的”,这没有表明这些命令实际上做了两个完全不同的事情(并且因此具有不同的输出)。
(第三个问题)为什么ISLR的作者会这样迷惑他们的读者?

4个回答

我相信答案不是关于数字稳定性(尽管它起作用),而是更多关于减少相关性。

从本质上讲 - 问题归结为这样一个事实,即当我们回归一堆高阶多项式时,我们回归的协变量变得高度相关。下面的示例代码:

x = rnorm(1000)
raw.poly = poly(x,6,raw=T)
orthogonal.poly = poly(x,6)
cor(raw.poly)
cor(orthogonal.poly)

这是非常重要的。随着协变量变得更加相关,我们确定哪些是重要的(以及它们的影响有多大)的能力会迅速下降。这通常被称为多重共线性问题。在极限情况下,如果我们有两个完全相关的变量,当我们对它们进行回归时,就不可能区分这两者——你可以认为这是问题的一个极端版本,但是这个问题会影响我们对相关程度也较低。因此,在真正意义上——即使数值不稳定性不是问题——来自高阶多项式的相关性对我们的推理程序造成了巨大的损害。这将表现为您会看到的更大的标准误差(因此更小的 t-stats)(参见下面的示例回归)。

y = x*2 + 5*x**3 - 3*x**2 + rnorm(1000)
raw.mod = lm(y~poly(x,6,raw=T))
orthogonal.mod = lm(y~poly(x,6))
summary(raw.mod)
summary(orthogonal.mod)

如果你运行这段代码,解释有点困难,因为系数都在变化,所以很难比较。不过,查看 T 统计数据,我们可以看到使用正交多项式确定系数的能力要大得多。对于 3 个相关系数,正交模型的 t-stats 为 (560,21,449),而原始多项式模型的 t-stats 为 (28,-38,121)。对于只有几个相对低阶多项式项很重要的简单模型来说,这是一个巨大的差异。

这并不是说这是没有成本的。有两个主要成本需要牢记。1)我们失去了正交多项式的一些可解释性。我们可能理解系数 on 的x**3含义,但解释系数 on x**3-3x(第三个 Hermite poly——不一定是你将使用的)可能要困难得多。其次——当我们说这些多项式是正交的——我们的意思是它们相对于某种距离度量是正交的。选择与您的情况相关的距离度量可能很困难。然而,话虽如此,我相信该poly函数的设计是为了选择它与协方差正交——这对于线性回归很有用。

我觉得其中几个答案都没有抓住重点。海涛的回答解决了拟合原始多项式的计算问题,但很明显,OP 是在询问两种方法之间的统计差异。也就是说,如果我们有一台可以准确表示所有值的完美计算机,为什么我们更喜欢一种方法而不是另一种呢?

user5957401 认为正交多项式降低了多项式函数之间的共线性,这使得它们的估计更加稳定。我同意 Jake Westfall 的批评;正交多项式中的系数表示与原始多项式系数完全不同的量。无论您使用正交多项式还是原始多项式,模型隐含的剂量反应函数的系数具有“的瞬时变化”的解释的正交多项式执行了边际效应程序R2XYX=0X=0,即使正交多项式回归中一阶项的系数和标准误差与其在原始多项式回归中的值完全不同,您也会得到完全相同的斜率和标准误差。也就是说,当试图从两个回归中获得相同的数量(即可以以相同方式解释的数量)时,估计值和标准误差将是相同的。使用正交多项式并不意味着您神奇地对在任何给定点的斜率有更多确定性。模型的稳定性是相同的。见下文:X

data("iris")

#Raw:
fit.raw <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
summary(fit.raw)

#> Coefficients:
#>                  Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)        1.1034     0.1304   8.464 2.50e-14 ***
#> Petal.Width        1.1527     0.5836   1.975  0.05013 .  
#> I(Petal.Width^2)   1.7100     0.5487   3.116  0.00221 ** 
#> I(Petal.Width^3)  -0.5788     0.1408  -4.110 6.57e-05 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Residual standard error: 0.3898 on 146 degrees of freedom
#> Multiple R-squared:  0.9522, Adjusted R-squared:  0.9512 
#> F-statistic: 969.9 on 3 and 146 DF,  p-value: < 2.2e-16

#Orthogonal
fit.orth <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), data = iris)

#Marginal effect of X at X=0 from orthogonal model
library(margins)
summary(margins(fit.orth, variables = "Petal.Width", 
                at = data.frame(Petal.Width = 0)))
#> Warning in check_values(data, at): A 'at' value for 'Petal.Width' is
#> outside observed data range (0.1,2.5)!
#>       factor Petal.Width    AME     SE      z      p  lower  upper
#>  Petal.Width      0.0000 1.1527 0.5836 1.9752 0.0482 0.0089 2.2965

reprex 包(v0.3.0)于 2019 年 10 月 25 日创建

正交拟合的 0 处的边际效应Petal.Width及其标准误差与原始多项式拟合的边际效应完全相同(即1.1527。使用正交多项式不会提高两个模型之间相同数量的估计精度。

关键如下:使用正交多项式允许您隔离每个项对解释结果方差的贡献,例如,通过平方半偏相关来衡量。如果您拟合 3 阶正交多项式,则每个项的平方半偏相关表示模型中该项解释的结果的方差。的方差有多少的线性分量YX?” 您可以拟合正交多项式回归,线性项上的平方半偏相关将表示该数量。原始多项式并非如此。如果您拟合相同阶的原始多项式模型,则线性项不代表的线性分量解释的中的方差比例。见下文。YX

library(jtools)
data("iris")

fit.raw3 <- lm(Petal.Length ~ Petal.Width + I(Petal.Width^2) +
                  I(Petal.Width^3), data = iris)
fit.raw1 <- lm(Petal.Length ~ Petal.Width, data = iris)

round(summ(fit.raw3, part.corr = T)$coef, 3)
#>                    Est.  S.E. t val.     p partial.r part.r
#> (Intercept)       1.103 0.130  8.464 0.000        NA     NA
#> Petal.Width       1.153 0.584  1.975 0.050     0.161  0.036
#> I(Petal.Width^2)  1.710 0.549  3.116 0.002     0.250  0.056
#> I(Petal.Width^3) -0.579 0.141 -4.110 0.000    -0.322 -0.074

round(summ(fit.raw1, part.corr = T)$coef, 3)
#>              Est.  S.E. t val. p partial.r part.r
#> (Intercept) 1.084 0.073 14.850 0        NA     NA
#> Petal.Width 2.230 0.051 43.387 0     0.963  0.963

fit.orth3 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3), 
               data = iris)
fit.orth1 <- lm(Petal.Length ~ stats::poly(Petal.Width, 3)[,1], 
               data = iris)

round(summ(fit.orth3, part.corr = T)$coef, 3)
#>                                Est.  S.E.  t val. p partial.r part.r
#> (Intercept)                   3.758 0.032 118.071 0        NA     NA
#> stats::poly(Petal.Width, 3)1 20.748 0.390  53.225 0     0.975  0.963
#> stats::poly(Petal.Width, 3)2 -3.015 0.390  -7.735 0    -0.539 -0.140
#> stats::poly(Petal.Width, 3)3 -1.602 0.390  -4.110 0    -0.322 -0.074

round(summ(fit.orth1, part.corr = T)$coef, 3)
#>                                    Est.  S.E. t val. p partial.r part.r
#> (Intercept)                       3.758 0.039 96.247 0        NA     NA
#> stats::poly(Petal.Width, 3)[, 1] 20.748 0.478 43.387 0     0.963  0.963

reprex 包(v0.3.0)于 2019 年 10 月 25 日创建

当 3 阶多项式拟合时,原始多项式的平方半偏相关为当仅拟合线性项时,平方半偏相关为拟合 3 阶多项式时,正交多项式的平方半偏相关为当只拟合线性项时,平方半偏相关仍为0.0010.0030.0050.9270.9270.0200.0050.927. 从正交多项式模型而不是原始多项式模型中,我们知道结果中解释的大部分方差是由于线性项,很少来自平方项,甚至更少来自三次项。原始多项式值并不能说明这个问题。

现在,如果您希望这种解释收益超过实际能够理解模型系数的解释收益,那么您应该使用正交多项式。如果您希望查看系数并确切知道它们的含义(尽管我怀疑通常会这样做),那么您应该使用原始多项式。如果您不在乎(即,您只想控制混淆或生成预测值),那么它真的没关系;两种表格都包含有关这些目标的相同信息。我还认为正交多项式在正则化(例如套索)中应该是首选,因为删除高阶项不会影响低阶项的系数,这对于原始多项式是不正确的,

为什么我不能只做一个“正常”回归来获得系数?

因为它在数值上不稳定。请记住,计算机使用固定位数来表示浮点数。查看IEEE754了解详细信息,您可能会感到惊讶,即使是简单的数字,计算机也需要将其存储为您可以在这里尝试其他号码0.40.4000000059604644775390625

使用原始多项式会导致问题,因为我们会有大量的数字。这是一个小证明:我们将矩阵条件数与原始和正交多项式进行比较。

> kappa(model.matrix(mpg~poly(wt,10),mtcars))
[1] 5.575962
> kappa(model.matrix(mpg~poly(wt,10, raw = T),mtcars))
[1] 2.119183e+13

您也可以在此处查看我的答案作为示例。

为什么高阶多项式的系数很大

我本来只是评论提到这一点,但我没有足够的代表,所以我会尝试扩展为一个答案。您可能有兴趣在“统计学习简介”(James 等人,2017 年,更正了第 8 次印刷)的实验室第 7.8.1 节中看到,他们确实讨论了是否使用正交多项式之间的一些差异,即使用raw=TRUEraw=FALSEpoly()函数中。例如,系数估计值会改变,但拟合值不会:

# using the Wage dataset in the ISLR library
fit1 <- lm(wage ~ poly(age, 4, raw=FALSE), data=Wage)
fit2 <- lm(wage ~ poly(age, 4, raw=TRUE), data=Wage)
print(coef(fit1)) # coefficient estimates differ
print(coef(fit2))
all.equal(predict(fit1), predict(fit2)) #returns TRUE    

本书还讨论了在使用正交多项式时,使用anova()嵌套 F 检验获得的 p 值(以探索可能保证多项式的次数)与使用标准 t 检验时获得的 p 值相同,输出为summary(fit). 这说明 F 统计量在某些情况下等于 t 统计量的平方。