如何从GAM中提取拟合样条(``mgcv：：GAM‘)

Question

我正在用GAM来模拟时间趋势的logistic回归。然而，我想从它中提取拟合样条，将它添加到另一种模型中，这在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中，有一种方法可以通过predict.gam方法和type = "lpmatrix"来实现(您的Q2)。

?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中完成的，那么只在样条协变量的值处计算样条就更容易了，因为你有数据，这是进入另一个模型的。要做到这一点，您可以创建一个值的数据框架，用于预测at，然后使用

predict(mod, newdata = newdat, type = "terms")

如果mod是适合的GAM模型(通过mgcv::gam)，则newdat是包含模型中每个变量的列的数据框架(包括参数项；设置不想更改到某些常量平均值的项，例如数据集中变量的平均值，或者如果某个因素的话，则为某一级别)。type = "terms"部分将为newdat中的每一行返回一个矩阵，并为模型中的每个项(包括样条项)的拟合值“贡献”。只需取这个矩阵的列，它对应于样条--同样，它是中心的。

也许我误解了你的Q1。如果要控制结，请参见knots参数到mgcv::gam。默认情况下，mgcv::gam在数据的极值处打结，然后剩余的“纽结”均匀地分布在间隔上。mgcv::gam没有找到结-它为您放置它们，您可以通过knots参数控制它们的位置。