Throughout the previous chapters we have estimated meta-analytic models with metafor using restricted maximum likelihood (REML). The model objects we obtained — overall effects, variance components, confidence intervals — are frequentist estimates.
A Bayesian meta-analysis instead returns a posterior distribution for every parameter in the model, by combining the data (through the likelihood) with prior beliefs about plausible parameter values. This re-framing brings several practical benefits when synthesising data in ecology, evolution and physiology:
Bayesian methods are well-established in ecological and evolutionary meta-analysis, not just a fashionable alternative, but are not used as often as frequentist methods. Pappalardo et al. (2020) compared traditional REML and Bayesian fits across a range of ecological meta-analytic scenarios and found that the two approaches generally recover the same overall effect, but that the Bayesian fit gives better-calibrated uncertainty when the number of studies is small or when variance components are near zero. Similarly, Nakagawa et al. (2022), in their re-analysis of the simulations of Song et al. (2020), included a Bayesian multilevel model alongside the frequentist alternatives and showed that it can help control the inflated Type I error rates that arise when effect sizes are non-independent (see the simulation figure in their Comment). Together, these papers position Bayesian multilevel models as a robust, often preferable, complement to REML for the kinds of complex, multi-strata data routinely synthesised in our field.
In this tutorial we use the brms R package (Bürkner, 2017). brms provides a lme4-style formula interface to the Stan probabilistic programming language, which means the syntax will look familiar if you have used lmer or metafor::rma.mv. We re-fit the phylogenetic multilevel meta-analysis from the previous chapter so you can directly compare the Bayesian and REML estimates on the same data.
Recall the phylogenetic multilevel meta-analytic model:
\[ \begin{aligned} y_{i} &= \mu + s_{j[i]} + spp_{k[i]} + a_{k[i]} + e_{i} + m_{i} \\ m_{i} &\sim N(0, v_{i}\textbf{I}) \\ s_{j} &\sim N(0, \tau^2\textbf{I}) \\ spp_{k} &\sim N(0, \sigma_{k}^2\textbf{I}) \\ a_{k} &\sim N(0, \sigma_{a}^2\textbf{A}) \\ e_{i} &\sim N(0, \sigma_{e}^2\textbf{I}) \end{aligned} \]
The model is structurally identical to the one we fit with metafor::rma.mv. To make it Bayesian we add priors on every estimated parameter:
\[ \begin{aligned} \mu &\sim N(0, 1) \\ \tau, \sigma_{k}, \sigma_{a}, \sigma_{e} &\sim t^{+}_{3}(0, 2.5) \end{aligned} \]
These are weakly informative priors. For an effect size on the Fisher’s \(Z_r\) scale, \(\mu \sim N(0, 1)\) says we expect the overall correlation to be near zero with a standard deviation that easily covers the range of biologically plausible effects. For each random-effect standard deviation we use a half-Student-t prior with 3 degrees of freedom and scale 2.5 — written \(t^{+}_{3}(0, 2.5)\). This is the conventional default for variance components in brms and Stan: it concentrates posterior mass on small, biologically plausible SDs while keeping heavier tails than a half-Normal, so the prior does not pull large variance components down toward zero when the data really do support them. We will discuss how to override these defaults below.
We re-use the dispersal data set (Z-transformed correlation coefficients, \(Z_r\)) from Wu and Seebacher (2022) introduced in the Effect Size tutorial and analysed with REML in the Phylogenetic Meta-analysis tutorial.
zr_data <- read.csv("https://raw.githubusercontent.com/daniel1noble/meta-workshop/gh-pages/data/ind_disp_raw_data.csv") %>%
dplyr::select(-dispersal_def) %>%
tibble::rowid_to_column("es_ID") %>%
dplyr::mutate(species_OTL = stringr::str_replace_all(species_OTL, "_", " "))
zr_data <- metafor::escalc(measure = "ZCOR", ri = corr_coeff, ni = sample_size,
data = zr_data, var.names = c("Zr", "V_Zr"))This step is identical to the Phylogenetic Meta-analysis chapter; we reproduce it here so the chapter is self-contained.
species <- sort(unique(as.character(zr_data$species_OTL)))
taxa <- rotl::tnrs_match_names(names = species)
tree <- rotl::tol_induced_subtree(ott_ids = ott_id(taxa), label_format = "name")
set.seed(1)
tree <- ape::compute.brlen(tree, method = "Grafen", power = 1)
tree <- ape::multi2di(tree, random = TRUE)
phylo_cov <- vcv(tree, cor = FALSE)
zr_data <- zr_data %>%
dplyr::mutate(species_OTL = stringr::str_replace_all(species_OTL, " ", "_"))brms will assign reasonable defaults for any parameter you do not name explicitly. To see what brms would use, ask for the prior skeleton implied by the model formula:
brms::get_prior(
Zr | se(sqrt(V_Zr)) ~ 1 +
(1 | study_ID) + (1 | es_ID) + (1 | species) + (1 | gr(species_OTL, cov = A)),
data = zr_data,
data2 = list(A = phylo_cov)
)We override the defaults with weakly-informative priors that match the formal definition above:
priors <- c(
brms::prior(normal(0, 1), class = "Intercept"),
brms::prior(student_t(3, 0, 2.5), class = "sd")
)class = "sd" applies the same prior to every group-level standard deviation. To put a different prior on, say, the phylogenetic SD only, use class = "sd", group = "species_OTL".
brmsbayes_phylo <- brms::brm(
Zr | se(sqrt(V_Zr)) ~ 1 +
(1 | study_ID) + (1 | es_ID) + (1 | species) + (1 | gr(species_OTL, cov = A)),
data = zr_data,
data2 = list(A = phylo_cov),
prior = priors,
chains = 4, iter = 3000, warmup = 1000, thin = 1, cores = 4, seed = 42,
control = list(adapt_delta = 0.95),
file = "models/bayes_phylo"
)The first render will compile and sample the model (a few minutes). Subsequent renders read the cached fit from models/bayes_phylo.rds and complete in seconds. Delete the .rds to force a refit (e.g. after editing the formula or priors).
A few things to note:
First, the | se(sqrt(V_Zr)) syntax tells brms to use the sampling variances to control for uncertainty in the effect sizes. The se stands for standard error, which is why we need to take the square root of the sampling variance. Now, there are a couple of ways we can use the se() function. Currently, we are not setting the sigma argument within the function to TRUE, i.e., se(sqrt(V_Zr), sigma = FALSE), which means that the model will not estimate a separate residual SD for the model. If we set sigma = TRUE we will allow the model to estimate a separate residual SD for the model. Instead, we are fitting an observation-level random effect (to keep consistent with metafor) which is a common way to model residual error in meta-analysis. If we had set sigma = TRUE as well the model would have estimated a separate residual SD for the model, but we would need to drop the observation-level random effect as these would be redundant. We can try that formulation below to show it’s roughly the same:
bayes_phylo2 <- brms::brm(
Zr | se(sqrt(V_Zr), sigma = TRUE) ~ 1 +
(1 | study_ID) + (1 | species) + (1 | gr(species_OTL, cov = A)),
data = zr_data,
data2 = list(A = phylo_cov),
prior = priors,
chains = 4, iter = 3000, warmup = 1000, thin = 1, cores = 4, seed = 42,
control = list(adapt_delta = 0.95),
file = "models/bayes_phylo2"
)We’ll explore below the similarities between these models.
Second, the gr() function around the phylogenetic random effect is a special syntax that tells brms to use the supplied covariance matrix A for that grouping factor. The data2 argument is where we pass the phylogenetic covariance matrix. Also note that the phylogenetic A matrix differs from metafor because it’s a covariance matrix whereas metafor expects a correlation matrix. In our specific case the distinction is purely terminological: because we built the tree with ape::compute.brlen(..., method = "Grafen", power = 1) it is ultrametric with root-to-tip distance equal to 1, so every diagonal of vcv(tree, cor = FALSE) equals 1 and the covariance matrix is numerically identical to the correlation matrix. The distinction only matters if you use a non-ultrametric tree or skip the Grafen normalisation.
Finally, the control argument is where we can tweak the MCMC sampling parameters. adapt_delta controls the target acceptance rate for the Hamiltonian Monte Carlo sampler; increasing it can help with convergence when the model is complex.
Before interpreting any posterior summary, it is important to confirm that the chains have converged. We can usually get a good sense of this from the printed summary of the model fit already but it’s also good to visually inspect the trace plots and posterior densities.
summary(bayes_phylo) Family: gaussian
Links: mu = identity
Formula: Zr | se(sqrt(V_Zr)) ~ 1 + (1 | study_ID) + (1 | es_ID) + (1 | species) + (1 | gr(species_OTL, cov = A))
Data: zr_data (Number of observations: 272)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~es_ID (Number of levels: 272)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.34 0.02 0.30 0.39 1.00 979 1716
~species (Number of levels: 79)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.12 0.06 0.01 0.24 1.01 318 649
~species_OTL (Number of levels: 74)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.07 0.06 0.00 0.21 1.00 1053 1636
~study_ID (Number of levels: 88)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.11 0.06 0.01 0.23 1.01 301 892
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.11 0.06 -0.01 0.25 1.00 1506 1621
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.00 0.00 0.00 0.00 NA NA NA
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Two things to look for in the printout:
iter, warmup, and thin more can help, but sometimes if the random effects are truly very low then ESS will remain low.A visual check of the trace plots and posterior densities is also really important to check to ensure proper mixing of the MCMC chains. They usually tell the same story as the printed diagnostics but can sometimes reveal issues that are not obvious from the numbers alone (e.g. multimodality, poor mixing, or a small number of effective samples).
plot(bayes_phylo, ask = FALSE)
Note that there appear to be some issues with mixing for some parameters (low ESS and chains are dragging), so we should play around with MCMC parameters to try to improve this, but this isn’t bad for demonstration purposes.
Recall that we also fit the same model with sigma = TRUE. We can look at how the output for that model compares to the one we have just fit:
summary(bayes_phylo2) Family: gaussian
Links: mu = identity
Formula: Zr | se(sqrt(V_Zr), sigma = TRUE) ~ 1 + (1 | study_ID) + (1 | species) + (1 | gr(species_OTL, cov = A))
Data: zr_data (Number of observations: 272)
Draws: 4 chains, each with iter = 3000; warmup = 1000; thin = 1;
total post-warmup draws = 8000
Multilevel Hyperparameters:
~species (Number of levels: 79)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.12 0.06 0.01 0.24 1.00 912 1698
~species_OTL (Number of levels: 74)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.07 0.06 0.00 0.21 1.00 2944 4164
~study_ID (Number of levels: 88)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.11 0.06 0.01 0.23 1.00 1163 3037
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.11 0.06 -0.01 0.24 1.00 4457 3475
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.34 0.02 0.29 0.39 1.00 5922 5079
Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Here, we can see that sigma is estimated and it matches the sd_es_ID__Intercept from the first model, which is the observation-level random effect. The two models are essentially the same, but the first one is more consistent with the way we have been modelling residual error in metafor. However, we can see that this is a more elegant way to model residual error for meta-analysis with brms (which we recommend) as you can see improvements in the MCMC diagnostics for the second model compared to the first one because of the different way the model is parameterised. So, we’ve solved our mixing issues!
The overall effect (intercept on the \(Z_r\) scale) and 95% credible interval:
brms::posterior_summary(bayes_phylo2, variable = "b_Intercept") Estimate Est.Error Q2.5 Q97.5
b_Intercept 0.107 0.0616 -0.0102 0.242Variance components — between-study, within-study residual (here estimated as sigma because we used se(..., sigma = TRUE) and dropped the OLRE), non-phylogenetic species, and phylogenetic species:
brms::posterior_summary(bayes_phylo2, variable = c("sd_study_ID__Intercept",
"sigma",
"sd_species__Intercept",
"sd_species_OTL__Intercept")) Estimate Est.Error Q2.5 Q97.5
sd_study_ID__Intercept 0.1065 0.0631 0.00593 0.231
sigma 0.3398 0.0241 0.29490 0.388
sd_species__Intercept 0.1237 0.0644 0.00824 0.245
sd_species_OTL__Intercept 0.0669 0.0564 0.00242 0.215Is there an overall positive or negative relationship between physiology and movement?
The posterior mean for the intercept is 0.107 with a 95% credible interval of -0.01 to 0.242. As in the REML analysis, the overall relationship is positive and the credible interval nearly excludes zero. Compare these numbers with the REML estimates from the previous chapter — they should be fairly similar. However, the Bayesian fit also gives us full posterior distributions for all parameters. For example, we can look at every variance component with
brms::mcmc_plot(bayes_phylo2, variable = "^sd_", regex = TRUE).
One of the key benefits of the Bayesian framework is that we can propagate uncertainty into any derived quantity simply by computing it draw-by-draw across the posterior. As an example, we will compute a full posterior distribution for \(I^2_{total}\) — the proportion of total variance attributable to genuine heterogeneity (rather than sampling noise) — using the same definition as in the Heterogeneity chapter:
\[ I^2_{total} = \frac{\tau^2 + \sigma_{k}^2 + \sigma_{a}^2 + \sigma_{e}^2}{\tau^2 + \sigma_{k}^2 + \sigma_{a}^2 + \sigma_{e}^2 + \sigma_{m}^2} \]
The denominator includes the typical sampling variance \(\sigma^2_m\), which is a fixed function of the inverse-variance weights \(w_i = 1/v_i\) and is computed once from the data (not estimated by the model):
\[ \sigma_{m}^2 = \frac{(k-1)\sum w_i}{(\sum w_i)^2 - \sum w_i^2} \]
Omitting \(\sigma^2_m\) would inflate \(I^2_{total}\) because the denominator would no longer represent total variance. With \(\sigma^2_m\) in hand, we compute the posterior of \(I^2_{total}\) by squaring each posterior draw of the random-effect standard deviations, summing them, and dividing by the total. We first want to extract the posterior draws from the model object. We know we should have 8000 draws (rows) in total (4 chains x 2000 iterations) but the effective number of draws may be lower due to autocorrelation. We can extract the posterior draws as a data frame with brms::as_draws_df(), which gives us one column per parameter and one row per draw.
# Typical sampling variance
w <- 1 / zr_data$V_Zr
k <- nrow(zr_data)
sigma2_m <- (k - 1) * sum(w) / (sum(w)^2 - sum(w^2))
# Posterior draws of every random-effect SD; square to get variances!
# Remember, for bayes_phylo2 the within-study residual variance is estimated as
# `sigma` (we used se(..., sigma = TRUE) and dropped the (1 | es_ID)
# observation-level random effect).
draws <- brms::as_draws_df(bayes_phylo2)
sigma2_study <- draws$sd_study_ID__Intercept^2
sigma2_es <- draws$sigma^2
sigma2_spp <- draws$sd_species__Intercept^2
sigma2_phylo <- draws$sd_species_OTL__Intercept^2
var_total <- sigma2_study + sigma2_es + sigma2_spp + sigma2_phylo
I2_total <- var_total / (var_total + sigma2_m)
# Posterior summary (median + 95% credible interval)
quantile(I2_total, c(0.025, 0.5, 0.975)) 2.5% 50% 97.5%
0.863 0.890 0.917 Because we computed \(I^2_{total}\) per draw, we have the posterior distribution of \(I^2_{total}\) rather than a single point estimate. The 95% interval reported above is a credible interval that fully reflects uncertainty in every variance component simultaneously — something that is awkward to obtain from the REML fit without bootstrapping.
We can also now nicely visualise the posterior distribution of \(I^2_{total}\):
ggplot(data.frame(I2_total), aes(x = I2_total)) +
geom_density(fill = "skyblue", alpha = 0.5) +
geom_vline(xintercept = quantile(I2_total, c(0.025, 0.5, 0.975)), linetype = "dashed") +
labs(x = expression(I^2[total]), y = "Density") +
theme_classic()
Adapt the code above to compute the posterior of \(I^2_{study}\) — the fraction of total variance attributable to between-study heterogeneity alone. Report the posterior median and 95% credible interval, and contrast it with \(I^2_{total}\).
Only the numerator changes; we re-use
sigma2_study,var_totalandsigma2_mfrom the chunk above. Formally:\[ I^2_{study} = \frac{\tau^2}{\tau^2 + \sigma_{k}^2 + \sigma_{a}^2 + \sigma_{e}^2 + \sigma_{m}^2} \]
Code
I2_study <- sigma2_study / (var_total + sigma2_m) quantile(I2_study, c(0.025, 0.5, 0.975))2.5% 50% 97.5% 0.000209 0.063309 0.281499
The posterior median for \(I^2_{study}\) is 0.063 with a 95% credible interval of 0 to 0.281. Comparing this with \(I^2_{total}\) (median 0.89) tells you how much of the real heterogeneity is concentrated at the study level versus spread across the other variance components (within-study residual, non-phylogenetic species, phylogenetic). You can compute \(I^2\) for any single component the same way — just put its squared SD in the numerator.
This is one of the powerful features of Bayesian meta-analysis; the ability to derive full posterior distributions for any quantity of interest by computing it draw-by-draw across the posterior. You can use this same approach to get a posterior distribution for the overall effect size, prediction intervals, or any other derived quantity you like. When you have moderators, it also means you can get overall means for each group or even contrasts for specific comparisons between levels of a moderator that maybe the model doesn’t present by default.
R version 4.5.3 (2026-03-11)
Platform: aarch64-apple-darwin20
attached base packages: stats, graphics, grDevices, utils, datasets, methods and base
other attached packages: bayesplot(v.1.15.0), mathjaxr(v.2.0-0), pander(v.0.6.6), flextable(v.0.9.9), orchaRd(v.2.2.0), lubridate(v.1.9.5), forcats(v.1.0.1), stringr(v.1.6.0), dplyr(v.1.2.1), purrr(v.1.2.2), readr(v.2.2.0), tidyr(v.1.3.2), tibble(v.3.3.1), ggplot2(v.4.0.3), tidyverse(v.2.0.0), phytools(v.2.5-2), maps(v.3.4.3), rotl(v.3.1.1), ape(v.5.8-1), brms(v.2.23.0), Rcpp(v.1.1.1-1.1), metafor(v.5.0-1), numDeriv(v.2016.8-1.1), metadat(v.1.4-0) and Matrix(v.1.7-4)
loaded via a namespace (and not attached): RColorBrewer(v.1.1-3), tensorA(v.0.36.2.1), jsonlite(v.2.0.0), magrittr(v.2.0.5), TH.data(v.1.1-3), estimability(v.1.5.1), farver(v.2.1.2), rmarkdown(v.2.31), ragg(v.1.5.2), vctrs(v.0.7.3), askpass(v.1.2.1), htmltools(v.0.5.9), progress(v.1.2.3), distributional(v.0.6.0), curl(v.7.1.0), DEoptim(v.2.2-8), StanHeaders(v.2.32.10), htmlwidgets(v.1.6.4), plyr(v.1.8.9), sandwich(v.3.1-1), emmeans(v.2.0.3), zoo(v.1.8-14), uuid(v.1.2-2), igraph(v.2.2.1), lifecycle(v.1.0.5), iterators(v.1.0.14), pkgconfig(v.2.0.3), R6(v.2.6.1), fastmap(v.1.2.0), digest(v.0.6.39), colorspace(v.2.1-1), textshaping(v.1.0.5), labeling(v.0.4.3), clusterGeneration(v.1.3.8), timechange(v.0.4.0), httr(v.1.4.8), abind(v.1.4-8), compiler(v.4.5.3), fontquiver(v.0.2.1), withr(v.3.0.2), doParallel(v.1.0.17), S7(v.0.2.2), backports(v.1.5.1), inline(v.0.3.21), optimParallel(v.1.0-2), QuickJSR(v.1.9.0), pkgbuild(v.1.4.8), MASS(v.7.3-65), openssl(v.2.4.0), scatterplot3d(v.0.3-44), loo(v.2.9.0), tools(v.4.5.3), rncl(v.0.8.9), otel(v.0.2.0), zip(v.2.3.3), rentrez(v.1.2.4), glue(v.1.8.1), quadprog(v.1.5-8), nlme(v.3.1-168), grid(v.4.5.3), checkmate(v.2.3.4), reshape2(v.1.4.5), generics(v.0.1.4), gtable(v.0.3.6), tzdb(v.0.5.0), data.table(v.1.18.2.1), hms(v.1.1.4), xml2(v.1.5.2), foreach(v.1.5.2), pillar(v.1.11.1), posterior(v.1.6.1), splines(v.4.5.3), lattice(v.0.22-9), survival(v.3.8-6), tidyselect(v.1.2.1), fontLiberation(v.0.1.0), knitr(v.1.51), fontBitstreamVera(v.0.1.1), gridExtra(v.2.3), V8(v.6.0.4), stats4(v.4.5.3), xfun(v.0.57), expm(v.1.0-0), bridgesampling(v.1.2-1), matrixStats(v.1.5.0), rstan(v.2.32.7), stringi(v.1.8.7), yaml(v.2.3.12), pacman(v.0.5.1), evaluate(v.1.0.5), codetools(v.0.2-20), officer(v.0.6.10), gdtools(v.0.4.2), cli(v.3.6.6), RcppParallel(v.5.1.11-2), xtable(v.1.8-8), systemfonts(v.1.3.2), coda(v.0.19-4.1), XML(v.3.99-0.22), parallel(v.4.5.3), rstantools(v.2.6.0), prettyunits(v.1.2.0), Brobdingnag(v.1.2-9), phangorn(v.2.12.1), mvtnorm(v.1.3-7), scales(v.1.4.0), crayon(v.1.5.3), combinat(v.0.0-8), rlang(v.1.2.0), fastmatch(v.1.1-8), multcomp(v.1.4-28) and mnormt(v.2.1.2)