Support _metafor_ package (predictions and plots)

[arthur-albuquerque]

Hi,

It seems that model_parameters() isn't calculating the correct point estimate (and 95% CI) for all subgroups in a metafor::rma() object. Below, I show that predict() and emmeans() yield the same results for all subgroups, while model_parameters() only recovers the "alternate" subgroup correctly.

# remotes::install_github("wviechtb/metafor") # install version 3.9-5
library(metafor)
library(emmeans)
library(parameters)

dat <- metadat::dat.bcg 

datRR <- metafor::escalc(measure="RR",
                         data=dat,
                         ai=tpos, bi=tneg, ci=cpos, di=cneg)

res <- metafor::rma(data=datRR,
                         method = "REML",
                         mods = ~ alloc,
                         yi = yi, vi = vi,
                         slab = author)

predict(res, newmods = rbind(c(1,0), # Random effects
                             c(0,0), # alternate
                             c(0,1)), # systematic
        digits=3) 
#> 
#>     pred    se  ci.lb  ci.ub  pi.lb pi.ub 
#> 1 -0.966 0.267 -1.490 -0.442 -2.255 0.324 
#> 2 -0.518 0.441 -1.383  0.347 -1.980 0.944 
#> 3 -0.429 0.345 -1.105  0.247 -1.788 0.930

gr <- metafor::emmprep(res) # version 3.9-5 required

emmeans::emmeans(gr, specs = ~alloc) 
#>  alloc      emmean    SE  df asymp.LCL asymp.UCL
#>  alternate  -0.518 0.441 Inf     -1.38     0.347
#>  random     -0.966 0.267 Inf     -1.49    -0.442
#>  systematic -0.429 0.345 Inf     -1.10     0.247
#> 
#> Results are given on the log (not the response) scale. 
#> Confidence level used: 0.95

parameters::model_parameters(res)
#> # Random Effects
#> 
#> Parameter         | Coefficient |   SE |         95% CI |      z |      p | Weight
#> ----------------------------------------------------------------------------------
#> Aronson           |       -0.89 | 0.57 | [-2.01,  0.23] |  -1.56 | 0.119  |   3.07
#> Ferguson & Simes  |       -1.59 | 0.44 | [-2.45, -0.72] |  -3.59 | < .001 |   5.14
#> Rosenthal et al.1 |       -1.35 | 0.64 | [-2.61, -0.08] |  -2.09 | 0.036  |   2.41
#> Hart & Sutherland |       -1.44 | 0.14 | [-1.72, -1.16] | -10.19 | < .001 |  49.97
#> Vandiviere et al  |       -1.62 | 0.47 | [-2.55, -0.70] |  -3.43 | < .001 |   4.48
#> TPT Madras        |        0.01 | 0.06 | [-0.11,  0.14] |   0.19 | 0.849  | 252.42
#> Coetzee & Berjak  |       -0.47 | 0.24 | [-0.94,  0.00] |  -1.98 | 0.048  |  17.72
#> Overall           |       -0.45 | 0.52 | [-1.46,  0.56] |  -0.87 | 0.385  |       
#> 
#> # alternate
#> 
#> Parameter            | Coefficient |   SE |         95% CI |     z |      p | Weight
#> ------------------------------------------------------------------------------------
#> Frimodt-Moller et al |       -0.22 | 0.23 | [-0.66,  0.23] | -0.96 | 0.336  |  19.53
#> Stein & Aronson      |       -0.79 | 0.08 | [-0.95, -0.62] | -9.46 | < .001 | 144.81
#> Overall              |       -0.52 | 0.44 | [-1.38,  0.35] | -1.17 | 0.240  |       
#> 
#> # systematic
#> 
#> Parameter          | Coefficient |   SE |         95% CI |     z |      p | Weight
#> ----------------------------------------------------------------------------------
#> Rosenthal et al.2  |       -1.37 | 0.27 | [-1.90, -0.84] | -5.07 | < .001 |  13.69
#> Comstock et al.1   |       -0.34 | 0.11 | [-0.56, -0.12] | -3.05 | 0.002  |  80.57
#> Comstock & Webster |        0.45 | 0.73 | [-0.98,  1.88] |  0.61 | 0.541  |   1.88
#> Comstock et al.2   |       -0.02 | 0.27 | [-0.54,  0.51] | -0.06 | 0.948  |  14.00
#> Overall            |        0.09 | 0.56 | [-1.01,  1.19] |  0.16 | 0.874  |
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.

sessionInfo()
#> R version 4.2.1 (2022-06-23)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Catalina 10.15.7
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] parameters_0.18.2 emmeans_1.8.0     metafor_3.9-5     metadat_1.2-0    
#> [5] Matrix_1.4-1     
#> 
#> loaded via a namespace (and not attached):
#>  [1] rstudioapi_0.14    knitr_1.40         mathjaxr_1.6-0     magrittr_2.0.3    
#>  [5] insight_0.18.2     xtable_1.8-4       lattice_0.20-45    rlang_1.0.4       
#>  [9] fastmap_1.1.0      stringr_1.4.1      tools_4.2.1        grid_4.2.1        
#> [13] nlme_3.1-157       xfun_0.32          coda_0.19-4        cli_3.3.0         
#> [17] bayestestR_0.12.1  withr_2.5.0        htmltools_0.5.3    datawizard_0.5.1  
#> [21] yaml_2.3.5         digest_0.6.29      lifecycle_1.0.1    fs_1.5.2          
#> [25] glue_1.6.2         evaluate_0.16      rmarkdown_2.16     reprex_2.0.2      
#> [29] stringi_1.7.8      estimability_1.4.1 compiler_4.2.1     mvtnorm_1.1-3

^{Created on 2022-09-03 with reprex v2.0.2}

[strengejacke]

model_parameters() does not report predictions (that should be in the modelbased package, but I guess not implemented yet). It reports the estimates, SE etc. from the summary. See:

library(metafor)
#> Loading required package: Matrix
#> Loading required package: metadat
#> 
#> Loading the 'metafor' package (version 3.8-1). For an
#> introduction to the package please type: help(metafor)
library(parameters)

dat <- metadat::dat.bcg 

datRR <- metafor::escalc(measure="RR",
                         data=dat,
                         ai=tpos, bi=tneg, ci=cpos, di=cneg)

res <- metafor::rma(data=datRR,
                         method = "REML",
                         mods = ~ alloc,
                         yi = yi, vi = vi,
                         slab = author)

summary(res)
#> 
#> Mixed-Effects Model (k = 13; tau^2 estimator: REML)
#> 
#>   logLik  deviance       AIC       BIC      AICc   
#> -10.3301   20.6602   28.6602   29.8705   36.6602   
#> 
#> tau^2 (estimated amount of residual heterogeneity):     0.3615 (SE = 0.2111)
#> tau (square root of estimated tau^2 value):             0.6013
#> I^2 (residual heterogeneity / unaccounted variability): 88.77%
#> H^2 (unaccounted variability / sampling variability):   8.91
#> R^2 (amount of heterogeneity accounted for):            0.00%
#> 
#> Test for Residual Heterogeneity:
#> QE(df = 10) = 132.3676, p-val < .0001
#> 
#> Test of Moderators (coefficients 2:3):
#> QM(df = 2) = 1.7675, p-val = 0.4132
#> 
#> Model Results:
#> 
#>                  estimate      se     zval    pval    ci.lb   ci.ub    
#> intrcpt           -0.5180  0.4412  -1.1740  0.2404  -1.3827  0.3468    
#> allocrandom       -0.4478  0.5158  -0.8682  0.3853  -1.4588  0.5632    
#> allocsystematic    0.0890  0.5600   0.1590  0.8737  -1.0086  1.1867    
#> 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
model_parameters(res)
#> # Random Effects
#> 
#> Parameter         | Coefficient |   SE |         95% CI |      z |      p | Weight
#> ----------------------------------------------------------------------------------
#> Aronson           |       -0.89 | 0.57 | [-2.01,  0.23] |  -1.56 | 0.119  |   3.07
#> Ferguson & Simes  |       -1.59 | 0.44 | [-2.45, -0.72] |  -3.59 | < .001 |   5.14
#> Rosenthal et al.1 |       -1.35 | 0.64 | [-2.61, -0.08] |  -2.09 | 0.036  |   2.41
#> Hart & Sutherland |       -1.44 | 0.14 | [-1.72, -1.16] | -10.19 | < .001 |  49.97
#> Vandiviere et al  |       -1.62 | 0.47 | [-2.55, -0.70] |  -3.43 | < .001 |   4.48
#> TPT Madras        |        0.01 | 0.06 | [-0.11,  0.14] |   0.19 | 0.849  | 252.42
#> Coetzee & Berjak  |       -0.47 | 0.24 | [-0.94,  0.00] |  -1.98 | 0.048  |  17.72
#> Overall           |       -0.45 | 0.52 | [-1.46,  0.56] |  -0.87 | 0.385  |       
#> 
#> # alternate
#> 
#> Parameter            | Coefficient |   SE |         95% CI |     z |      p | Weight
#> ------------------------------------------------------------------------------------
#> Frimodt-Moller et al |       -0.22 | 0.23 | [-0.66,  0.23] | -0.96 | 0.336  |  19.53
#> Stein & Aronson      |       -0.79 | 0.08 | [-0.95, -0.62] | -9.46 | < .001 | 144.81
#> Overall              |       -0.52 | 0.44 | [-1.38,  0.35] | -1.17 | 0.240  |       
#> 
#> # systematic
#> 
#> Parameter          | Coefficient |   SE |         95% CI |     z |      p | Weight
#> ----------------------------------------------------------------------------------
#> Rosenthal et al.2  |       -1.37 | 0.27 | [-1.90, -0.84] | -5.07 | < .001 |  13.69
#> Comstock et al.1   |       -0.34 | 0.11 | [-0.56, -0.12] | -3.05 | 0.002  |  80.57
#> Comstock & Webster |        0.45 | 0.73 | [-0.98,  1.88] |  0.61 | 0.541  |   1.88
#> Comstock et al.2   |       -0.02 | 0.27 | [-0.54,  0.51] | -0.06 | 0.948  |  14.00
#> Overall            |        0.09 | 0.56 | [-1.01,  1.19] |  0.16 | 0.874  |
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald z-distribution approximation.

^{Created on 2022-09-03 with reprex v2.0.2}

The "Overall" rows correspond to the three estimates from the summary.

[arthur-albuquerque]

I see, thank you for the quick reply!

I don't think that users expect to get results from summary(), especially not when plotting results with see.

Don't you agree this forest plot looks odd by depicting coefficients, instead of predictions by subgroup?

[strengejacke]

Well, it's similar to forestplot():

library(metafor)
#> Loading required package: Matrix
#> Loading required package: metadat
#> 
#> Loading the 'metafor' package (version 3.8-1). For an
#> introduction to the package please type: help(metafor)
library(parameters)

dat <- metadat::dat.bcg 

datRR <- metafor::escalc(measure="RR",
                         data=dat,
                         ai=tpos, bi=tneg, ci=cpos, di=cneg)

res <- metafor::rma(data=datRR,
                         method = "REML",
                         mods = ~ alloc,
                         yi = yi, vi = vi,
                         slab = author)

model_parameters(res) |> plot(show_labels = TRUE)

forest(res)

^{Created on 2022-09-03 with reprex v2.0.2}

[arthur-albuquerque]

metafor is an excellent package, but there isn't a built-in easy way to plot results from subgroups.

By contrast, the R package meta has a nice forest() function that quickly plots by subgroup:

I think parameters/see could fill this gap in metafor's functionalities by showing model predictions instead of summary results.

What do you think?

[strengejacke]

I'd say we should have both. But this would indeed be something for modelbased.

[strengejacke]

@bwiernik @DominiqueMakowski just to let you know. I don't know how to implement a predict-method for rma, there's also a years-old issue for ggeffects: strengejacke/ggeffects#30

[arthur-albuquerque]

Perhaps the new function metafor::emmprep can help. Its features were described by @wviechtb in this thread.

[bwiernik]

I have work in progress on improving support for metafor in easystats: easystats/insight#514

[vincentarelbundock]

Should we have insight::get_predicted.rma? (I know nothing about these models, so I probably can't help. Sorry.)

[bwiernik]

The major headache currently is that metafor doesn't currently accept data frames as input in its predict method. wviechtb/metafor#48