【问题标题】:How to extract fitted splines from a GAM (`mgcv::gam`)如何从 GAM 中提取拟合样条曲线(`mgcv::gam`)
【发布时间】:2013-03-23 07:47:07
【问题描述】:

我正在使用 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 的大量文档,但我仍不确定。 任何建议都非常感谢。

【问题讨论】:

  • “提取样条线”并不容易。虽然我很高兴被证明是错误的。出于目的 2)您可以在网格上使用 predict。”我使用 package::rms 是因为它可以让您执行所有这些操作。
  • 谢谢,但是您将如何使用 rms 做到这一点?
  • 缩短了相当多的准备工作并对变量结构进行了一些猜测:fit &lt;- 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)
  • 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)只给我节点位置,每个区间的多项式。
  • 您可以指定结的位置,也可以使用Function() 结果查看最适合的位置。它可能比仅仅运行模型要复杂一些。我不明白为什么你会认为specs() 可以与 rms 模型一起使用。也许我一开始就不应该提供一个切线的替代方案。

标签: r gam mgcv


【解决方案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 向量/矩阵选择适当的值。这些对应于模型中nth 项的值。因此,将您想要固定的值设置为某个平均值,然后改变与样条关联的值。

如果您在 R 中执行所有这些操作,则仅在您拥有数据的样条协变量值处评估样条曲线会更容易进入另一个模型。您可以通过创建一个数据框来进行预测,然后使用

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

其中mod 是拟合的 GAM 模型(通过 mgcv::gam),newdat 是包含模型中每个变量的列的数据框(包括参数项;设置您不想要的项变化到某个恒定的平均值[比如数据集中变量的平均值]或某个水平(如果是一个因素)。 type = "terms" 部分将为newdat 中的每一行返回一个矩阵,其中包含模型中每个项的拟合值的“贡献”,包括样条项。只需取该矩阵中与样条线对应的列 - 再次居中。

也许我误解了您的 Q1。如果要控制结,请参阅mgcv::gamknots 参数。默认情况下,mgcv::gam 在数据的极端处放置一个结,然后剩余的“结”均匀分布在该区间上。 mgcv::gam 没有找到结 - 它会为您放置它们,您可以通过 knots 参数控制它们放置的位置。

【讨论】:

  • 这是一个非常有帮助的答案。由于我不能轻易捐赠额外积分,因此我将看看是否可以找到您的一些早期答案以进行投票。应该不会太难。 Gavin,您是一位知识渊博的优秀教师。
  • 这是一个非常好的解释。我的问题确实不清楚。我想做一个混合程序。我想放置一个或两个结而不是在特定位置让程序放置任何需要的剩余结;有可能的?谢谢
  • @AntonioPedroRamos 就像我说的,mgcv::gam 所做的唯一一件事就是将结放置在端点之间并均匀地放置在它们之间。如果要选择一些结位置,则需要自己定位所有结。 IIRC 这些惩罚回归模型对节点位置不是很敏感。
猜你喜欢
  • 2016-01-07
  • 1970-01-01
  • 1970-01-01
  • 2020-03-21
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多