如何在glmmTMB中获取残差方差-协方差矩阵？

针对glmmTMB模型，我们无法直接复用lme4的getME(model, "Lambdat")获取随机效应的相对协方差因子，但可以通过glmmTMB提供的theta参数和随机效应结构信息实现相同功能。以下是适配后的代码及说明：

基础版函数（适配无复杂残差结构的模型） #

这个函数对应lme4版本的逻辑，适用于仅含随机效应、残差为独立同分布的模型：

rescov_glmmTMB <- function(model, data) {
  # 获取随机效应设计矩阵并转置
  Z <- getME(model, "Z")
  Zt <- t(Z)
  
  # 提取残差方差
  vr <- sigma(model)^2
  
  # 构建随机效应协方差矩阵D（对应lme4中的crossprod(Lambdat)）
  reTrms <- getME(model, "reTrms")
  theta <- getME(model, "theta")
  # 利用glmmTMB内部函数生成Cholesky因子L，再计算D = L %*% t(L)
  L <- glmmTMB:::mkReTrms(reTrms$flist, reTrms$cnms, theta)$L
  D <- tcrossprod(L)
  
  # 计算随机效应贡献的方差部分
  var.b <- vr * (Zt %*% D %*% Z)
  
  # 残差的对角方差矩阵（独立同分布假设）
  sI <- vr * Matrix::Diagonal(nrow(data))
  
  # 总残差方差-协方差矩阵
  var.y <- var.b + sI
  invisible(var.y)
}

测试代码 #

用你提供的Pastes数据集验证：

library(lme4)
library(glmmTMB)
library(Matrix)

data("Pastes")
# lme4模型对比
lme.fit <- lmer(strength ~ (1|batch/cask), data=Pastes)
image(rescov(lme.fit, Pastes))

# glmmTMB模型测试
glmm.fit <- glmmTMB(strength ~ (1|batch/cask), data=Pastes)
image(rescov_glmmTMB(glmm.fit, Pastes))

进阶版函数（支持复杂残差相关结构） #

glmmTMB的优势之一是支持自定义残差相关结构（如AR1、复合对称等），此时残差不再是对角矩阵，需要调整函数以适配：

rescov_glmmTMB_full <- function(model, data) {
  Z <- getME(model, "Z")
  Zt <- t(Z)
  vr <- sigma(model)^2
  
  # 构建随机效应协方差矩阵D
  reTrms <- getME(model, "reTrms")
  theta <- getME(model, "theta")
  L <- glmmTMB:::mkReTrms(reTrms$flist, reTrms$cnms, theta)$L
  D <- tcrossprod(L)
  var.b <- vr * (Zt %*% D %*% Z)
  
  # 处理自定义残差协方差结构
  if (!is.null(getME(model, "residualCov"))) {
    # 提取模型预设的残差协方差矩阵，乘以残差方差
    sCov <- getME(model, "residualCov") * vr
  } else {
    # 默认独立同分布残差
    sCov <- vr * Matrix::Diagonal(nrow(data))
  }
  
  var.y <- var.b + sCov
  invisible(var.y)
}

测试复杂残差结构 #

比如构建带AR1残差结构的模型：

glmm.fit_ar1 <- glmmTMB(strength ~ (1|batch/cask), data=Pastes, 
                        dispformula = ~0,
                        correlation = corAR1(form = ~1|batch))
image(rescov_glmmTMB_full(glmm.fit_ar1, Pastes))

关键说明 #

• glmmTMB用theta参数来参数化随机效应的Cholesky因子，通过mkReTrms函数可以将theta转换为对应的Cholesky矩阵L，进而得到随机效应协方差矩阵D。
• 对于自定义残差结构，glmmTMB会生成residualCov矩阵，直接提取后乘以残差方差即可得到残差的协方差部分。

内容的提问来源于stack exchange，提问作者Reid