Meta-análise multivariada mixmeta vs metafor

Gostaria de reproduzir uma meta-análise multivariada de mixmeta (Gasparrini et al., 2012) usando metafor (reconheço que houve uma pergunta semelhante antes, com resultados diferentes, em uma meta-análise multinível com dois pacotes semelhantes em R (metafor e mixmeta), que, no entanto, não foi respondida, pois nenhum exemplo reproduzível foi fornecido). Especificamente, espero reproduzir o modelo de meta-análise produzido por mixmeta() usando rma.mv().

Análises univariadas desses dois pacotes (mixmeta::mixmeta() e metafor::rma.uni()) podem ser facilmente comparadas. No entanto, os modelos multivariados produzidos por esses dois pacotes apresentaram diversas discrepâncias. Usando o conjunto de dados mixmeta::berkey98 como exemplo (veja o exemplo reproduzível abaixo), embora mixmeta() e rma.mv() tenham especificado REML como método de estimativa e CS (simetria composta) como estrutura de correlação (os resultados PD e AL do mesmo estudo são correlacionados), os modelos resultantes apresentam diferentes estimativas de efeito fixo, logLik dos modelos, etc. Embora as diferenças sejam muito pequenas (principalmente na segunda casa decimal), isso ajuda a entender por que essas discrepâncias ocorrem.

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

^{Criado em 2025-04-08 pelo pacote reprex (v2.0.1)}

Eu apreciaria qualquer conselho sábio sobre o porquê disso acontecer e como produzir modelos mais próximos (se não idênticos) a partir desses dois pacotes.

Ref.: Gasparrini A, Armstrong B, Kenward MG (2012). Meta-análise multivariada para associações não lineares e outras associações multiparâmetros. Estatística em Medicina. 31(29):3821–3839.