## 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.

```
install.packages("brms")
library(brms)
```

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.

`library(dmetar)`

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.

`str(ThirdWave[,1:3])`

```
## 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.

`pp_check(m.brm)`

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`

.

`summary(m.brm)`

```
## 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:

`ranef(m.brm)`

```
## $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"))
names(post.samples)
```

`## [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:

```
library(ggplot2)
# 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()) +
theme_minimal()
# 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()) +
theme_minimal()
```

**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)
smd.ecdf(0.3)
```

`## [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.

### References

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.