在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]估计的惩罚？

我尝试通过操纵样条函数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) # this is always true

- Lacey Etzkorn

---

在lm中正确的方法可能是使用poly(x,2)，然后与poly(x,3)进行比较，但这可能不适用于mgcv方法。您可以使用I(x^2)和I(x^3)将二次和三次项强制输入公式，但这样会失去poly的正交特性。在您所在领域的“标准方法”是否了解当二次和三次项输入回归公式时出现的统计问题？（注意：在R中，您不能使用x^2和x^3来实现此目的。） - IRTFM

Answer 1

它不是一个立方回归样条，但它确实清晰地分离了您问题的两个组成部分。在这里，我展示如何使用薄板样条来实现这一点，它是mgcv::gam()中的默认基础。

使用您的设置加上我的用于处理广义加性模型（GAM）的软件包

library("mgcv")
library("gratia")
library("ggplot2")

set.seed(1)
df <- data.frame(x = seq(0,1, length = 100),
                 y = 0.5*x + x^2 + rnorm(100))

我们首先看一下带有第三阶导数惩罚的薄板样条的基础：

basis(s(x, m = 3), data = df) |>
  draw() +
  facet_wrap(~ bf)

这导致

最后三个基函数（9、10和11）位于惩罚的零空间中；它们不受惩罚的影响，因为它们在任何地方都具有0的三阶导数。函数11是一个二次项。函数9与模型截距混淆，并将通过对基础设置总和为零的约束来从基础中删除；这是默认约束条件，当拟合GAM时，gam()默认执行此操作。

假设薄板样条可以使用，则您要拟合的GAM如下：

m <- gam(y ~ s(x, m = 3), data = df, method = "REML")
summary(m)

根据数据模拟的方式，该模型使用了2个EDF，这是我们所期望的。

> summary(m)                                                                  

Family: gaussian 
Link function: identity 

Formula:
y ~ s(x, m = 3)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.6939     0.0907   7.651 1.47e-11 ***
---
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)   2      2 11.1 4.56e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =   0.17   Deviance explained = 18.6%
-REML = 134.45  Scale est. = 0.82259   n = 100

如果您想对样条比二次多项式进行正式测试，那么我们可以重新拟合上述模型，但是通过poly(x, 2)将所需的项作为参数项包含在内，但修改薄板样条的基础以消除惩罚空间中的所有函数。我们通过设置m = c(3,0)来移除零空间：

basis(s(x, m = c(3, 0)), data = df) |>
  draw() +
  facet_wrap(~ bf)

注意函数9、10和11不再在基础中找到。这将使我们能够隔离样条中超过二次的波动性组件，将二阶多项式波动性留给模型中的参数项。

m0 <- gam(y ~ poly(x, 2) + s(x, m = c(3, 0)), data = df, method = "REML")
summary(m0)

最后一行产生了

> summary(m0)                                                                 

Family: gaussian 
Link function: identity 

Formula:
y ~ poly(x, 2) + s(x, m = c(3, 0))

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.6939     0.0907   7.651 1.47e-11 ***
poly(x, 2)1   4.2437     0.9113   4.657 1.02e-05 ***
poly(x, 2)2   0.5129     0.9074   0.565    0.573    
---
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) 0.0001194      8 0    0.72

R-sq.(adj) =   0.17   Deviance explained = 18.6%
-REML = 130.47  Scale est. = 0.82259   n = 100

正如我们之前所看到的，超二次波动性非常小，我们未能拒绝零假设，即 f(x) === 0。

---

我认为您可以使用B样条基函数来实现此操作，这样可以控制多项式的阶数和惩罚项的阶数（不像立方回归样条基函数，它固定为由分段三次多项式组成，并带有二阶导数惩罚）：

# make it mildly cubic
set.seed(1)
df2 <- data.frame(x = seq(0,1, length = 100),
                  y = 0.5*x + x^2 + (0.1 * (x^3)) + rnorm(100))

m1 <- gam(y ~ s(x, bs = "bs", m = c(3,3)),
          data = df2, method = "REML")

这将惩罚具有非零三阶导数的函数，因此不应该惩罚基础的线性和二次函数。只是它如何做到这一点就不那么清楚了，而且你肯定不能像对CRS那样解释系数（对于bs =“cr”基础也不能这样解释-绘制s（x，bs =“cr”）的基础函数以了解原因）。这是使用B样条基础的m1拟合基础函数的样子。

basis(m1) |>
  draw()

拟合的函数必须看起来像，因为估计的样条是：

draw(m1)

- Gavin Simpson

Answer 2

这是一种向二次模型添加立方项的方法：

> 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 

- IRTFM