13.1 Bayesian Meta-Analysis in R using the brms package

Now that we have defined the Bayesian model for our meta-analysis, it is time to implement it in R. Here, we will use the brms package (Bürkner 2017, 2018) to fit our model. The brms package is a very versatile and powerful tool to fit Bayesian regression models. It can be used for huge range of applications, including multilevel (mixed-effects) models, generalized linear models, multivariate models, nonlinear, and generalized additive models, to name just a few. Most of these applications use person-level data, but brms can also be easily used for meta-analysis, where we deal with (weighted) study-level data. We chose the brms here because it offers a vast amount of functionalities, and because it is very user friendly. Furthermore, the brms package is, at least primarily, a package for multilevel models; this means that it may help us further understand the multilevel nature of Bayesian Meta-Analysis models “by doing”. It should be noted that there are also other great resources for Bayesian Meta-Analysis, for example the bayesmeta package (Röver 2017), which may be worth looking into.

13.1.1 Fitting a Bayesian Meta-Analysis Model

Before we start fitting the model, we first have to install and load the brms package. If you are having trouble installing or running functions in the package, you might find this forum helpful.


In this example, I will use my ThirdWave dataset, which contains data of a real-world meta-analysis investigating the effects of “Third-Wave” psychotherapies in college students. The data is identical to the madata dataset we used in Chapter 4. The dataset is part of the dmetar package. If you have the package installed already, you have to load it into your library first.


If you don’t want to use the dmetar package, you can download the madata dataset in the Datasets section of this guide.

Before we fit the model, let us first define the priors for the true pooled effect size \(\mu\) and the between-study heterogeneity \(\tau\). Previously, we specified that \(\mu \sim \mathcal{N}(0,1)\) and \(\tau \sim \mathcal{HC}(0,0.5)\). We can use the prior function to specify priors. The function takes two arguments: in the first argument, we specify the distribution we want to assume for our prior, including the distribution parameters. In the second argument, we have to define the class of the prior. For \(\mu\), the appropriate class is Intercept, since it is a fixed population-level effect. For \(\tau\), the class is sd, because it is a variance. We can define both priors using the prior function, then concatenate them using c(), and save the resulting object as priors. The code then looks like this:

priors <- c(prior(normal(0,1), class = Intercept),
            prior(cauchy(0,0.5), class = sd))

We can then proceed and fit the model. To do this, we can use the brm function in the brms package. The function has many parameters, but only a few are relevant for us:

  • formula: In this argument, the formula for the model is specified. brms uses a regression formula notation, in which an outcome (in our case, the effect size) y is predicted by one or more variables x. A tilde (~) is used to specify that there is a predictive relationship: y ~ x. Meta-analyses are somewhat special, because we do not have a predictor variable of the effect size (unless when performing meta-regression), so x has to be replaced with 1, indicating an intercept only model. Furthermore, we cannot simply use the effect size of each study in isolation as the predicted response y; we also have to give studies with a larger sample size (and thus greater precision of the effect size estimate) a greather weight. This can be done by using y|se(se_y) instead of only y, where the se(se_y) part stands for the standard error of each effect size y in our data. If we want to use a random-effects model, the last step is to add a random-effects term (1|study) to the predictor part of the formula. This specifies that we assume that the effect sizes in y are nested within studies, the true effects of which are themselves random draws from an overarching population of true effect sizes. If we want to use a fixed-effect model, we can simply omit this term. The full formula for a random-effects model may therefore look something like this: y|se(se_y) ~ 1 + (1|random). For more information on the formula setup for brm models, you can type ?brmsformula in your Console and then hit Enter to open the documentation.
  • data: The dataset to be used for our meta-analysis.
  • prior: The priors we want to use for our model. In our case, we can simply plug in the priors object we created previously here.
  • iter: The number of iterations the MCMC algorithm should run. The more complex your model, the higher this number should be. However, more iterations also mean that the function will take longer to finish.

Before we go on, let us first have a look at the structure of my ThirdWave dataset.

## tibble [18 × 3] (S3: tbl_df/tbl/data.frame)
##  $ Author: chr [1:18] "Call et al." "Cavanagh et al." "DanitzOrsillo" "de Vibe et al." ...
##  $ TE    : num [1:18] 0.709 0.355 1.791 0.182 0.422 ...
##  $ seTE  : num [1:18] 0.261 0.196 0.346 0.118 0.145 ...

We see that there are three relevant columns in the data.frame: TE, which contains the calculated effect size of each study, expressed as the Standardized Mean Difference (SMD), seTE, the standard error corresponding to each effect size, and Author, a unique identifier for each study/effect size.

I want to save my Bayesian Meta-Analysis model as m. Using this data, the code for my model looks like this:

m.brm <- brm(TE|se(seTE) ~ 1 + (1|Author),
             data = ThirdWave,
             prior = priors,
             iter = 4000)

Please be aware that Bayesian methods are much more computationally intensive compared to the standard meta-analytic techniques we covered before; it may therefore take a few minutes until the sampling is completed.

13.1.2 Assessing Convergence

Before we start analyzing the results, we have to make sure that the model has converged (i.e., if the MCMC algorithm found the optimal solution). If this is not the case, the parameters are not trustworthy and should not be interpreted. Non-convergence can happen very frequently in Bayesian models and can often be solved by rerunning the model with more iterations (iter). Generally, we should do two things: first, conduct posterior predictive checks, and secondly, check the \(\hat{R}\) values of the parameter estimates.

In posterior predictive checks, data are simulated through random draws from the posterior predictive distribution, which are then compared to the observed data. If a model has converged, we would expect that the density of the replications are roughly similar to the ones of the observed data. This can easily be controlled by using the pp_check function.


The \(\hat{R}\) value, on the other hand, represents the Potential Scale Reduction Factor (PSRF) we already covered when discussing Bayesian Network Meta-Analysis. Here, the value of \(\hat{R}\) should be smaller than \(1.01\). To check this, we can produce a summary of our m.brm object, which also includes the \(\hat{R}\) under Rhat.

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: TE | se(seTE) ~ 1 + (1 | Author) 
##    Data: ThirdWave (Number of observations: 18) 
## Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
##          total post-warmup samples = 8000
## Group-Level Effects: 
## ~Author (Number of levels: 18) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.30      0.10     0.12     0.53 1.00     2095     3600
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     0.57      0.09     0.40     0.77 1.00     3638     3998
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

As we can see, the Rhat value for both parameters is 1, signifying convergence. This means that we can start interpreting the results.

13.1.3 Interpreting the Results

We can start interpreting the results by looking at the Group-Level Effects in our summary output first. This section is reserved for random effects we defined in our formula (+ (1|Author)). Since we are dealing with a random-effects meta-analysis here, the variable ~Author, signifying the individual studies, has been modeled with a random intercept. As we described before, this represents our assumption on level 2 that each study has its own “true” effect size, which has been sampled from an overarching distribution of true effect sizes. We also see that our group-level effect has 18 levels, corresponding with the \(k=18\) studies in our data. We see that the Estimate of our between-study heterogeneity (sd(Intercept)) is \(\tau = 0.30\), closely resembling our initial “best guess” when setting the priors. By using the ranef function, we can also extract the estimated deviation of each study’s “true” effect size from the pooled effect:

## $Author
## , , Intercept
##                           Estimate Est.Error         Q2.5       Q97.5
## Call et al.             0.07181028 0.1993251 -0.319820598  0.47649403
## Cavanagh et al.        -0.14767075 0.1789226 -0.517921502  0.19531573
## DanitzOrsillo           0.50460835 0.2909379  0.008174737  1.12549786
## de Vibe et al.         -0.32597554 0.1466809 -0.619205195 -0.04989067
## Frazier et al.         -0.12048853 0.1536381 -0.433784299  0.17057881
## Frogeli et al.          0.03365368 0.1727736 -0.306290768  0.37758921
## Gallego et al.          0.08856291 0.1875901 -0.269807296  0.47297695
## Hazlett-Stevens & Oren -0.03170411 0.1804180 -0.398259361  0.31158212
## Hintz et al.           -0.21150417 0.1680367 -0.560700839  0.10479442
## Kang et al.             0.29672080 0.2455310 -0.130471985  0.83139982
## Kuhlmann et al.        -0.31742818 0.1906900 -0.713570495  0.02557018
## Lever Taylor et al.    -0.11168048 0.1938160 -0.510486158  0.25080754
## Phang et al.           -0.02093020 0.1970425 -0.416875655  0.36790605
## Rasanen et al.         -0.08342763 0.2009354 -0.494478087  0.29966525
## Ratanasiripong         -0.02812745 0.2292966 -0.492509640  0.42087368
## Shapiro et al.          0.41080179 0.2610767 -0.043115453  0.96235589
## Song & Lindquist        0.01986460 0.1854680 -0.343915792  0.38516697
## Warnecke et al.         0.01422909 0.1940359 -0.367426711  0.39775144

The next part of the output we want to interpret are the Population-Level Effects. This section represents the “fixed” population effects we modeled. In our case, the population we are dealing with is the overarching distribution of effect sizes in our meta-analysis, and the fixed population effect is its mean \(\mu\): the pooled effect size.

In the output, we see that the Estimate of the pooled effect size is \(SMD = 0.57\), with the 95% credibility interval (not confidence interval!) ranging from \(95\% CrI: 0.40-0.77\). This indicates that there is in fact a moderate-sized overall effect of the interventions studied in this meta-analysis. Because this is a Bayesian Model after all, you won’t find any \(p\)-values here, but this example should underline that we can also make reasonable inferences without having to resort to significance testing. A great thing we can do in Bayesian Meta-Analysis, but not in frequentist meta-analysis, is model the parameters we want to estimate probabilistically. The Bayesian model not only estimates the parameters of interest, but a whole posterior distribution for \(\tau\) and \(\mu\), which we can access quite easily. To do this, we can use the posterior_samples function.

post.samples <- posterior_samples(m.brm, c("^b", "^sd"))
## [1] "b_Intercept"          "sd_Author__Intercept"

The resulting data.frame contains two columns: b_Intercept, the posterior sample data for the pooled effect size, and sd_Author_Intercept, the one for the between-study heterogeneity \(\tau\). We rename the columns smd and tau to make the name more informative.

names(post.samples) <- c("smd", "tau")

Using the ggplot2 package, we can make a density plot of the posterior distributions, using the data in post.samples. The code to do this looks like this:


# Plot for SMD
ggplot(aes(x = smd), data = post.samples) +
  geom_density(fill = "lightblue", color = "lightblue", alpha = 0.7) +
  geom_point(y = 0, x = mean(post.samples$smd)) +
  labs(x = expression(italic(SMD)),
       y = element_blank()) +

# Plot for tau
ggplot(aes(x = tau), data = post.samples) +
  geom_density(fill = "lightgreen", color = "lightgreen", alpha = 0.7) +
  geom_point(y = 0, x = mean(post.samples$tau)) +
    labs(x = expression(tau),
       y = element_blank()) +

This is the output we get:

We see that the posterior distributions follow a unimodal, and roughly normal distribution, peaking around the values for \(\mu\) and \(\tau\) we saw in the output. The fact that Bayesian methods create an actual sampling distribution for our parameters of interest means that we can calculate exact probabilities that \(\mu\) or \(\tau\) will be larger than some specific value. For example, we may have found in previous literature that if effects of an intervention are below \(SMD=0.30\), they are not meaningful anymore for the type of patients we are studying. We could therefore calculate the probability that the pooled effect is in fact \(SMD=0.30\) or smaller, based on our model. This can be done by looking at the Empirical Cumulative Distribution Function (ECDF), which lets us select one specific value \(X\), and returns the probability of some value \(x\) being smaller than \(X\) based on the data. The ECDF of the posterior distribution for the pooled effect size can be seen below:

We can use the ecdf function to implement the ECDF in R, and then check the probability of our pooled effect being smaller than 0.30. The code looks like this.

smd.ecdf <- ecdf(post.samples$smd)
## [1] 0.002125

We see that with 0.21%, the probability of our pooled effect being smaller than \(SMD = 0.30\) is very, very low. In this scenario, this would mean that the effects of the interventions we find in this meta-analysis are very likely to be meaningful.


Bürkner, Paul-Christian. 2017. “brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.

Bürkner, Paul-Christian. 2018. “Advanced Bayesian Multilevel Modeling with the R Package brms.” The R Journal 10 (1): 395–411. https://doi.org/10.32614/RJ-2018-017.

Röver, Christian. 2017. “Bayesian Random-Effects Meta-Analysis Using the Bayesmeta R Package.” arXiv Preprint 1711.08683, November. http://www.arxiv.org/abs/1711.08683.