R 中的高（或非常高）阶多项式回归（或替代方案？）

Question

我想对 R 中的一组数据进行（非常）高阶回归，但是 poly() 函数的阶数限制为 25。

对于此应用程序，我需要一个范围为 100 到 120 的订单。

model <- lm(noisy.y ~ poly(q,50))
# Error in poly(q, 50) : 'degree' must be less than number of unique points
model <- lm(noisy.y ~ poly(q,30))
# Error in poly(q, 30) : 'degree' must be less than number of unique points
model <- lm(noisy.y ~ poly(q,25))
# OK

Answer 1

多项式和正交多项式

poly(x) 对 degree 没有硬编码限制。然而，在实践中有两个数值约束。

基函数是在 x 值的唯一位置上构建的。 k 次多项式有 k + 1 基和系数。 poly 生成没有截距项的基础，因此 degree = k 意味着 k 基础和 k 系数。如果有 n 个唯一的 x 值，则必须满足 k <= n，否则根本没有足够的信息来构造多项式。在 poly() 内，以下行检查此条件：
```
if (degree >= length(unique(x))) 
    stop("'degree' must be less than number of unique points")
```
随着 k 的增加，x ^ k 和 x ^ (k+1) 之间的相关性越来越接近 1。这种接近速度当然取决于x值。 poly先生成普通多项式基，然后进行QR分解求正交跨度。如果 x ^ k 和 x ^ (k+1) 之间出现数值排名不足，poly 也会停止并抱怨：
```
if (QR$rank < degree) 
    stop("'degree' must be less than number of unique points")
```
但在这种情况下，错误消息并不提供信息。此外，这不一定是错误；这可能是一个警告，然后 poly 可以将 degree 重置为 rank 以继续。也许 R 核心可以改进这一点？？

您的反复试验表明您无法构造超过 25 次的多项式。可以先查看length(unique(q))。如果你的学位比这个小但仍然触发错误，你肯定知道这是由于数字排名不足。

但是我想说的是3-5次以上的多项式是没有用的！关键的原因是Runge's phenomenon。在统计术语方面：高阶多项式总是严重过度拟合数据！。不要天真地认为因为正交多项式在数值上比原始多项式更稳定，龙格效应就可以消除。不，一个k次多项式形成一个向量space，所以无论你用什么基础来表示，它们都有相同的跨度！

样条曲线：分段三次多项式及其在回归中的应用

多项式回归确实有用，但我们往往需要分段多项式。最受欢迎的选择是三次样条。就像多项式有不同的表示法一样，样条曲线有很多表示法：

截断幂基础
自然三次样条基
B 样条基础

B 样条基在数值上最稳定，因为它具有紧凑的支持。因此，协方差矩阵 X'X 是带状的，因此求解正规方程 (X'X) b = (X'y) 非常稳定。

在 R 中，我们可以使用 splines 包（R 基础包之一）中的 bs 函数来获取 B 样条基。对于 bs(x)，自由度 df 的唯一数值限制是我们不能有比 length(unique(x)).

更多的基础

我不确定你的数据是什么样的，但也许你可以试试

library(splines)
model <- lm(noisy.y ~ bs(q, df = 10))

惩罚平滑/回归样条

如果您不断增加自由度，回归样条仍有可能过度拟合您的数据。最好的模型似乎就是选择最好的自由度。

一个很好的方法是使用惩罚平滑样条或惩罚回归样条，以便将模型估计和自由度选择（即"smoothness"）集成在一起。

stats 包中的 smooth.spline 函数可以做到这两点。不像它的名字所暗示的那样，在大多数情况下，它只是拟合一个惩罚回归样条而不是平滑样条。阅读 ?smooth.spline 了解更多信息。对于您的数据，您可以尝试

fit <- smooth.spline(q, noisy.y)

（注意，smooth.spline没有公式界面。）

加性惩罚样条和广义加性模型 (GAM)

一旦我们有多个协变量，我们就需要加法模型来克服 "curse of dimensionality" 同时保持合理。根据平滑函数的表示，GAM 可以有多种形式。在我看来，最受欢迎的是 R 推荐的 mgcv 包。

您仍然可以使用 mgcv 拟合单变量惩罚回归样条：

library(mgcv)
fit <- gam(noisy.y ~ s(q, bs = "cr", k = 10))

R 中的高（或非常高）阶多项式回归（或替代方案？）

High (or very high) order polynomial regression in R (or alternatives?)

regression

r

linear-regression

lm

polynomials