search

easystats / modelbased

:chart_with_upwards_trend: Estimate effects, contrasts and means based on statistical models
https://easystats.github.io/modelbased/
GNU General Public License v3.0
232 stars 19 forks source link

Support _metafor_ package (predictions and plots) #212

Open arthur-albuquerque opened 1 year ago

arthur-albuquerque commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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 commented 1 year ago

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:

image

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 commented 1 year ago

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

strengejacke commented 1 year ago

@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: https://github.com/strengejacke/ggeffects/issues/30

arthur-albuquerque commented 1 year ago

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

bwiernik commented 1 year ago

I have work in progress on improving support for metafor in easystats: https://github.com/easystats/insight/pull/514

vincentarelbundock commented 1 year ago

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

bwiernik commented 1 year ago

The major headache currently is that metafor doesn't currently accept data frames as input in its predict method. https://github.com/wviechtb/metafor/issues/48

wviechtb commented 1 year ago

Will try to tackle this soon.

bwiernik commented 1 year ago

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

That's the PR I link above

vincentarelbundock commented 1 year ago

Awesome

strengejacke commented 1 year ago

How can we work on that PR? The branch isn't published, and I cannot pull?

bwiernik commented 1 year ago

Nothing special about it? What happens when you try to checkout origin/metafor?

strengejacke commented 1 year ago

It says I need to fetch, but I can't because the branch isn't published. When I want to publish, it says I can't because I need to fetch latest files. Doesn't work from within vscode nor GitHub desktop app.

bwiernik commented 1 year ago

I'm not sure what that means because the branch is definitely published. I'll take a look when I am at a computer