I'd like to reproduce multivariate meta-analysis from mixmeta (Gasparrini et al, 2012) using metafor (I acknowledge that there had been one similar question before different results in multilevel meta-analysis with two similar packages in R (metafor and mixmeta) which was, however, not answered as no reproducible example was provided). Specifically, I am hoping to reproduce meta-analysis model produced by mixmeta() using rma.mv().
Univariate analyses from these two packages(mixmeta::mixmeta() and metafor::rma.uni()) can be easily matched. However, multivariate models produced by these two packages had several discrepancies. Using mixmeta::berkey98 dataset as an example (see below reproducible example), although both mixmeta() and rma.mv() specified REML as the estimation method and CS (compound symmetry) as the correlation structure (the outcomes PD, AL from the same study are correlated), the resulting models have different fix effect estimates, logLik of the models etc.... While the differences are very small (in the second decimal place mostly), it helps to understand why these discrepancies occur.
library(mixmeta)
#> This is mixmeta 1.2.0. For an overview type: help('mixmeta-package').
library(metafor)
#> Loading required package: Matrix
#>
#> Loading the 'metafor' package (version 3.0-2). For an
#> introduction to the package please type: help(metafor)
#>
#> Attaching package: 'metafor'
#> The following object is masked from 'package:mixmeta':
#>
#> blup
library(data.table)
# Multivariate
summary(mixmeta(cbind(PD,AL), S=berkey98[5:7], data=berkey98, bscov='cs')) ### , method='reml' is default, cov str = unstr default
#> Call: mixmeta(formula = cbind(PD, AL), S = berkey98[5:7], data = berkey98,
#> bscov = "cs")
#>
#> Multivariate random-effects meta-analysis
#> Dimension: 2
#> Estimation method: REML
#>
#> Fixed-effects coefficients
#> Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
#> PD 0.3636 0.0788 4.6163 0.0000 0.2092 0.5180 ***
#> AL -0.3380 0.0782 -4.3229 0.0000 -0.4912 -0.1847 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#>
#> Random-effects (co)variance components
#> Structure: Compound symmetry
#> Std. Dev Corr
#> PD 0.1582 PD
#> AL 0.1582 0.529
#>
#> Multivariate Cochran Q-test for heterogeneity:
#> Q = 128.2267 (df = 8), p-value = 0.0000
#> I-square statistic = 93.8%
#>
#> 5 units, 2 outcomes, 10 observations, 2 fixed and 2 random-effects parameters
#> logLik AIC BIC
#> 3.3106 1.3789 1.6966
#prep data for metafor
yi=data.table(cbind(PD=berkey98$PD, AL=berkey98$AL, STUDY=rownames(berkey98), pubyear=berkey98$pubyear))
yi_long=melt(yi, id.vars = 'STUDY', measure.vars = c('PD', 'AL'),
value.name = c('Estimate'), variable.name = 'Measures')
v_list=list()
for (i in 1:nrow(berkey98)){
V <- matrix(0, nrow=2, ncol=2)
V=matrix(c(berkey98$var_PD[i], berkey98$cov_PD_AL[i], berkey98$cov_PD_AL[i], berkey98$var_AL[i]), ncol=2 ); V
v_list[[i]]=V
}
summary(rma.mv(yi=as.numeric(yi_long$Estimate), V=bldiag(v_list), mods = ~ 0+Measures, struct = 'CS',
data=yi_long, random=~Measures | STUDY )) ### ###
#>
#> Multivariate Meta-Analysis Model (k = 10; method: REML)
#>
#> logLik Deviance AIC BIC AICc
#> 3.0997 -6.1994 1.8006 2.1184 15.1340
#>
#> Variance Components:
#>
#> outer factor: STUDY (nlvls = 5)
#> inner factor: Measures (nlvls = 2)
#>
#> estim sqrt fixed
#> tau^2 0.0273 0.1651 no
#> rho 0.5447 no
#>
#> Test for Residual Heterogeneity:
#> QE(df = 8) = 142.8194, p-val < .0001
#>
#> Test of Moderators (coefficients 1:2):
#> QM(df = 2) = 64.3617, p-val < .0001
#>
#> Model Results:
#>
#> estimate se zval pval ci.lb ci.ub
#> MeasuresPD 0.3786 0.0813 4.6558 <.0001 0.2192 0.5379 ***
#> MeasuresAL -0.3361 0.0848 -3.9630 <.0001 -0.5023 -0.1699 ***
#>
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Created on 2025-04-08 by the reprex package (v2.0.1)
I would appreciate anyone's wise advice on why this is the case and producing closer (if not the identical) models from these two packages.
Ref: Gasparrini A, Armstrong B, Kenward MG (2012). Multivariate meta-analysis for non-linear and other multi-parameter associations. Statistics in Medicine. 31(29):3821–3839.
The ordering of the rows in yi_long does not match up with v_list. If you use yi_long=yi_long[c(1,6,2,7,3,8,4,9,5,10)] before fitting the model, then you will get identical results.