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.