在 MGCV 中移除惩罚三次样条中二次项的惩罚？

Question

我想用 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]

估计中的惩罚？

我曾尝试在样条函数

m

中使用参数

s()

，因为

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) # this is always true

Answer 1

这是一个将三次项添加到二次模型的方法：

> x <- seq(0,1, length = 100)
> y <- 0.5*x + x^2 + rnorm(100)
> mod1 <- gam(
+     y~s(x, fx = F, k = 3, bs = "cr")+ I(x^3)
+ )
> mod1$coefficients
(Intercept)      I(x^3)      s(x).1      s(x).2 
-1.70442708  9.32342373  0.03978659 -5.49355285 
> summary(mod1)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(x, fx = F, k = 3, bs = "cr") + I(x^3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)  
(Intercept)   -1.704      1.014  -1.680   0.0961 .
I(x^3)         9.323      3.999   2.331   0.0218 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
     edf Ref.df     F p-value
s(x) 1.8   1.96 1.925   0.125

R-sq.(adj) =  0.255   Deviance explained = 27.6%
GCV = 0.94177  Scale est. = 0.90598   n = 100

这是情节的输出：

png( ); plot(mod1) ; dev.off()

比较：

> mod2 <- gam(
+     y~s(x, fx = F, k = 3, bs = "cr")
+ )
> summary(mod2)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(x, fx = F, k = 3, bs = "cr")

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.64997    0.09781   6.646 1.74e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
       edf Ref.df     F  p-value    
s(x) 1.433  1.679 14.62 1.35e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.213   Deviance explained = 22.4%
GCV = 0.98046  Scale est. = 0.9566    n = 100
> png( ); plot(mod2) ; dev.off()
quartz 
     2