如何获取metafor::rma.mv()对应lme4::VarCorr()结果以计算伪R²

核心思路 #

rma.mv()生成的模型对象中，m1$vvc存储了随机效应的方差-协方差分量（多水平随机斜率模型下，会包含各水平的随机截距、斜率方差及二者协方差），我们可以基于这个矩阵提取所需信息，还原出类似lme4/nlme包VarCorr()的输出格式。

实操代码 #

假设你的模型是包含随机截距+随机斜率的多水平模型，示例模型定义如下：

library(metafor)
# 示例：两水平（study/site）含treatment随机斜率的元回归模型
m1 <- rma.mv(yi, vi, 
             mods = ~ treatment,
             random = ~ treatment | study/site,
             data = dat)

运行以下代码提取并格式化随机效应方差信息：

# 提取随机效应方差-协方差矩阵
vc_matrix <- m1$vvc
# 获取各效应的名称（含水平和效应类型）
effect_names <- rownames(vc_matrix)

# 拆分水平与效应类型，整理基础方差/标准差数据
varcorr_result <- data.frame(
  分组 = gsub("\\:.*", "", effect_names),
  效应 = gsub(".*\\:", "", effect_names),
  方差 = diag(vc_matrix),
  标准差 = sqrt(diag(vc_matrix))
)

# 提取并添加协方差信息（如果存在）
if(nrow(vc_matrix) > 1) {
  cov_values <- vc_matrix[lower.tri(vc_matrix)]
  if(length(cov_values) > 0) {
    cov_pairs <- apply(which(lower.tri(vc_matrix), arr.ind = TRUE), 1, 
                       function(x) paste(effect_names[x[1]], effect_names[x[2]], sep = "-"))
    cov_df <- data.frame(
      分组 = rep(varcorr_result$分组[1], length(cov_values)),
      效应 = cov_pairs,
      协方差 = cov_values
    )
    varcorr_result <- merge(varcorr_result, cov_df, by = c("分组", "效应"), all.x = TRUE)
  }
}

# 输出类似VarCorr的格式
cat("随机效应方差-协方差矩阵\n")
for(group in unique(varcorr_result$分组)) {
  cat(sprintf("\n分组: %s\n", group))
  group_data <- varcorr_result[varcorr_result$分组 == group, ]
  print(group_data[, c("效应", "方差", "标准差", "协方差")], na.print = "")
}

关键说明 #

• m1$vvc是对称矩阵：对角线元素对应各随机效应的方差，非对角线元素是不同随机效应（如截距与斜率）的协方差
• 代码通过正则拆分效应名称，区分不同水平（如study、site）和效应类型（如(Intercept)、treatment）
• 最终输出按水平分组展示，完全适配VarCorr()的信息结构，可直接用于文献中的伪R²计算逻辑

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