如何在mgcv的惩罚三次样条中移除二次项的惩罚？

问题 #

我想用R包mgcv拟合惩罚三次样条模型，要求不对截距、线性项和二次项施加任何惩罚，仅对三次项及样条基中的其他项施加惩罚。这么做是因为我所在领域的标准做法是用二次项调整x，类似lm(y~x+x^2)的形式。但我认为我的数据可能和这个模型有适度偏离，所以想拟合一个更灵活但不过度波动的模型，因此选择惩罚样条。

我目前的理解是：mgcv会自动不对截距和线性项施加惩罚，但会惩罚二次项。

我的工作模型代码如下：

x <- seq(0,1, length = 100)
y <- 0.5*x + x^2 + rnorm(100)
mod1 <- gam(
    y~s(x, fx = F, k = 5, bs = "cr")
)

调用mod1$coefficients会得到长度为5的向量，我认为这对应截距、线性项、二次项、三次项和一个三次样条项，其中mod1$coefficients[1:2]不受惩罚，mod1$coefficients[3:5]会被惩罚。请问我的理解是否正确？如果正确，该怎么修改代码以移除mod1$coefficients[3]的惩罚？

我曾尝试调整s()中的m参数（mgcv文档说这个参数会改变施加惩罚的样条函数导数），但拟合的样条完全没变化：

mod1 <- gam(
    y~s(x, fx = F, k = 10, bs = "cr")
)
mod2 <- gam(
    y~s(x, fx = F, k = 10, bs = "cr", m = c(3,3))
)
all(mod1$fitted.values == mod2$fitted.values) # 结果始终为TRUE

回答 #

1. 你的理解存在错误 #

mgcv中默认的自然三次样条（bs="cr"），惩罚的是样条的二阶导数，对应的是约束三次及更高阶的波动。而样条基中的前3个分量（对应常数、线性、二次项）本来就不在惩罚范围内——也就是说，mod1$coefficients[1:3]（截距、线性、二次项）都是无惩罚的，只有从第4个系数开始的高阶样条项才会被惩罚。你误以为二次项会被惩罚，这是对mgcv惩罚机制的误解。

2. 实现需求的正确代码 #

如果你想严格对应「固定截距、线性、二次项，仅惩罚三次及以上项」的需求，最清晰的做法是把线性和二次项单独作为固定效应纳入模型，再用惩罚样条拟合剩余的平滑部分：

x <- seq(0,1, length = 100)
y <- 0.5*x + x^2 + rnorm(100)
# 拆分固定项与惩罚平滑项
mod_correct <- gam(
    y ~ x + I(x^2) + s(x, bs = "cr", k = 5, fx = FALSE)
)

这里的s(x)会自动生成中心化的样条基，不包含常数、线性、二次项，因此它的系数仅对应三次及以上的波动部分，且这些系数会被惩罚。而单独的x和I(x^2)作为固定效应，完全不受惩罚，完美匹配你「先遵循领域标准二次模型，再用惩罚样条捕捉额外偏离」的需求。

3. 关于m参数无效果的解释 #

• m参数控制惩罚的导数阶数：默认m=c(2,0)是惩罚二阶导数（约束三次样条的平滑性），m=c(3,3)是惩罚三阶导数（约束四次样条的平滑性）。
• 你的模拟数据是严格二次函数加噪声，模型不需要额外的高阶波动就能完美拟合，因此调整m参数不会改变拟合结果。如果换成带有三次或更高阶趋势的数据，调整m参数就能看到明显的拟合差异。

内容的提问来源于stack exchange，提问作者Lacey Etzkorn