[R-meta] Forestplot for Bayesian Three-Level-Model

[Viechtbauer, Wolfgang (NP)]

Dear Hanna,

Please see below for my responses.

Best,
Wolfgang

> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of Hanna Mütze via R-sig-meta-analysis
> Sent: Friday, January 24, 2025 16:19
> To: r-sig-meta-analysis using r-project.org
> Cc: Hanna Mütze <muetzeh using uni-bremen.de>
> Subject: [R-meta] Forestplot for Bayesian Three-Level-Model
>
> Dear all,
>
> I am attempting to create a forest plot for a Bayesian three-level
> meta-analysis. The plot should visualize the estimated true effect sizes
> along with their corresponding 95% credible intervals (CrIs).
> Optionally, I would like to include density distributions for each
> observed effect size. The mathematical formulation is described here:
> https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/multilevel-ma.html
> <https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/multilevel-
> ma.html#multilevel-ma>.
>
> Specifically, each observed effect size θ_ij can be expressed as:
>
> Θ_ij = μ + ζ_ij(2) + ζ_j(3) + ϵ_ij
>
> Where:
>
>   * μ: the pooled true effect size,
>   * ζ_ij(2): study-specific random deviation,
>   * ζ_j(3): effect size-specific random deviation, and
>   * ϵ_ij​: residual error.

ζ_ij(2) is the effect size specific random effect and ζ_j(3) the study specific random effect.

> Using the *|Chernobyl|* dataset from the |dmetar| package and a Bayesian
> model fitted via |brms|, I am trying to map the statistical components
> in this equation to the output of the |brms| model. Below is the code I
> have used to define the model:
>
> # Load libraries
> library(pacman)
> p_load(dmetar, brms, rstan, posterior)
> # Set Stan options
> options(mc.cores = parallel::detectCores())
> rstan_options(auto_write = TRUE)
> # Load data data("Chernobyl", package = "dmetar")
> head(Chernobyl)
> # Define priors
> prior <- c( prior(normal(0, 1), class = Intercept), prior(cauchy(0, 0.5), class = sd))
> # Fit the three-level model
> l3.model <- brm( z | se(se.z) ~ 1 + (1 | author/es.id), data = Chernobyl, prior
> = prior, seed = 20121220, control = list(max_treedepth = 15, adapt_delta =
> 0.999), iter = 4000, warmup = 1000, sample_prior = "yes", save_pars =
> save_pars(all = TRUE))
> # Summarize posterior draws
> l3.model_sum_draws <- summarise_draws(l3.model)
>
> To compute the modeled effect size for, for instance, es.id = 1,
> I believe the formula should be:
>
> θ11= b_Intercept + r_author[Aghajanyan.&Suskuv.(2009),Intercept] +
> r_author:es.id[Aghajanyan.&Suskuv.(2009)_id_1,Intercept]

Correct.

> I would calculate this based on the mean column in the
> |l3.model_sum_draws| data frame. My interpretation is that this would
> extend the two-level forest plot described here:
> https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/bayesian-ma.html
>
> However, I am uncertain about the following:
>
>  1. Is this interpretation of the hierarchical structure and the
>     calculation of θ11 correct?

Yes. In essence:

with(l3.model_sum_draws,
l3.model_sum_draws[variable == "b_Intercept","mean"] +
l3.model_sum_draws[variable == "r_author[Aghajanyan.&.Suskov.(2009),Intercept]", "mean"] +
l3.model_sum_draws[variable == "r_author:es.id[Aghajanyan.&.Suskov.(2009)_id_1,Intercept]", "mean"])

Note that this takes the mean for each posterior distribution. But if you want to whole posterior distribution, then you need to add the individual draws:

l3.draws <- as_draws_df(l3.model)

which you can then again summarize:

mean(l3.draws$b_Intercept + l3.draws$"r_author[Aghajanyan.&.Suskov.(2009),Intercept]" + l3.draws$"r_author:es.id[Aghajanyan.&.Suskov.(2009)_id_1,Intercept]")
sd(l3.draws$b_Intercept + l3.draws$"r_author[Aghajanyan.&.Suskov.(2009),Intercept]" + l3.draws$"r_author:es.id[Aghajanyan.&.Suskov.(2009)_id_1,Intercept]")

or compute quantiles of:

quantile(l3.draws$b_Intercept + l3.draws$"r_author[Aghajanyan.&.Suskov.(2009),Intercept]" + l3.draws$"r_author:es.id[Aghajanyan.&.Suskov.(2009)_id_1,Intercept]", prob=c(.025, .975))

>  2. What is the role or meaning of the |z_1[1,1]| output in the
>     |l3.model_sum_draws|? Should it be included in the forest plot
>     calculation?

This has to do with some technical aspects of how brm() translates the model into Stan code. In essence, you can ignore this.

>  3. Is it generally appropriate to visualize a three-level model in a
>     forest plot, or would an alternative approach be more suitable?

Sure, why not? You may want to indicate which estimates 'belong' to the same study. See the last examples at:

https://wviechtb.github.io/metafor/reference/forest.rma.html

Of course, if the number of estimates is very large, then the forest plot also gets very large (in terms of its height).

> I appreciate your guidance and suggestions.
>
> Kind regards,
> *Hanna Mütze*