带有 mgcv 包的 GAM 随机效应的方差-协方差矩阵

Variance-covariance matrix for random effect from GAM with mgcv package

提取lme4包的随机效应和随机效应的方差-协方差矩阵如下:

library(lme4)
fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fm1.rr <- ranef(fm1,condVar=TRUE)
fm1.pv <- attr(rr[[1]],"postVar")

我想知道如何使用 mgcv 执行此操作? 'gam.vcomp'函数确实提取了估计的方差分量,但不是针对每个水平的随机效应。

library(mgcv)
fm2 <- gam(Reaction ~ Days + s (Subject, bs="re"), data = sleepstudy, method = "REML")
gam.vcomp(fm2)
library(lme4)
data(sleepstudy)

fm1 <- lmer(Reaction ~ Days + (1|Subject), sleepstudy)
fm1.rr <- ranef(fm1,condVar=TRUE)$Subject[,1]
fm1.pv <- sqrt(attr(ranef(fm1,condVar=TRUE) [['Subject']],"postVar")[1,1,])

library(mgcv)
fm2 <- gam(Reaction ~ Days + s (Subject, bs="re"), 
data = sleepstudy,   method = "REML")

为每个 Subject

提取随机效果
idx <-grep("Subject", names(coef(fm2)))
fm2.rr<-coef(fm2)[idx]
attributes(fm2.rr)<-NULL

我们可以看到两个模型中的随机效应与预期相同。

为了提取随机效应的方差-协方差矩阵并计算误差,我们使用参数 Vp,它是贝叶斯后验协方差矩阵:

fm2.pv <-sqrt(diag(fm2$Vp))[idx]

或者频率估计的协方差矩阵Ve

fm2.pv <-sqrt(diag(fm2$Ve))[idx]

我们可以看到用 mgcv 估计的随机效应误差与用 lme4 模型估计的随机效应误差略有不同。基于贝叶斯后验协方差矩阵的误差较大,而基于频率矩阵的误差较小。

您还可以使用包 gamm4,它基于 gamm 包,但在下面使用 lme4。该模型将拟合为:

fm3 <- gamm4(Reaction ~ Days, random = ~ (1|Subject), data = sleepstudy)

可以按照正常的lme4程序获得随机效应和随机效应的方差-协方差矩阵。

fm3.rr <- ranef(fm3$mer,condVar=TRUE)
fm3.pv <- attr(fm3.rr[[1]],"postVar")[1,1,]

但是 gamm4 可能比 gam 慢得多,因此请阅读帮助文件以了解何时最适合您的需要。