Lab 2: Using Stan via brms

Applied Bayesian Modeling (ICPSR Summer Program 2025)

Author

Affiliation

Ursinus College

Published

July 20, 2025

The purpose of this tutorial is to show a complete workflow for estimating Bayesian models in R using the brms package Bürkner (2021) as an interface to Stan, as an alternative to rstanarm. The main difference between these two packages is that brms compiles each Stan model for each model fit, which will take a few extra seconds compared to models estimated in rstanarm. However, brms models may offer more flexibility, as you will see later in this course.

If you are new to R, please refer to the tutorial for Lab 1 first.
For suggestions for model presentation, processing MCMC output, or using Stan directly, stay tuned for the remainder of this course.

If you find any errors in this document or have suggestions for improvement, please email me.

How to use this tutorial

This tutorial focuses on using brms for fitting Bayesian models via R.

There are other options (foremost using Stan directly, but also JAGS and WinBUGS) for fitting Bayesian models that we will briefly discuss during the workshop. The brms package is similar in its ease of use and particularly useful for multilevel models. You can find more information about all these packages at the on the website for my ICPSR workshop. Some sections on that website are relevant for Mac or Windows users only as indicated.

Note

Refer to the rstanarm version of this tutorial for a more detailed explanation of some of the steps below.

What is not in this tutorial

This tutorial does not address the following topics, which will instead show up in future labs:

Convergence diagnostics (\(\rightarrow\) Lab 3)
Processing and working with MCMC output (\(\rightarrow\) Lab 4)
Writing a customized model directly in Stan (\(\rightarrow\) Lab 5)
Model presentation (\(\rightarrow\) Lab 6)

Using R as frontend

A convenient way to fit Bayesian models using Stan (or JAGS or WinBUGS or OpenBUGS) is to use R packages that function as frontends for Stan. These packages make it easy to do all your Bayesian data analysis in R, including:

importing and preparing the data
writing the empirical model
estimate the model using MCMC
process the output of Bayesian models
present output in publication-ready form

brms allows you to estimate Bayesian models almost without any additional steps compared to estimating standard frequentist models in R.

Installing Stan and brms

Please refer to Lab 0 for instructions to install the software needed for our course.

Using brms

The vignette for the brms package is exceptionally detailed and well structured, so there is no need to reproduce it here. Please read it. Below, I show what I consider the key steps for fitting a Bayesian linear regression model and producing output and convergence diagnostics.

Preparing the data and model

brms uses rectangular data frames as input, as do most other modeling functions in R (but not other Bayesian tools such as Stan or JAGS). For an example dataset, we simulate our own data in R. For this tutorial, we aim to fit a linear model, so we create a continuous outcome variable \(y\) as a function of two predictors \(x_1\) and \(x_2\) and a disturbance term \(e\). We simulate a dataset with 100 observations.

First, we create the predictors:

Code

n.sim <- 100; set.seed(123)
x1 <- rnorm(n = n.sim, mean = 5, sd = 2)
x2 <- rbinom(n.sim, size = 1, prob = 0.3)
e <- rnorm(n = n.sim, mean = 0, sd = 2)

Next, we create the outcome \(y\) based on coefficients \(b_1\) and \(b_2\) for the respective predictors and an intercept \(a\):

Code

b1 <- 1.2
b2 <- -3.1
a <- 1.5
y <- a + b1 * x1 + b2 * x2 + e

Now, we combine the variables into one data frame for processing later:

Code

sim.dat <- data.frame(y, x1, x2)

And we create and summarize a (frequentist) linear model (lm)¹ fit on these data:

¹ This is equivalent to fitting a generalized linear model using glm with the Gaussian (normal) link function; see more below.

Code

freq.mod <- lm(y ~ x1 + x2, 
               data = sim.dat)
summary(freq.mod)


Call:
lm(formula = y ~ x1 + x2, data = sim.dat)

Residuals:
   Min     1Q Median     3Q    Max 
-2.686 -1.359 -0.222  1.073  6.461 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    2.199      0.576    3.82  0.00024 ***
x1             1.070      0.103   10.37  < 2e-16 ***
x2            -3.087      0.413   -7.48  3.4e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.87 on 97 degrees of freedom
Multiple R-squared:  0.622, Adjusted R-squared:  0.614 
F-statistic: 79.9 on 2 and 97 DF,  p-value: <2e-16

Fitting the Bayesian model with all default arguments

Now, we create a Bayesian linear model using brms Conveniently, the relevant function is simply called brm(). This function allows fitting the model by specifying the prior distribution for coefficient estimates; all other arguments have default values (that we will modify later). You can find more about how to specify priors under ?brms::prior. For this first take, we use default priors and do not specify them.

Code

library("brms")
bayes.mod <- brm(y ~ x1 + x2, 
                     data = sim.dat,
                     family = gaussian()
                 )

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
using SDK: ‘MacOSX15.5.sdk’
clang -arch arm64 -std=gnu2x -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1


SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 3.3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 0.017 seconds (Warm-up)
Chain 1:                0.012 seconds (Sampling)
Chain 1:                0.029 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 2e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 0.017 seconds (Warm-up)
Chain 2:                0.013 seconds (Sampling)
Chain 2:                0.03 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 5e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 0.018 seconds (Warm-up)
Chain 3:                0.013 seconds (Sampling)
Chain 3:                0.031 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 2e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.02 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 0.016 seconds (Warm-up)
Chain 4:                0.013 seconds (Sampling)
Chain 4:                0.029 seconds (Total)
Chain 4:

Code

summary(bayes.mod)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y ~ x1 + x2 
   Data: sim.dat (Number of observations: 100) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     2.20      0.58     1.05     3.31 1.00     4197     3328
x1            1.07      0.11     0.87     1.28 1.00     4312     3197
x2           -3.07      0.42    -3.92    -2.26 1.00     4354     2872

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.90      0.14     1.65     2.20 1.00     4132     2744

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

You have now estimated your first Bayesian regression model! Note some more information in this output:

The posterior distribution contains 4000 draws - four chains with 1000 draws each.
The regression output contains quantiles of the posterior distribution for each parameter.
In addition to coefficients, the output also contains the posterior distribution for sigma (the (estimated) standard deviation of the outcome variable y)

In addition to this barebones call for brm, you can (and should) explicitly specify a few more arguments for full transparency and reproducibility. A call to ?brm and ?stan will reveal some of the additional arguments that you can specify.

prior: the prior distribution for the other regression coefficient(s). You can specify these using the set_prior function.
chains is the number of Markov chains to be used; the default is currently 4.
iter is the number of iterations for each chain (including warmup). The default is currently 2000.
warmup is the number of warmup (aka burnin) iterations per chain. The default is currently 50% of iter, i.e. 1000 if iter is set to its default.
thin specifies the which iterations should be saved; 1 (the default) saves each iteration, 5 the draw from every 5th iteration, etc.
init specifies the starting values.
- The default is set to use random starting values, generated by Stan: init = "random". This is not always preferable, especially in more complex models.
- Instead, you can give specific starting values for each parameter by providing a list or by providing the same starting value for all parameters: init = 1.
seed is the seed for random number generation; providing this explicitly is useful for reproducible estimates.
cores allows you to use multiple CPU cores on your computer, for parallel processing and therefore faster estimation.
backend allows you to choose between RStan and CmdStanR as backend for estimation. CmdStanR can be much faster; see Lab 5 for more information on CmdStanR.

With this information, you can re-specify the linear model above with explicit arguments:

Code

bayes.priors <- prior(normal(0, 0.1), class = "coef")
bayes.mod <- brm(y ~ x1 + x2, 
                     data = sim.dat,
                     family = gaussian(),
                     chains = 4,
                     iter = 5000,
                     warmup = 2500,
                     thin = 1,
                     init = 1,
                     seed = 123,
                     cores = 4,
                     backend = "cmdstanr")

Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 5000 [  0%]  (Warmup) 
Chain 1 Iteration:  100 / 5000 [  2%]  (Warmup) 
Chain 1 Iteration:  200 / 5000 [  4%]  (Warmup) 
Chain 1 Iteration:  300 / 5000 [  6%]  (Warmup) 
Chain 1 Iteration:  400 / 5000 [  8%]  (Warmup) 
Chain 1 Iteration:  500 / 5000 [ 10%]  (Warmup) 
Chain 1 Iteration:  600 / 5000 [ 12%]  (Warmup) 
Chain 1 Iteration:  700 / 5000 [ 14%]  (Warmup) 
Chain 1 Iteration:  800 / 5000 [ 16%]  (Warmup) 
Chain 1 Iteration:  900 / 5000 [ 18%]  (Warmup) 
Chain 1 Iteration: 1000 / 5000 [ 20%]  (Warmup) 
Chain 1 Iteration: 1100 / 5000 [ 22%]  (Warmup) 
Chain 1 Iteration: 1200 / 5000 [ 24%]  (Warmup) 
Chain 1 Iteration: 1300 / 5000 [ 26%]  (Warmup) 
Chain 1 Iteration: 1400 / 5000 [ 28%]  (Warmup) 
Chain 1 Iteration: 1500 / 5000 [ 30%]  (Warmup) 
Chain 1 Iteration: 1600 / 5000 [ 32%]  (Warmup) 
Chain 1 Iteration: 1700 / 5000 [ 34%]  (Warmup) 
Chain 1 Iteration: 1800 / 5000 [ 36%]  (Warmup) 
Chain 1 Iteration: 1900 / 5000 [ 38%]  (Warmup) 
Chain 1 Iteration: 2000 / 5000 [ 40%]  (Warmup) 
Chain 1 Iteration: 2100 / 5000 [ 42%]  (Warmup) 
Chain 1 Iteration: 2200 / 5000 [ 44%]  (Warmup) 
Chain 1 Iteration: 2300 / 5000 [ 46%]  (Warmup) 
Chain 1 Iteration: 2400 / 5000 [ 48%]  (Warmup) 
Chain 1 Iteration: 2500 / 5000 [ 50%]  (Warmup) 
Chain 1 Iteration: 2501 / 5000 [ 50%]  (Sampling) 
Chain 1 Iteration: 2600 / 5000 [ 52%]  (Sampling) 
Chain 1 Iteration: 2700 / 5000 [ 54%]  (Sampling) 
Chain 1 Iteration: 2800 / 5000 [ 56%]  (Sampling) 
Chain 1 Iteration: 2900 / 5000 [ 58%]  (Sampling) 
Chain 1 Iteration: 3000 / 5000 [ 60%]  (Sampling) 
Chain 1 Iteration: 3100 / 5000 [ 62%]  (Sampling) 
Chain 1 Iteration: 3200 / 5000 [ 64%]  (Sampling) 
Chain 1 Iteration: 3300 / 5000 [ 66%]  (Sampling) 
Chain 1 Iteration: 3400 / 5000 [ 68%]  (Sampling) 
Chain 1 Iteration: 3500 / 5000 [ 70%]  (Sampling) 
Chain 1 Iteration: 3600 / 5000 [ 72%]  (Sampling) 
Chain 1 Iteration: 3700 / 5000 [ 74%]  (Sampling) 
Chain 1 Iteration: 3800 / 5000 [ 76%]  (Sampling) 
Chain 1 Iteration: 3900 / 5000 [ 78%]  (Sampling) 
Chain 1 Iteration: 4000 / 5000 [ 80%]  (Sampling) 
Chain 1 Iteration: 4100 / 5000 [ 82%]  (Sampling) 
Chain 1 Iteration: 4200 / 5000 [ 84%]  (Sampling) 
Chain 1 Iteration: 4300 / 5000 [ 86%]  (Sampling) 
Chain 1 Iteration: 4400 / 5000 [ 88%]  (Sampling) 
Chain 1 Iteration: 4500 / 5000 [ 90%]  (Sampling) 
Chain 1 Iteration: 4600 / 5000 [ 92%]  (Sampling) 
Chain 1 Iteration: 4700 / 5000 [ 94%]  (Sampling) 
Chain 1 Iteration: 4800 / 5000 [ 96%]  (Sampling) 
Chain 1 Iteration: 4900 / 5000 [ 98%]  (Sampling) 
Chain 1 Iteration: 5000 / 5000 [100%]  (Sampling) 
Chain 2 Iteration:    1 / 5000 [  0%]  (Warmup) 
Chain 2 Iteration:  100 / 5000 [  2%]  (Warmup) 
Chain 2 Iteration:  200 / 5000 [  4%]  (Warmup) 
Chain 2 Iteration:  300 / 5000 [  6%]  (Warmup) 
Chain 2 Iteration:  400 / 5000 [  8%]  (Warmup) 
Chain 2 Iteration:  500 / 5000 [ 10%]  (Warmup) 
Chain 2 Iteration:  600 / 5000 [ 12%]  (Warmup) 
Chain 2 Iteration:  700 / 5000 [ 14%]  (Warmup) 
Chain 2 Iteration:  800 / 5000 [ 16%]  (Warmup) 
Chain 2 Iteration:  900 / 5000 [ 18%]  (Warmup) 
Chain 2 Iteration: 1000 / 5000 [ 20%]  (Warmup) 
Chain 2 Iteration: 1100 / 5000 [ 22%]  (Warmup) 
Chain 2 Iteration: 1200 / 5000 [ 24%]  (Warmup) 
Chain 2 Iteration: 1300 / 5000 [ 26%]  (Warmup) 
Chain 2 Iteration: 1400 / 5000 [ 28%]  (Warmup) 
Chain 2 Iteration: 1500 / 5000 [ 30%]  (Warmup) 
Chain 2 Iteration: 1600 / 5000 [ 32%]  (Warmup) 
Chain 2 Iteration: 1700 / 5000 [ 34%]  (Warmup) 
Chain 2 Iteration: 1800 / 5000 [ 36%]  (Warmup) 
Chain 2 Iteration: 1900 / 5000 [ 38%]  (Warmup) 
Chain 2 Iteration: 2000 / 5000 [ 40%]  (Warmup) 
Chain 2 Iteration: 2100 / 5000 [ 42%]  (Warmup) 
Chain 2 Iteration: 2200 / 5000 [ 44%]  (Warmup) 
Chain 2 Iteration: 2300 / 5000 [ 46%]  (Warmup) 
Chain 2 Iteration: 2400 / 5000 [ 48%]  (Warmup) 
Chain 2 Iteration: 2500 / 5000 [ 50%]  (Warmup) 
Chain 2 Iteration: 2501 / 5000 [ 50%]  (Sampling) 
Chain 2 Iteration: 2600 / 5000 [ 52%]  (Sampling) 
Chain 2 Iteration: 2700 / 5000 [ 54%]  (Sampling) 
Chain 2 Iteration: 2800 / 5000 [ 56%]  (Sampling) 
Chain 2 Iteration: 2900 / 5000 [ 58%]  (Sampling) 
Chain 2 Iteration: 3000 / 5000 [ 60%]  (Sampling) 
Chain 2 Iteration: 3100 / 5000 [ 62%]  (Sampling) 
Chain 2 Iteration: 3200 / 5000 [ 64%]  (Sampling) 
Chain 2 Iteration: 3300 / 5000 [ 66%]  (Sampling) 
Chain 2 Iteration: 3400 / 5000 [ 68%]  (Sampling) 
Chain 2 Iteration: 3500 / 5000 [ 70%]  (Sampling) 
Chain 2 Iteration: 3600 / 5000 [ 72%]  (Sampling) 
Chain 2 Iteration: 3700 / 5000 [ 74%]  (Sampling) 
Chain 2 Iteration: 3800 / 5000 [ 76%]  (Sampling) 
Chain 2 Iteration: 3900 / 5000 [ 78%]  (Sampling) 
Chain 2 Iteration: 4000 / 5000 [ 80%]  (Sampling) 
Chain 2 Iteration: 4100 / 5000 [ 82%]  (Sampling) 
Chain 2 Iteration: 4200 / 5000 [ 84%]  (Sampling) 
Chain 2 Iteration: 4300 / 5000 [ 86%]  (Sampling) 
Chain 2 Iteration: 4400 / 5000 [ 88%]  (Sampling) 
Chain 2 Iteration: 4500 / 5000 [ 90%]  (Sampling) 
Chain 2 Iteration: 4600 / 5000 [ 92%]  (Sampling) 
Chain 2 Iteration: 4700 / 5000 [ 94%]  (Sampling) 
Chain 2 Iteration: 4800 / 5000 [ 96%]  (Sampling) 
Chain 2 Iteration: 4900 / 5000 [ 98%]  (Sampling) 
Chain 2 Iteration: 5000 / 5000 [100%]  (Sampling) 
Chain 3 Iteration:    1 / 5000 [  0%]  (Warmup) 
Chain 3 Iteration:  100 / 5000 [  2%]  (Warmup) 
Chain 3 Iteration:  200 / 5000 [  4%]  (Warmup) 
Chain 3 Iteration:  300 / 5000 [  6%]  (Warmup) 
Chain 3 Iteration:  400 / 5000 [  8%]  (Warmup) 
Chain 3 Iteration:  500 / 5000 [ 10%]  (Warmup) 
Chain 3 Iteration:  600 / 5000 [ 12%]  (Warmup) 
Chain 3 Iteration:  700 / 5000 [ 14%]  (Warmup) 
Chain 3 Iteration:  800 / 5000 [ 16%]  (Warmup) 
Chain 3 Iteration:  900 / 5000 [ 18%]  (Warmup) 
Chain 3 Iteration: 1000 / 5000 [ 20%]  (Warmup) 
Chain 3 Iteration: 1100 / 5000 [ 22%]  (Warmup) 
Chain 3 Iteration: 1200 / 5000 [ 24%]  (Warmup) 
Chain 3 Iteration: 1300 / 5000 [ 26%]  (Warmup) 
Chain 3 Iteration: 1400 / 5000 [ 28%]  (Warmup) 
Chain 3 Iteration: 1500 / 5000 [ 30%]  (Warmup) 
Chain 3 Iteration: 1600 / 5000 [ 32%]  (Warmup) 
Chain 3 Iteration: 1700 / 5000 [ 34%]  (Warmup) 
Chain 3 Iteration: 1800 / 5000 [ 36%]  (Warmup) 
Chain 3 Iteration: 1900 / 5000 [ 38%]  (Warmup) 
Chain 3 Iteration: 2000 / 5000 [ 40%]  (Warmup) 
Chain 3 Iteration: 2100 / 5000 [ 42%]  (Warmup) 
Chain 3 Iteration: 2200 / 5000 [ 44%]  (Warmup) 
Chain 3 Iteration: 2300 / 5000 [ 46%]  (Warmup) 
Chain 3 Iteration: 2400 / 5000 [ 48%]  (Warmup) 
Chain 3 Iteration: 2500 / 5000 [ 50%]  (Warmup) 
Chain 3 Iteration: 2501 / 5000 [ 50%]  (Sampling) 
Chain 3 Iteration: 2600 / 5000 [ 52%]  (Sampling) 
Chain 3 Iteration: 2700 / 5000 [ 54%]  (Sampling) 
Chain 3 Iteration: 2800 / 5000 [ 56%]  (Sampling) 
Chain 3 Iteration: 2900 / 5000 [ 58%]  (Sampling) 
Chain 3 Iteration: 3000 / 5000 [ 60%]  (Sampling) 
Chain 3 Iteration: 3100 / 5000 [ 62%]  (Sampling) 
Chain 3 Iteration: 3200 / 5000 [ 64%]  (Sampling) 
Chain 3 Iteration: 3300 / 5000 [ 66%]  (Sampling) 
Chain 3 Iteration: 3400 / 5000 [ 68%]  (Sampling) 
Chain 3 Iteration: 3500 / 5000 [ 70%]  (Sampling) 
Chain 3 Iteration: 3600 / 5000 [ 72%]  (Sampling) 
Chain 3 Iteration: 3700 / 5000 [ 74%]  (Sampling) 
Chain 3 Iteration: 3800 / 5000 [ 76%]  (Sampling) 
Chain 3 Iteration: 3900 / 5000 [ 78%]  (Sampling) 
Chain 3 Iteration: 4000 / 5000 [ 80%]  (Sampling) 
Chain 3 Iteration: 4100 / 5000 [ 82%]  (Sampling) 
Chain 3 Iteration: 4200 / 5000 [ 84%]  (Sampling) 
Chain 3 Iteration: 4300 / 5000 [ 86%]  (Sampling) 
Chain 3 Iteration: 4400 / 5000 [ 88%]  (Sampling) 
Chain 3 Iteration: 4500 / 5000 [ 90%]  (Sampling) 
Chain 3 Iteration: 4600 / 5000 [ 92%]  (Sampling) 
Chain 3 Iteration: 4700 / 5000 [ 94%]  (Sampling) 
Chain 3 Iteration: 4800 / 5000 [ 96%]  (Sampling) 
Chain 3 Iteration: 4900 / 5000 [ 98%]  (Sampling) 
Chain 3 Iteration: 5000 / 5000 [100%]  (Sampling) 
Chain 4 Iteration:    1 / 5000 [  0%]  (Warmup) 
Chain 4 Iteration:  100 / 5000 [  2%]  (Warmup) 
Chain 4 Iteration:  200 / 5000 [  4%]  (Warmup) 
Chain 4 Iteration:  300 / 5000 [  6%]  (Warmup) 
Chain 4 Iteration:  400 / 5000 [  8%]  (Warmup) 
Chain 4 Iteration:  500 / 5000 [ 10%]  (Warmup) 
Chain 4 Iteration:  600 / 5000 [ 12%]  (Warmup) 
Chain 4 Iteration:  700 / 5000 [ 14%]  (Warmup) 
Chain 4 Iteration:  800 / 5000 [ 16%]  (Warmup) 
Chain 4 Iteration:  900 / 5000 [ 18%]  (Warmup) 
Chain 4 Iteration: 1000 / 5000 [ 20%]  (Warmup) 
Chain 4 Iteration: 1100 / 5000 [ 22%]  (Warmup) 
Chain 4 Iteration: 1200 / 5000 [ 24%]  (Warmup) 
Chain 4 Iteration: 1300 / 5000 [ 26%]  (Warmup) 
Chain 4 Iteration: 1400 / 5000 [ 28%]  (Warmup) 
Chain 4 Iteration: 1500 / 5000 [ 30%]  (Warmup) 
Chain 4 Iteration: 1600 / 5000 [ 32%]  (Warmup) 
Chain 4 Iteration: 1700 / 5000 [ 34%]  (Warmup) 
Chain 4 Iteration: 1800 / 5000 [ 36%]  (Warmup) 
Chain 4 Iteration: 1900 / 5000 [ 38%]  (Warmup) 
Chain 4 Iteration: 2000 / 5000 [ 40%]  (Warmup) 
Chain 4 Iteration: 2100 / 5000 [ 42%]  (Warmup) 
Chain 4 Iteration: 2200 / 5000 [ 44%]  (Warmup) 
Chain 4 Iteration: 2300 / 5000 [ 46%]  (Warmup) 
Chain 4 Iteration: 2400 / 5000 [ 48%]  (Warmup) 
Chain 4 Iteration: 2500 / 5000 [ 50%]  (Warmup) 
Chain 4 Iteration: 2501 / 5000 [ 50%]  (Sampling) 
Chain 4 Iteration: 2600 / 5000 [ 52%]  (Sampling) 
Chain 4 Iteration: 2700 / 5000 [ 54%]  (Sampling) 
Chain 4 Iteration: 2800 / 5000 [ 56%]  (Sampling) 
Chain 4 Iteration: 2900 / 5000 [ 58%]  (Sampling) 
Chain 4 Iteration: 3000 / 5000 [ 60%]  (Sampling) 
Chain 4 Iteration: 3100 / 5000 [ 62%]  (Sampling) 
Chain 4 Iteration: 3200 / 5000 [ 64%]  (Sampling) 
Chain 4 Iteration: 3300 / 5000 [ 66%]  (Sampling) 
Chain 4 Iteration: 3400 / 5000 [ 68%]  (Sampling) 
Chain 4 Iteration: 3500 / 5000 [ 70%]  (Sampling) 
Chain 4 Iteration: 3600 / 5000 [ 72%]  (Sampling) 
Chain 4 Iteration: 3700 / 5000 [ 74%]  (Sampling) 
Chain 4 Iteration: 3800 / 5000 [ 76%]  (Sampling) 
Chain 4 Iteration: 3900 / 5000 [ 78%]  (Sampling) 
Chain 4 Iteration: 4000 / 5000 [ 80%]  (Sampling) 
Chain 4 Iteration: 4100 / 5000 [ 82%]  (Sampling) 
Chain 4 Iteration: 4200 / 5000 [ 84%]  (Sampling) 
Chain 4 Iteration: 4300 / 5000 [ 86%]  (Sampling) 
Chain 4 Iteration: 4400 / 5000 [ 88%]  (Sampling) 
Chain 4 Iteration: 4500 / 5000 [ 90%]  (Sampling) 
Chain 4 Iteration: 4600 / 5000 [ 92%]  (Sampling) 
Chain 4 Iteration: 4700 / 5000 [ 94%]  (Sampling) 
Chain 4 Iteration: 4800 / 5000 [ 96%]  (Sampling) 
Chain 4 Iteration: 4900 / 5000 [ 98%]  (Sampling) 
Chain 4 Iteration: 5000 / 5000 [100%]  (Sampling) 
Chain 1 finished in 0.1 seconds.
Chain 2 finished in 0.1 seconds.
Chain 3 finished in 0.1 seconds.
Chain 4 finished in 0.1 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.1 seconds.
Total execution time: 0.3 seconds.

Code

summary(bayes.mod)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: y ~ x1 + x2 
   Data: sim.dat (Number of observations: 100) 
  Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
         total post-warmup draws = 10000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     2.21      0.57     1.06     3.34 1.00    11220     7552
x1            1.07      0.10     0.87     1.27 1.00    10769     7646
x2           -3.09      0.42    -3.91    -2.27 1.00     9978     7189

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     1.89      0.14     1.65     2.18 1.00    10744     7563

Draws were sampled using sample(hmc). 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).

Creating a regression table

To summarize the key results from the model, you will usually create a table showing estimates for each parameter and a measure of uncertainty around them. With a few colleagues, I’ve built the R package “BayesPostEst” that contains two functions doing just that: creating a regression table for you.² You have to load the package first (since you already installed it at the beginning of our workshop, no need to install it again)

² For more information, see the package website or the companion article in the Journal of Open Source Software.

Code

library("BayesPostEst")

Now, you can use mcmcTab() to create your summary table:

Code

mcmcTab(bayes.mod)

     Variable  Median    SD   Lower   Upper
1 b_Intercept    2.21 0.574    1.05    3.34
2        b_x1    1.07 0.103    0.87    1.27
3        b_x2   -3.09 0.417   -3.91   -2.27
4       sigma    1.89 0.138    1.65    2.18
5   Intercept    6.85 0.190    6.48    7.22
6      lprior   -3.76 0.036   -3.84   -3.70
7        lp__ -207.99 1.439 -212.03 -206.52

You can also select specific parameters of interest only:

Code

mcmcTab(bayes.mod, pars = c("b_x1", "b_x2"))

  Variable Median    SD Lower Upper
1     b_x1   1.07 0.103  0.87  1.27
2     b_x2  -3.09 0.417 -3.91 -2.27

If you wish to print only some of the columns, you can select these, too:

Code

mcmcTab(bayes.mod, pars = c("b_x1", "b_x2"))[, c("Median", "SD")]

  Median    SD
1   1.07 0.103
2  -3.09 0.417

And you can export this table to RMarkdown using the knitr::kable() function:

Code

bayes.tab <- mcmcTab(bayes.mod, pars = c("b_x1", "b_x2"))
library("knitr")
kable(bayes.tab, 
      digits = 2,
      caption = "Summary of posterior estimates.",
      col.names = c("Variable", "Median", "Std. Dev.", "Lower 95% CI", "Upper 95% CI"))

Summary of posterior estimates.
Variable	Median	Std. Dev.	Lower 95% CI	Upper 95% CI
b_x1	1.07	0.10	0.87	1.27
b_x2	-3.09	0.42	-3.91	-2.27

An alternative that’s more useful for displaying multiple models and conforming to journal style requirements is the mcmcReg() function, also from the BayesPostEst package. It works similar to the texreg() function from the package with the same name:

Code

mcmcReg(mod = bayes.mod,
        pars = c("b_x1", "b_x2"))


\begin{table}
\begin{center}
\begin{tabular}{l c}
\hline
 & Model 1 \\
\hline
b\_x1 & $1.07^{*}$        \\
      & $ [ 0.87;  1.27]$ \\
b\_x2 & $-3.09^{*}$       \\
      & $ [-3.91; -2.27]$ \\
\hline
\multicolumn{2}{l}{\scriptsize{$^*$ 0 outside 95\% credible interval.}}
\end{tabular}
\caption{Statistical models}
\label{table:coefficients}
\end{center}
\end{table}

mcmcReg() produces a table for use in a LaTeX file by default, but can also write an HTML file that can be opened and processed in Word and similar word processors.

Code

mcmcReg(mod = bayes.mod, 
        pars = c("b_x1", "b_x2"),
        format = "html", 
        file = "bayes_m1_table2.docx")

The HTML output can also be passed on directly to the RMarkdown document:

Code

mcmcReg(mod = bayes.mod, 
        pars = c("b_x1", "b_x2"),
        format = "html")

Statistical models
	Model 1
b_x1	1.07^*
	[ 0.87; 1.27]
b_x2	-3.09^*
	[-3.91; -2.27]
^* 0 outside 95% credible interval.

For the many ways to customize a call to mcmcReg(), see the help file. This function is under active development and will likely see some enhancements soon. Leave us feature requests on Github - just open an issue!³

³ You can find the Github site for opening issues here: https://github.com/ShanaScogin/BayesPostEst/issues.

A more generalized package for tables from many regression models (Bayesian or frequentist) is the amazing “modelsummary” package. It will create tables for brms objects as well:

Code

library("modelsummary")
modelsummary(list(bayes.mod))

	(1)
b_Intercept	2.212
b_x1	1.070
b_x2	-3.086
sigma	1.885
Num.Obs.	100
R2	0.621
R2 Adj.	0.608
ELPD	-207.1
ELPD s.e.	8.1
LOOIC	414.3
LOOIC s.e.	16.2
WAIC	414.3
RMSE	1.85

but may require a bit of customization to create simple, straightforward tables.

Visualizing results directly from the brms object

brms and rstan have a convenient suite of commands to summarize and visualize the results from a fitted model. You already saw the results from summary() above. In addition, you can plot a variety of quantities using mcmc_plot (see ?mcmc_plot for more information).

A call to plot without arguments will produce a coefficient dotplot for all parameters:

Code

mcmc_plot(bayes.mod)

You can restrict the plot to parameters of interest using the pars argument:

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"))

In addition, you can produce a multitude of plots to visualize the posterior distribution, including:

Autocorrelation plots for each parameter

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = c("acf"))

Histograms

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "hist")

Density plots by parameter

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "dens")

Density plots distinguishing between chains

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "dens_overlay")

Violin plots distinguishing between chains

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "violin")

Density plots added to a coefficient dot plot setup

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "areas_ridges")

Density plots added to a coefficient dot plot setup, with credible intervals highlighted

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "areas")

As with any use of ggplot2, you can create an object and then save the plot to your folder:

Code

library("ggplot2")
p <- mcmc_plot(bayes.mod, 
     pars = c("x1", "x2"),
     type = "areas_ridges")
ggsave(p, file = "bayes_m1_ridgeplot.pdf", width = 5, height = 2.5)

And because the plotting function is based on ggplot2, you can customize the appearance of the plots:

Code

mcmc_plot(bayes.mod, 
     pars = c("x1", "x2")) + 
  theme_bw() + 
  xlab("Posterior estimate") + 
  labs(title = "Regression results", caption = "N = 100; estimates based on 4 chains with 2500 draws each.")

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.

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