如何用metafor复现mixmeta多变量Meta分析？解析结果差异

问题概述 #

尝试用metafor包的rma.mv()函数复现mixmeta包的多变量Meta分析模型，单变量分析（mixmeta::mixmeta()与metafor::rma.uni()）结果可完全匹配，但多变量模型在固定效应估计值、模型logLik等指标上存在细微差异。以mixmeta::berkey98数据集为例，即使均指定REML估计方法、复合对称（CS）相关结构，结果仍有小数点后第二位的差异，需明确差异原因并实现结果对齐。

差异原因分析 #

• REML计算细节差异：两个包对REML估计的自由度调整、对数似然函数的常数项处理存在细微差别，mixmeta直接基于多变量联合似然优化，而metafor的rma.mv()在处理长格式数据时，优化路径的参数化方式略有不同。
• 优化算法与收敛阈值：默认情况下，mixmeta和metafor使用的优化算法（如BFGS、Nelder-Mead）或收敛判定阈值不同，可能导致收敛到相近但不完全一致的局部最优解。
• 协方差矩阵参数化差异：虽然均指定CS结构，但mixmeta的bscov='cs'是直接对多变量随机效应协方差矩阵施加CS约束，而metafor的struct='CS'针对分组内的效应关联，参数化初始值设定可能不同。

实现结果对齐的方法 #

• 统一优化参数：在rma.mv()中指定与mixmeta一致的优化算法，调整收敛阈值，缩小优化结果差异：

  # 调整rma.mv的优化控制参数，匹配mixmeta的默认设置
  rma.mv(yi=as.numeric(yi_long$Estimate), V=bldiag(v_list),  
         mods = ~ 0+Measures, struct = 'CS',  
         data=yi_long, random=~Measures | STUDY,
         control = list(optimizer="BFGS", rel.tol=1e-10, iter.max=1000))
  

• 明确指定REML方法：虽然两个包默认都是REML，但显式指定method="REML"避免默认行为差异：
  • mixmeta：mixmeta(..., method="reml")
  • metafor：rma.mv(..., method="REML")
• 检查数据格式一致性：确保长格式数据的分组（STUDY）、效应量（yi）和协方差矩阵（V）完全对应，避免数据整理过程中出现的细微偏差（如矩阵维度、顺序错误）。

可复现代码 #

mixmeta模型代码 #

library(mixmeta)
library(metafor)
library(data.table)

# mixmeta多变量模型
summary(mixmeta(cbind(PD,AL), S=berkey98[5:7], data=berkey98, bscov='cs', method="reml"))

metafor模型代码（含优化参数调整） #

# 数据整理为长格式
yi <- data.table(cbind(PD=berkey98$PD, AL=berkey98$AL, STUDY=rownames(berkey98)))
yi_long <- melt(yi, id.vars = 'STUDY', measure.vars = c('PD', 'AL'), 
                value.name = 'Estimate', variable.name = 'Measures')

# 构建协方差矩阵列表
v_list <- list()
for (i in 1:nrow(berkey98)){
  V <- matrix(c(berkey98$var_PD[i], berkey98$cov_PD_AL[i], 
                berkey98$cov_PD_AL[i], berkey98$var_AL[i]), ncol=2)
  v_list[[i]] <- V
}

# 调整优化参数后的rma.mv模型
summary(rma.mv(yi=as.numeric(yi_long$Estimate), V=bldiag(v_list),  
               mods = ~ 0+Measures, struct = 'CS',  
               data=yi_long, random=~Measures | STUDY,
               method="REML",
               control = list(optimizer="BFGS", rel.tol=1e-10, iter.max=1000)))

参考文献 #

加斯帕里尼A、阿姆斯特朗B、肯沃德MG（2012）。非线性及其他多参数关联的多变量Meta分析。《医学统计学》，31(29):3821–3839。

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