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

Question

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

rgammgcv

18

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

因此，我有两个问题：

如何拟合平滑曲线以便强制一个结点位于特定位置，同时让模型找到其他结点？
如何从拟合的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的广泛文档，但仍不确定。非常感谢您的任何建议。

- Tom

“提取样条曲线”并不那么容易，虽然我很乐意被证明是错误的。对于目的2），您可以在网格上使用predict。我使用package::rms因为它可以让您执行所有这些操作。 - IRTFM

谢谢，但您如何使用rms来实现这一点？ - Tom

短路法可以省去很多准备工作，并对变量结构进行一些猜测：

fit <- lrm(mortality.under.2 ~ rcs(maternal_age_c, 3) + rcs(birth_year, 3) %ia% rcs(wealth2, 3) + sex + residence + maternal_educ + birth_order, data=colombia2)); Function(fit)

- IRTFM

lrm(formula = mortality.under.2 ~ rcs(birth_year, 8) + rcs(maternal_age, 3) + +wealth2 + sex + residence + maternal_educ + birth_order, data = colombia2)

可以工作，但 specs(gam.2) 只给出每个区间中多项式的节点位置。 - Tom

您可以指定节点位置，或者使用Function()的结果来查看最佳拟合情况。这可能比仅运行模型要复杂一些。我不明白为什么您认为specs()可以与rms模型一起使用。也许我一开始不应该提供一个离题的替代方案。 - IRTFM

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Gavin Simpson · Accepted Answer

在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是通过mgcv::gam拟合的GAM模型，newdat是包含模型中每个变量的列的数据帧（包括参数项；将您不想变化的项设置为某个固定平均值[例如数据集中该变量的平均值]或某个级别，如果是因子）。type = "terms"部分将返回一个矩阵，该矩阵对于newdat中的每一行，都会显示模型中每个项对拟合值的“贡献”，包括样条项。只需取此矩阵对应于样条的列-再次进行中心化。

也许我误解了您的Q1。如果您想要控制节点，请参见mgcv::gam的knots参数。默认情况下，mgcv::gam在数据的极端处放置一个节点，然后将其余的“节点”均匀分布在区间上。 mgcv::gam不会“找到”节点-它会为您放置它们，并且您可以通过knots参数控制它们放置的位置。