## 13.2 Forest Plots for Bayesian Meta-Analysis

As you have seen before, Bayesian models allow us to extract their sampled posterior distribution, which can be extremely helpful to directly assess the probability of specific values given our model. We can exploit this feature of Bayesian models to create enhanced Forest Plots which are both very informative and pleasing to the eye (see blog posts by A. Salomon Kurz and Matti Vuorre for more information). Unfortunately, there is currently no maintained “pre-packaged” function to create such forest plots, but it is possible to build one oneself using functions of the tidybayes package. So, let us first load this package along with a few other ones before we proceed.

library(tidybayes)
library(dplyr)
library(ggplot2)
library(ggridges)
library(glue)
library(stringr)
library(forcats)

Before we can generate the plot, we have to prepare the data. In particular, we need to extract the posterior distribution for each study individually (since forest plots also depict the specific effect size of each study). To achieve this, we can use the spread_draws function in the tidybayes package. The function needs three arguments as input: our fitted brms model, the random-effects factor by which the results should be indexed, and the parameter we want to extract (here b_Intercept, since we want to extract the fixed term, the effect size). Using the pipe operator %>%, we can directly manipulate the output; using the mutate function in dplyr, we calculate the actual effect size of each study by adding the pooled effect size b_Intercept to the estimated deviation for each study. We save the result as study.draws.

study.draws <- spread_draws(m.brm, r_Author[Author,], b_Intercept) %>%
mutate(b_Intercept = r_Author + b_Intercept)

Next, we want to generate the distribution of the pooled effect in a similar way (since in forest plots, the summary effect is usally displayed in the last row). We therefore slightly adapt the code from before, dropping the second argument to only get the pooled effect. The call to mutate only adds an extra column called “Author”. We save the result as pooled.effect.draws.

pooled.effect.draws <- spread_draws(m.brm, b_Intercept) %>%
mutate(Author = "Pooled Effect")

Next, we bind study.draws and pooled.effect.draws together in one data frame. We then start a pipe again, calling ungroup first, and then using mutate to (1) clean the study labels (i.e., replace dots with spaces), and (2) reorder the study factor levels by effect size. This is the data we will need for plotting, so we save it as forest.data.

forest.data <- bind_rows(study.draws, pooled.effect.draws) %>%
ungroup() %>%
mutate(Author = str_replace_all(Author, "[.]", " ")) %>%
mutate(Author = reorder(Author, b_Intercept))

Lastly, we also need summarized data (the mean and credibility interval) of each study. To do this, we use our newly generated forest.data dataset, group it by Author, and then use the mean_qi function to calculate the values. We save the output as forest.data.summary.

forest.data.summary <- group_by(forest.data, Author) %>%
mean_qi(b_Intercept)

$~$

Here is the complete code for you to copy, paste, and adapt:

study.draws <- spread_draws(m.brm, r_Author[Author,], b_Intercept) %>%
mutate(b_Intercept = r_Author + b_Intercept)

pooled.effect.draws <- spread_draws(m.brm, b_Intercept) %>%
mutate(Author = "Pooled Effect")

forest.data <- bind_rows(study.draws, pooled.effect.draws) %>%
ungroup() %>%
mutate(Author = str_replace_all(Author, "[.]", " ")) %>%
mutate(Author = reorder(Author, b_Intercept))

forest.data.summary <- group_by(forest.data, Author) %>%
mean_qi(b_Intercept)

$~$

We are now ready to generate the forest plot using the ggplot2 package. The code to generate the plot looks like this:

ggplot(aes(b_Intercept, relevel(Author, "Pooled Effect", after = Inf)),
data = forest.data) +
geom_vline(xintercept = fixef(m.brm)[1, 1], color = "grey", size = 1) +
geom_vline(xintercept = fixef(m.brm)[1, 3:4], color = "grey", linetype = 2) +
geom_vline(xintercept = 0, color = "black", size = 1) +
geom_density_ridges(fill = "blue", rel_min_height = 0.01, col = NA, scale = 1,
alpha = 0.8) +
geom_pointintervalh(data = forest.data.summary, size = 1) +
geom_text(data = mutate_if(forest.data.summary, is.numeric, round, 2),
aes(label = glue("{b_Intercept} [{.lower}, {.upper}]"), x = Inf), hjust = "inward") +
labs(x = "Standardized Mean Difference",
y = element_blank()) +
theme_minimal()

Looks good!

One thing is very important to mention here: the effect sizes displayed in the forest plot do not represent the observed effect sizes (i.e., the effect sizes we find reported in the original studies), but the estimate of the “true” effect size ($$\theta_k$$) of a study based on the Bayesian model. The effect sizes shown in the forest plot are equivalent to the study-wise estimates we saw when extracting the random effects using ranef(m.brm) (except that these were centered around the pooled effect size).