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

[Hanna Mütze]

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\theta_{ij}can be expressed as:

θij=μ+ζij(2)+ζj(3)+ϵij\theta_{ij} = \mu + \zeta^{(2)}_{ij} + 
\zeta^{(3)}_j + \epsilon_{ij}

Where:

  * μ\m: the pooled true effect size,
  * ζij(2)\zeta^{(2)}_{ij}​: study-specific random deviation,
  * ζj(3)\zeta^{(3)}_j​: effect size-specific random deviation, and
  * ϵij\epsilon_{ij}​: residual error.

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=1es.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]

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\theta_{11}correct?
 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?
 3. Is it generally appropriate to visualize a three-level model in a
    forest plot, or would an alternative approach be more suitable?

I appreciate your guidance and suggestions.

Kind regards,
*Hanna Mütze*

------------------------------------------------------------------------

	[[alternative HTML version deleted]]