我正在使用GAM在Logistic回归中对时间趋势进行建模。但是我想从中提取拟合的样条,以将其添加到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的详尽文档,但仍不确定。
任何建议都非常感谢。

最佳答案

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向量/矩阵选择适当的值。这些对应于模型中第cc项的值。因此,将要固定的值设置为某个平均值,然后更改与样条曲线关联的值。

如果您在R中执行所有这些操作,则仅使用样条协变量的值来评估样条就容易得多,因为该样条协变量的数据已进入另一个模型。您可以通过创建一个数值框架来进行预测,然后使用

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


其中,n是拟合的GAM模型(通过mod),mgcv::gam是数据帧,其中包含该模型中每个变量的列(包括参数项;设置您不想更改的项)常数平均值(例如,数据集中变量的平均值)或一定水平(如果有因素的话)。 newdat部分将为type = "terms"中的每一行返回一个矩阵,并对模型中每个项(包括样条项)的拟合值“贡献”。只需取该矩阵对应于样条线的列-再次将其居中即可。

也许我误解了你的Q1。如果要控制结,请参见newdatknots参数。默认情况下,mgcv::gam在数据的极端处打一个结,然后剩余的“结”在时间间隔内均匀分布。 mgcv::gam找不到结-它为您放置了结,您可以通过mgcv::gam参数控制将结放置在何处。

关于r - 如何从GAM(`mgcv::gam`)中提取拟合的样条曲线,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/15584541/

10-15 03:49