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)具有用于没有硬编码限制。但是，在实践中有两个数值限制。

• 基础函数是在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") 
  

• 相关及x^(k+1)之间x^k被越来越接近为1，k增加。这种接近速度当然取决于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的多项式形成一个向量空间，所以无论您用于表示的基础是什么，它们都具有相同的跨度！

---

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

多项式回归的确是有益的，但我们经常要分段多项式。最流行的选择是立方样条。像多项式，有不同的表示，有很多表示的为花键：

• 截短功率基础
• 自然三次样条基
• B样条基

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

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

我不知道你的数据是什么样子，但也许你可以试试

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

---

处罚平滑/回归样条

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

一个很好的方法是使用惩罚平滑样条或惩罚回归样条，以便对模型估计和自由度选择（即“平滑度”）进行整合。

smooth.spline函数在stats包中可以兼得。与其名称似乎暗示的不同，大部分时间只是拟合惩罚回归样条而不是平滑样条。请阅读?smooth.spline了解更多信息。为了您的资料，您可以尝试

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

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

---

加惩罚的键槽和广义加法模型（GAM）

一旦我们有不止一个协变量，我们需要可加模型来克服“维度的诅咒”，同时又是合理的。根据平滑功能的表示，GAM可以有各种形式。在我看来，最受欢迎的是由R推荐的mgcv包。

，您仍然可以安装一个单变量回归惩罚样条mgcv：

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