在
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")
xn <- c(.341,.122,.476,.981)
x0 <- 1
dx <- 1/30
for (j in 0:2) {
cols <- 1+j*9 +1:9
i <- floor(xn[j+1]*30)
w1 <- (xn[j+1]-i*dx)/dx
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
se <- sqrt(x0%*%b$Vp%*%t(x0));se
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
参数控制它们放置的位置。
predict
。我使用package::rms
因为它可以让您执行所有这些操作。 - IRTFMfit <- 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)
- IRTFMlrm(formula = mortality.under.2 ~ rcs(birth_year, 8) + rcs(maternal_age, 3) + +wealth2 + sex + residence + maternal_educ + birth_order, data = colombia2)
可以工作,但specs(gam.2)
只给出每个区间中多项式的节点位置。 - TomFunction()
的结果来查看最佳拟合情况。这可能比仅运行模型要复杂一些。我不明白为什么您认为specs()
可以与rms模型一起使用。也许我一开始不应该提供一个离题的替代方案。 - IRTFM