它不是三次回归样条，但它确实清楚地分离了问题的两个组成部分。在这里，我将展示如何使用薄板样条（mgcv::gam()中的默认基础）来实现这一点。
使用您的设置和我的软件包来处理GAMs

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

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

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

注意函数9、10和11在基函数中是如何消失的。这将允许我们隔离样条中的二次分量之外的摆动，将2阶多项式摆动留给模型中的参数项。

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)

以下是向二次模型添加三次项的方法：

> 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

下面是plot的输出：

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