Bayesian multilevel mediation with brms

library(knitr)
library(bmlm)
library(tinytable)
library(ggdist)
library(posterior)
library(patchwork)
library(brms)
library(tidyverse)

dir.create("cache", FALSE)
options(
  brms.backend = Sys.getenv("BRMS_BACKEND", "rstan"),
  brms.threads = as.numeric(Sys.getenv("BRMS_THREADS", 1)),
  mc.cores = as.numeric(Sys.getenv("MAX_CORES", 4)),
  brms.short_summary = TRUE,
  tinytable_tt_digits = 2,
  tinytable_format_num_fmt = "decimal",
  tinytable_format_num_zero = TRUE,
  tinytable_tt_theme = "spacing"
)
theme_set(
  theme_linedraw(base_size = if_else(knitr::is_html_output(), 10, 12)) +
    theme(
      panel.grid = element_blank(),
      strip.background = element_blank(),
      strip.text = element_text(hjust = 0, colour = "black")
    )
)
# Omit MCMC info in brmsfit.summary
.summary <- function(x) {
  out <- summary(x)
  out$random$id <- out$random$id[, 1:4]
  out$fixed <- out$fixed[, 1:4]
  out$spec_pars <- out$spec_pars[, 1:4]
  out
}

bmlm

Back in 2016 I wrote an R package for bayesian estimation of a multivariate multilevel model for assessing a three-variable causal mediation model (Vuorre [2016] 2024). In the abstract to the article that discussed the methodology, Niall Bolger and I wrote

“Statistical mediation allows researchers to investigate potential causal effects of experimental manipulations through intervening variables. It is a powerful tool for assessing the presence and strength of postulated causal mechanisms. Although mediation is used in certain areas of psychology, it is rarely applied in cognitive psychology and neuroscience. One reason for the scarcity of applications is that these areas of psychology commonly employ within-subjects designs, and mediation models for within-subjects data is considerably more complicated than for between-subjects data. Here, we draw attention to the importance and ubiquity of mediational hypotheses in within-subjects designs, and we present a general and flexible software package for conducting Bayesian within-subjects mediation analyses in the R programming environment. We use experimental data from cognitive psychology to illustrate the benefits of within-subject mediation for theory testing and comparison.” (Vuorre and Bolger 2017)

I wrote the R package as an interface to a Stan (Team 2024) model because brms (Bürkner 2017a)—still in its early days—did not implement the kind of multivariate structure required by the model. Shortly afterwards, probably within a few months actually, Paul updated brms to fit the required model structure and bmlm as a standalone package lost much of its value. So whenever people email me about bmlm, I keep suggesting them to estimate their models with brms instead because it can do this, and so much more. (I’ve had a tutorial for this up at https://vuorre.com/brms-workshop/posts/mediation/ but it’s difficult to find.)

So in this post I’ll briefly show how to fit bmlm’s multilevel mediation model with brms, along with the required post-processing for computing the indirect effects, figures, etc.

Mediation—a word of caution

Mediation models are used to make causal claims from observational data. This is a complex and difficult endeavor, and all uses of mediation must appropriately wrestle with the implications and assumptions behind the models and claims that are being made. See Green, Ha, and Bullock (2010) and Rohrer et al. (2022).

Analysis

In one of bmlm’s vignettes we analyse an example data set from Intensive Longitudinal Methods: An Introduction to Diary and Experience Sampling Research (Bolger and Laurenceau 2013). The data, shown in Table 1 and included in the bmlm package, indicate several (hypothetical) participants’ (id) work stressors (fwkstrs), work dissatisfaction (fwkdis), and relationship dissatisfaction (freldis) over several days of study.

dat <- tibble(BLch9)[, c(1, 3:5)]
tt(head(dat))

Table 1: Six rows of example data.

id	fwkstrs	fwkdis	freldis
101	3	5.59	3.03
101	3	5.54	4.62
101	3	3.89	2.85
101	4	5.35	6.40
101	1	4.48	2.54
101	2	3.34	5.16

To see how the full analysis is conducted with bmlm, please see the vignette (https://vuorre.com/bmlm/articles/bmlm-blch9/bmlm-blch9.html). Here, we implement the analysis without the use of bmlm’s functions.

Data preparation

First, we must isolate the within-person deviations of the key variables. We consider work stressors to be the independent variable, work dissatisfaction the mediator, and relationship dissatisfaction the outcome variable, and label them accordingly as x, m, and y for brevity.

dat <- dat |>
  mutate(
    x = fwkstrs - mean(fwkstrs, na.rm = TRUE),
    m = fwkdis - mean(fwkdis, na.rm = TRUE),
    y = freldis - mean(freldis, na.rm = TRUE),
    .by = id,
    .keep = "unused"
  )
tt(head(dat))

Table 2: Subject-mean centered variables.

id	x	m	y
101	0.33	0.98	-1.44
101	0.33	0.93	0.14
101	0.33	-0.72	-1.63
101	1.33	0.75	1.92
101	-1.67	-0.12	-1.93
101	-0.67	-1.27	0.69

Model fitting

The model comprises of two regression formulas that share a variance-covariance matrix for the random effects (the | p | syntax). Notice that I have optimized my MCMC sampler options in an environment files and using options() as shown here. In addition to directly connecting this model to the underlying regressions, brms estimates the model faster than bmlm (and allows using better priors, omitted here).

path_m <- bf(
  m ~ x + (x | p | id)
)
path_y <- bf(
  y ~ x + m + (x + m | p | id)
)
fit <- brm(
  path_m + path_y + set_rescor(FALSE),
  data = dat,
  file = "cache/fit"
)

Model summary

The model parameters directly gives us the x -> m (m_x, called a in bmlm), x -> y (y_m, called c' or cp in bmlm) and m -> y (y_m, called b in bmlm) path coefficients. These precisely match bmlm’s estimates.

 Family: MV(gaussian, gaussian) 
  Links: mu = identity; sigma = identity
         mu = identity; sigma = identity 
Formula: m ~ x + (x | p | id) 
         y ~ x + m + (x + m | p | id) 
   Data: dat (Number of observations: 2100) 

Multilevel Hyperparameters:
~id (Number of levels: 100) 
                             Estimate Est.Error l-95% CI u-95% CI
sd(m_Intercept)                  0.02      0.01     0.00     0.05
sd(m_x)                          0.26      0.04     0.19     0.34
sd(y_Intercept)                  0.02      0.01     0.00     0.05
sd(y_x)                          0.08      0.03     0.01     0.14
sd(y_m)                          0.22      0.03     0.17     0.28
cor(m_Intercept,m_x)             0.02      0.39    -0.72     0.77
cor(m_Intercept,y_Intercept)     0.02      0.41    -0.75     0.77
cor(m_x,y_Intercept)            -0.01      0.41    -0.77     0.74
cor(m_Intercept,y_x)            -0.01      0.41    -0.77     0.75
cor(m_x,y_x)                     0.24      0.29    -0.39     0.77
cor(y_Intercept,y_x)            -0.01      0.42    -0.77     0.76
cor(m_Intercept,y_m)            -0.00      0.40    -0.75     0.72
cor(m_x,y_m)                     0.48      0.15     0.17     0.75
cor(y_Intercept,y_m)            -0.01      0.42    -0.79     0.76
cor(y_x,y_m)                     0.56      0.26    -0.11     0.91

Regression Coefficients:
            Estimate Est.Error l-95% CI u-95% CI
m_Intercept    -0.00      0.02    -0.05     0.05
y_Intercept    -0.00      0.02    -0.04     0.04
m_x             0.19      0.04     0.12     0.26
y_x             0.10      0.02     0.06     0.15
y_m             0.15      0.03     0.09     0.21

Further Distributional Parameters:
        Estimate Est.Error l-95% CI u-95% CI
sigma_m     1.09      0.02     1.06     1.13
sigma_y     0.93      0.01     0.90     0.96

Mediation parameters

To get the additional mediation metrics (the indirect effect, the total effect, and others), we simply wrangle the model’s posterior samples. The indirect effect, or “mediated effect” is \(me = ab + \sigma_{{a_j}{b_j}}\), or the population-level m_x times the population-level y_m plus the covariance of the person-level m_x and y_ms. The total effect is then \(c = me + c'\). The proportion of the total effect that is mediated is \(me / c\). While wrangling the variables below, I rename them for clarity.

draws <- as_draws_df(
  fit,
  variable = c(
    "b_m_x",
    "b_y_x",
    "b_y_m",
    "sd_id__m_x",
    "sd_id__y_m",
    "cor_id__m_x__y_m"
  )
) |>
  mutate(
    a = b_m_x,
    b = b_y_m,
    cp = b_y_x,
    covab = cor_id__m_x__y_m * sd_id__m_x * sd_id__y_m,
    me = a * b + covab,
    c = me + cp,
    pme = me / c,
    .keep = "unused"
  )

draws now contains the posterior draws of the key population-level effects, which we summarize below.

draws |>
  summarise_draws(
    mean,
    sd,
    ~ quantile2(.x, probs = c(.025, .975))
  ) |>
  tt()

Table 3: Summaries of key quantities posterior draws.

variable	mean	sd	q2.5	q97.5
a	0.19	0.04	0.12	0.26
b	0.15	0.03	0.09	0.21
cp	0.10	0.02	0.06	0.15
covab	0.03	0.01	0.01	0.05
me	0.06	0.01	0.03	0.09
c	0.16	0.03	0.11	0.21
pme	0.36	0.08	0.21	0.53

These parameter estimates directly reproduce the values obtained by bmlm.

Graphics

brms has a number of excellent built-in visualization facilities, such as drawing conditional effects with ggplot2 (Wickham 2016) (Figure 1).

ce1 <- conditional_effects(
  fit,
  effects = "x",
  resp = "m",
  robust = FALSE
)
ce1 <- plot(ce1, plot = FALSE)[[1]]
ce2 <- conditional_effects(
  fit,
  effects = "m",
  resp = "y",
  robust = FALSE
)
ce2 <- plot(ce2, plot = FALSE)[[1]]
ce1 | ce2

Figure 1: Conditional population-level regression of m on x (left) and y on m (right).

And the resulting objects and posterior samples are easily visualized with e.g. functions from the ggdist package (Kay 2024a) (Figure 2).

draws |>
  select(!starts_with(".")) |>
  pivot_longer(everything()) |>
  ggplot(aes(x = value, y = name)) +
  stat_slabinterval(normalize = "xy")

Figure 2: Approximate posterior distributions of key mediation parameters.

Conclusion

If you are going to fit the model described in (Vuorre and Bolger 2017), I recommend using the brms R package (Bürkner 2017a) because of its flexibility, how it clearly connects to R regression syntax, and its estimation efficiency.

Feedback & comments

I’d appreciate any feedback or comments you might have. Feel free to le me know what you think either using the comments field (below) or on Bluesky:

R environment

library(grateful)
cite_packages(output = "paragraph", pkgs = "Session", out.dir = getwd())

We used R version 4.5.0 (R Core Team 2025) and the following R packages: bmlm v. 1.3.15 (Vuorre and Bolger 2018; Vuorre 2023), brms v. 2.22.0 (Bürkner 2017b, 2018, 2021), distributional v. 0.5.0 (O’Hara-Wild et al. 2024), ggdist v. 3.3.2 (Kay 2024c, 2024b), ggnewscale v. 0.5.1 (Campitelli 2025), knitr v. 1.50 (Xie 2014, 2015, 2025), latex2exp v. 0.9.6 (Meschiari 2022), patchwork v. 1.3.0 (Pedersen 2024), posterior v. 1.6.1 (Vehtari et al. 2021, 2024; Lambert and Vehtari 2022; Margossian et al. 2024; Bürkner et al. 2025), Rcpp v. 1.0.14 (Eddelbuettel and François 2011; Eddelbuettel 2013; Eddelbuettel and Balamuta 2018; Eddelbuettel et al. 2025), scales v. 1.3.0 (Wickham, Pedersen, and Seidel 2023), tidybayes v. 3.0.7 (Kay 2024d), tidyverse v. 2.0.0 (Wickham et al. 2019), tinytable v. 0.8.0 (Arel-Bundock 2025).

References

Arel-Bundock, Vincent. 2025. tinytable: Simple and Configurable Tables in “HTML,” “LaTeX,” “Markdown,” “Word,” “PNG,” “PDF,” and “Typst” Formats. https://doi.org/10.32614/CRAN.package.tinytable.

Bolger, Niall, and Jean-Philippe Laurenceau. 2013. Intensive Longitudinal Methods: An Introduction to Diary and Experience Sampling Research. New York, NY: Guilford Press. http://www.intensivelongitudinal.com/.

Bürkner, Paul-Christian. 2017a. “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.

———. 2017b. “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.

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

———. 2021. “Bayesian Item Response Modeling in R with brms and Stan.” Journal of Statistical Software 100 (5): 1–54. https://doi.org/10.18637/jss.v100.i05.

Bürkner, Paul-Christian, Jonah Gabry, Matthew Kay, and Aki Vehtari. 2025. “posterior: Tools for Working with Posterior Distributions.” https://mc-stan.org/posterior/.

Campitelli, Elio. 2025. ggnewscale: Multiple Fill and Colour Scales in “ggplot2”. https://doi.org/10.32614/CRAN.package.ggnewscale.

Eddelbuettel, Dirk. 2013. Seamless R and C++ Integration with Rcpp. New York: Springer. https://doi.org/10.1007/978-1-4614-6868-4.

Eddelbuettel, Dirk, and James Joseph Balamuta. 2018. “Extending R with C++: A Brief Introduction to Rcpp.” The American Statistician 72 (1): 28–36. https://doi.org/10.1080/00031305.2017.1375990.

Eddelbuettel, Dirk, Romain Francois, JJ Allaire, Kevin Ushey, Qiang Kou, Nathan Russell, Iñaki Ucar, Doug Bates, and John Chambers. 2025. Rcpp: Seamless r and c++ Integration. https://doi.org/10.32614/CRAN.package.Rcpp.

Eddelbuettel, Dirk, and Romain François. 2011. “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software 40 (8): 1–18. https://doi.org/10.18637/jss.v040.i08.

Green, Donald P., Shang E. Ha, and John G. Bullock. 2010. “Enough Already about ‘Black Box’ Experiments: Studying Mediation Is More Difficult Than Most Scholars Suppose.” The ANNALS of the American Academy of Political and Social Science 628 (1): 200–208. https://doi.org/10.1177/0002716209351526.

Kay, Matthew. 2024a. ggdist: Visualizations of Distributions and Uncertainty. https://doi.org/10.5281/zenodo.3879620.

———. 2024b. ggdist: Visualizations of Distributions and Uncertainty. https://doi.org/10.5281/zenodo.3879620.

———. 2024c. “ggdist: Visualizations of Distributions and Uncertainty in the Grammar of Graphics.” IEEE Transactions on Visualization and Computer Graphics 30 (1): 414–24. https://doi.org/10.1109/TVCG.2023.3327195.

———. 2024d. tidybayes: Tidy Data and Geoms for Bayesian Models. https://doi.org/10.5281/zenodo.1308151.

Lambert, Ben, and Aki Vehtari. 2022. “Rstar: A Robust MCMC Convergence Diagnostic with Uncertainty Using Decision Tree Classifiers.” Bayesian Analysis 17 (2): 353–79. https://doi.org/10.1214/20-BA1252.

Margossian, Charles C., Matthew D. Hoffman, Pavel Sountsov, Lionel Riou-Durand, Aki Vehtari, and Andrew Gelman. 2024. “Nested Rhat: Assessing the Convergence of Markov Chain Monte Carlo When Running Many Short Chains.” Bayesian Analysis. https://doi.org/10.1214/24-BA1453.

Meschiari, Stefano. 2022. Latex2exp: Use LaTeX Expressions in Plots. https://doi.org/10.32614/CRAN.package.latex2exp.

O’Hara-Wild, Mitchell, Matthew Kay, Alex Hayes, and Rob Hyndman. 2024. distributional: Vectorised Probability Distributions. https://doi.org/10.32614/CRAN.package.distributional.

Pedersen, Thomas Lin. 2024. patchwork: The Composer of Plots. https://doi.org/10.32614/CRAN.package.patchwork.

R Core Team. 2025. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Rohrer, Julia M., Paul Hünermund, Ruben C. Arslan, and Malte Elson. 2022. “That’s a Lot to Process! Pitfalls of Popular Path Models.” Advances in Methods and Practices in Psychological Science 5 (2): 25152459221095827. https://doi.org/10.1177/25152459221095827.

Team, Stan Development. 2024. “Stan Modeling Language Users Guide and Reference Manual, Version 2.36.” https://mc-stan.org.

Vehtari, Aki, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. 2021. “Rank-Normalization, Folding, and Localization: An Improved Rhat for Assessing Convergence of MCMC (with Discussion).” Bayesian Analysis 16 (2): 667–718.

Vehtari, Aki, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. 2024. “Pareto Smoothed Importance Sampling.” Journal of Machine Learning Research 25 (72): 1–58.

Vuorre, Matti. 2023. bmlm: Bayesian Multilevel Mediation. https://CRAN.R-project.org/package=bmlm.

———. (2016) 2024. “Bmlm: Bayesian Multilevel Mediation.” https://github.com/mvuorre/bmlm.

Vuorre, Matti, and Niall Bolger. 2017. “Within-Subject Mediation Analysis for Experimental Data in Cognitive Psychology and Neuroscience.” Behavior Research Methods, December, 1–19. https://doi.org/10.3758/s13428-017-0980-9.

———. 2018. “Within-Subject Mediation Analysis for Experimental Data in Cognitive Psychology and Neuroscience.” Behavior Research Methods. https://doi.org/10.3758/s13428-017-0980-9.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.

Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2023. scales: Scale Functions for Visualization. https://doi.org/10.32614/CRAN.package.scales.

Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.

———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.

———. 2025. knitr: A General-Purpose Package for Dynamic Report Generation in R. https://yihui.org/knitr/.