如何从 GAM 中提取拟合样条线（`mgcv::gam`）

Question

我正在使用 GAM 对逻辑回归中的时间趋势进行建模。然而我想从中提取拟合的样条线以将其添加到另一个模型，该模型无法在 GAM 或 GAMM 中拟合。

因此我有两个问题：

1. 如何随着时间的推移拟合平滑器，以便强制一个结位于特定位置，同时让模型找到其他结？

2. 如何从拟合的 GAM 中提取矩阵，以便可以将其用作不同模型的插补？

我运行的模型类型如下：

gam <- gam(mortality.under.2~ maternal_age_c+ I(maternal_age_c^2)+
           s(birth_year,by=wealth2) + wealth2 + sex +
           residence + maternal_educ + birth_order,
           data=colombia2, family="binomial")

我已经阅读了 GAM 的大量文档，但我仍然不确定。 任何建议都非常感激。

Answer 1

在

mgcv::gam

中有一种方法可以做到这一点（你的Q2），通过

predict.gam

方法和

type = "lpmatrix"

。

?predict.gam

甚至有一个例子，我在下面复制：

 library(mgcv)
 n <- 200
 sig <- 2
 dat <- gamSim(1,n=n,scale=sig)
 
 b <- gam(y ~ s(x0) + s(I(x1^2)) + s(x2) + offset(x3), data = dat)
 
 newd <- data.frame(x0=(0:30)/30, x1=(0:30)/30, x2=(0:30)/30, x3=(0:30)/30)
 
 Xp <- predict(b, newd, type="lpmatrix")

 ##################################################################
 ## The following shows how to use use an "lpmatrix" as a lookup 
 ## table for approximate prediction. The idea is to create 
 ## approximate prediction matrix rows by appropriate linear 
 ## interpolation of an existing prediction matrix. The additivity 
 ## of a GAM makes this possible. 
 ## There is no reason to ever do this in R, but the following 
 ## code provides a useful template for predicting from a fitted 
 ## gam *outside* R: all that is needed is the coefficient vector 
 ## and the prediction matrix. Use larger `Xp'/ smaller `dx' and/or 
 ## higher order interpolation for higher accuracy.  
 ###################################################################
 
 xn <- c(.341,.122,.476,.981) ## want prediction at these values
 x0 <- 1         ## intercept column
 dx <- 1/30      ## covariate spacing in `newd'
 for (j in 0:2) { ## loop through smooth terms
   cols <- 1+j*9 +1:9      ## relevant cols of Xp
   i <- floor(xn[j+1]*30)  ## find relevant rows of Xp
   w1 <- (xn[j+1]-i*dx)/dx ## interpolation weights
   ## find approx. predict matrix row portion, by interpolation
   x0 <- c(x0,Xp[i+2,cols]*w1 + Xp[i+1,cols]*(1-w1))
 }
 dim(x0)<-c(1,28) 
 fv <- x0%*%coef(b) + xn[4];fv    ## evaluate and add offset
 se <- sqrt(x0%*%b$Vp%*%t(x0));se ## get standard error
 ## compare to normal prediction
 predict(b,newdata=data.frame(x0=xn[1],x1=xn[2],
         x2=xn[3],x3=xn[4]),se=TRUE)

这贯穿了整个过程，甚至是在 R 或 GAM 模型之外完成的预测步骤。您必须稍微修改示例才能执行您想要的操作，因为该示例评估模型中的所有项，并且除了样条线之外还有其他两项 - 本质上您做同样的事情，但仅针对样条项，这涉及查找样条线的

Xp

矩阵的相关列和行。然后您还应该注意样条线是居中的，因此您可能也可能不想撤消它。

对于 Q1，为示例中的

xn

向量/矩阵选择适当的值。这些对应于模型中第

n

项的值。因此，将您想要固定的值设置为某个平均值，然后改变与样条线关联的值。

如果您在 R 中执行所有这些操作，则只需在您拥有要进入另一个模型的数据的样条协变量的值来评估样条就会更容易。您可以通过创建一个要预测的值的数据框来实现这一点，然后使用 predict(mod, newdata = newdat, type = "terms")

其中 

mod

是拟合的 GAM 模型（通过

mgcv::gam

），

newdat

是包含模型中每个变量的列的数据框（包括参数项；设置您不想更改的项）某个恒定平均值（即数据集中变量的平均值）或某个水平（如果是一个因素）。

type = "terms"

部分将为

newdat

中的每一行返回一个矩阵，其中包含对模型中每个项（包括样条项）的拟合值的“贡献”。只需取该矩阵中与样条线相对应的列即可 - 再次居中。也许我误解了你的问题1。如果您想控制结，请参阅

knots

的

mgcv::gam

参数。默认情况下，

mgcv::gam

在数据的极值处放置一个结，然后将剩余的“结”分布在数据的均匀间隔的分位数处。

mgcv::gam

不会

找到

结 - 它会为您放置它们，您可以通过 knots 参数控制它放置它们的位置。