Getting Started with brms

The brms package provides a powerful interface for fitting Bayesian regression models in R. Built on top of Stan, a state-of-the-art platform for Bayesian inference, brms allows you to specify complex statistical models using familiar R formula syntax. This tutorial walks you through fitting your first Bayesian models with brms.

What Is brms?

brms stands for Bayesian Regression Models using Stan. It translates R formula notation into Stan code, handles the compilation and sampling, and returns results in a format that integrates seamlessly with tidyverse workflows.

Key advantages of using brms:

• Familiar syntax: Use the same ~ formula notation as lm() and glm()
• Flexible models: Handle multilevel/hierarchical models, ordinal outcomes, censored data, and more
• Automatic MCMC: Stan handles Markov Chain Monte Carlo sampling under the hood
• Tidy output: Results play nicely with dplyr, tidybayes, and ggplot2

Installation

First, install brms and its dependencies:

# Install brms
install.packages("brms")

# Install Stan if not already installed
install.packages("rstan")

# Recommended packages for visualization
install.packages("bayesplot")
install.packages("tidybayes")

Loading the packages:

library(brms)
library(bayesplot)
library(tidybayes)
library(ggplot2)
library(dplyr)

Your First Bayesian Model

Let’s start with a simple linear regression. We’ll use the mtcars dataset to predict miles per gallon (mpg) from weight (wt) and number of cylinders (cyl):

# Fit a Bayesian linear regression
model1 <- brm(
  mpg ~ wt + cyl,
  data = mtcars,
  family = gaussian(),
  chains = 4,
  iter = 2000,
  seed = 123
)

Breaking down the arguments:

• mpg ~ wt + cyl: Formula specifying mpg as the outcome, with wt and cyl as predictors
• data: The data frame
• family: The probability distribution for the outcome (gaussian for continuous outcomes)
• chains: Number of Markov chains to run (4 is standard)
• iter: Total iterations per chain (including warmup)
• seed: For reproducibility

Examining Model Results

The summary() function displays the core results:

summary(model1)

Output includes:

• Population-level effects: Coefficients for your predictors
• Sigma: The residual standard deviation
• Rhat: Convergence diagnostic (should be close to 1.0)
• ESS: Effective sample size (higher is better)

For a more tidy output:

fixef(model1)

This shows the posterior estimates for each coefficient:

        Estimate Est.Error  Q2.5 Q97.5
Intercept  39.69      1.96 35.80 43.63
wt         -5.65      0.58 -6.78 -4.52
cyl        -1.51      0.78 -3.04 -0.02

Interpretation:

• Estimate: Posterior mean (the most likely value)
• Est.Error: Standard deviation of the posterior
• Q2.5, Q97.5: 95% credible interval bounds

Interpreting Credible Intervals

The 95% credible interval [-6.78, -4.52] for wt means there’s a 95% probability that the true coefficient falls in this range, given the data and prior. This is fundamentally different from a confidence interval—it’s a direct statement about the parameter, not about hypothetical repeated samples.

# Plot posterior distributions
mcmc_plot(model1, type = "areas", pars = "^b_")

Visualizing Posterior Distributions

bayesplot provides excellent visualizations:

# Posterior density plots
mcmc_areas(model1, pars = c("b_wt", "b_cyl"), 
           prob = 0.95) +
  labs(title = "Posterior Distributions of Coefficients")

Priors in brms

brms uses weakly informative priors by default. You can specify custom priors using the prior argument:

# Specify custom priors
model2 <- brm(
  mpg ~ wt + cyl,
  data = mtcars,
  family = gaussian(),
  prior = c(
    prior(normal(0, 10), class = Intercept),
    prior(normal(-3, 1), class = b, coef = "wt"),
    prior(normal(0, 2), class = b, coef = "cyl"),
    prior(exponential(1), class = sigma)
  ),
  chains = 4,
  seed = 123
)

Prior specifications:

• class = Intercept: Prior for the intercept
• class = b: Prior for all slope coefficients
• class = b, coef = "wt": Prior for specific coefficient
• class = sigma: Prior for residual standard deviation

Bayesian Logistic Regression

brms supports many families. Here’s logistic regression for a binary outcome:

# Create a binary outcome: high vs low efficiency
mtcars$high_eff <- ifelse(mtcars$mpg > median(mtcars$mpg), 1, 0)

# Fit Bayesian logistic regression
model3 <- brm(
  high_eff ~ wt + disp,
  data = mtcars,
  family = bernoulli(link = "logit"),
  chains = 4,
  seed = 123
)

summary(model3)

The interpretation changes slightly—coefficients are on the log-odds scale.

Model Comparison

Compare models using leave-one-out cross-validation:

# Fit models with different predictors
model_full <- brm(mpg ~ wt + cyl + disp, data = mtcars, chains = 4)
model_reduced <- brm(mpg ~ wt, data = mtcars, chains = 4)

# Compare using LOO
loo(model_full, model_reduced)

The output includes:

• elpd_loo: Expected log predictive density
• se: Standard error
• Difference: How much better the preferred model is

Multilevel/Hierarchical Models

One of brms’ strengths is handling multilevel data. Here’s a simple random-intercepts model:

# Simulate multilevel data
set.seed(123)
n_schools <- 8
n_students <- 30

school_data <- data.frame(
  school = rep(1:n_schools, each = n_students),
  student = 1:n_students
)
school_data$school_effect <- rnorm(n_schools, 0, 2)[school_data$school]
school_data$ability <- rnorm(n_students, 0, 1)
school_data$score <- 50 + school_data$school_effect + school_data$ability + rnorm(nrow(school_data), 0, 3)

# Fit multilevel model
model_school <- brm(
  score ~ ability + (1 | school),
  data = school_data,
  chains = 4,
  seed = 123
)

The (1 | school) syntax adds a random intercept for each school.

Posterior Predictions

Generate predictions from your model:

# Get posterior predictive samples
pp <- posterior_predict(model1, newdata = data.frame(wt = c(2, 3, 4), cyl = c(4, 6, 8)))

# pp is a matrix: rows = iterations, columns = new observations
dim(pp)
# [1] 4000    3

# Summarize predictions
apply(pp, 2, mean)
apply(pp, 2, quantile, probs = c(0.025, 0.975))

Checking Convergence

Always verify that your MCMC chains converged properly:

# Trace plots
mcmc_trace(model1, pars = c("b_wt", "b_cyl"))

# Rhat values (should be < 1.05)
rhat(model1)

# Effective sample sizes
neff_ratio(model1)

If Rhat > 1.1 or effective sample sizes are low, you may need more iterations or better priors.

Summary

You’ve learned the fundamentals of fitting Bayesian regression models with brms:

Concept	Description
Formula syntax	Familiar R notation with ~ for regression
Families	gaussian, bernoulli, poisson, and many more
Priors	Default weakly informative or custom-specified
Credible intervals	Direct probability statements about parameters
MCMC diagnostics	Check Rhat, ESS, and trace plots

The Bayesian approach gives you richer information than frequentist methods—you get full posterior distributions rather than single point estimates, allowing direct probability statements about your parameters.

See Also

• Introduction to Bayesian Thinking — Core concepts of Bayesian inference
• Linear Regression in R — Frequentist foundation for understanding regression concepts

Next Steps

Continue building your Bayesian skills:

• Explore more families (poisson, ordinal, zero-inflated)
• Add random effects for nested and crossed designs
• Learn about posterior predictive model checking with bayesplot
• Master prior specification for your specific domains