21  {brms}

In most cases, we do not have to write our own Stan code, because the {brms} package can do that for us. The brm() function in {brms} takes glm()-like arguments, and writes accurate and fast Stan code. If brm() works for your application, you should much prefer it to writing your own Stan models. It will (probably) be faster and is less likely to have errors.

# load packages
library(brms)

{brms} offers very flexible options for specifying your models and works with {marginaleffects}. We can use it for almost any application, though for simple lm() or glm() fits, it isn’t needed.

See ?brms::brm for the details. There are numerous tutorials available for specific models.

21.1 The turnout example

Here’s the syntax for our usual logistic regression model.

# load only the turnout data frame and hard-code rescaled variables
turnout <- ZeligData::turnout  |>
  mutate(across(age:income, arm::rescale, .names = "rs_{.col}"))

# formula
f  <- vote ~ rs_age + rs_educate + rs_income + race

# fit model with glm()
fit_glm <- glm(f, data = turnout, family = binomial)

# fit model with brm() via cmdstanr
fit_brm <- brm(f, data = turnout, family = bernoulli, 
           chains = 10, cores = 10, 
           backend = "cmdstanr")
  • [CTK: brm()’s default backend is still {rstan}; options(brms.backend = "cmdstanr") makes {cmdstanr} the default for the session] glm() and brm() obtain the same estimates.1

1 The two are not mathematically equal, but they are very close here because {brms}’s default prior on the coefficients is flat. With flat or weak priors and this much data, it is rare that a Bayesian estimate or interval differs from the ML estimate or interval, outside of cases where maximum likelihood fails.

Code 
library(modelsummary)

normalize_terms <- function(x) {
  x <- sub("^b_", "", x)  # brms to glm names
  x <- ifelse(x == "(Intercept)", "Intercept", x)  # unify intercept
  x
}

modelsummary(
  list("glm()" = fit_glm, "brm()" = fit_brm),
  coef_rename = normalize_terms,
  statistic = "conf.int",
  notes = "95% confidence intervals for glm(); 95% credible intervals for brm().",
  gof_map = NULL
)
glm() brm()
95% confidence intervals for glm(); 95% credible intervals for brm().
Intercept 1.058 1.063
[0.787, 1.335] [0.794, 1.344]
rs_age 0.993 0.999
[0.758, 1.234] [0.763, 1.244]
rs_educate 1.184 1.188
[0.918, 1.456] [0.933, 1.464]
rs_income 1.001 1.005
[0.706, 1.308] [0.708, 1.311]
racewhite 0.251 0.248
[-0.039, 0.536] [-0.045, 0.534]
Num.Obs. 2000 2000
R2 0.119
AIC 2034.0
BIC 2062.0
Log.Lik. -1011.991
F 47.751
ELPD -1017.2
ELPD s.e. 23.0
LOOIC 2034.4
LOOIC s.e. 46.0
WAIC 2034.4
RMSE 0.41 0.41

{brms} warns only when an R-hat is larger than 1.05, so we check R-hat, ESS, and divergent transitions ourselves.

library(posterior)

# R-hat and ESS for each parameter
summarise_draws(fit_brm, "rhat", "ess_bulk", "ess_tail")
# A tibble: 8 × 4
  variable      rhat ess_bulk ess_tail
  <chr>        <dbl>    <dbl>    <dbl>
1 b_Intercept   1.00    9977.    7302.
2 b_rs_age      1.00    9177.    7265.
3 b_rs_educate  1.00    9036.    6706.
4 b_rs_income   1.00   10362.    7909.
5 b_racewhite   1.00   10487.    7381.
6 Intercept     1.00    9378.    7688.
7 lprior        1.00    9378.    7688.
8 lp__          1.00    4900.    6602.
# number of divergent transitions after warmup
sum(nuts_params(fit_brm, pars = "divergent__")$Value)
[1] 0

And {brms} plays very nicely with {marginaleffects}.

library(marginaleffects)

# compute first difference
comparisons(fit_brm, variables = list(rs_age = c(-0.5, 0.5)), 
            newdata = datagrid(grid_type = "mean_or_mode"))

 Estimate 2.5 % 97.5 %
    0.167 0.128  0.207

Term: rs_age
Type: response
Comparison: 0.5 - -0.5
# compute and plot expected values (i.e., predicted probabilities)
p <- predictions(fit_brm, variables = list(rs_age = seq(-1, 1, by = 0.1)), 
            newdata = datagrid(grid_type = "mean_or_mode"))
ggplot(p, aes(x = rs_age, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_ribbon(fill = "grey") + 
  geom_line()

21.2 The flexibility of MCMC and brm()

21.2.1 Modeling “smooth” relationships

To illustrate the power of MCMC and brm(), we can imagine fitting a smooth curve for age. The details are not important for now, but this isn’t an easy task in our usual framework. But it’s straightforward with Bayesian MCMC via brm().

# smooth relationship between age and voting--notice the s()
f  <- vote ~ s(rs_age) + rs_educate + rs_income + race
smooth_fit <- brm(f, data = turnout, family = bernoulli, 
           chains = 10, cores = 10, 
           backend = "cmdstanr",
           control = list(adapt_delta = 0.99),  # avoids divergent transitions
           seed = 1234)
# compute and plot expected values (i.e., predicted probabilities)
p <- predictions(smooth_fit, variables = list(rs_age = seq(-1, 1, by = 0.1)), 
            newdata = datagrid(grid_type = "mean_or_mode"))
ggplot(p, aes(x = rs_age, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_ribbon(fill = "grey") + 
  geom_line()

  • [CTK: smooths often cause divergent transitions; with the default adapt_delta = 0.8 these two fits had 60 and 156 (of 10,000; CmdStan 2.40, seed = 1234); adapt_delta = 0.99 gives none with seed = 1234 (other seeds occasionally leave 1 of 10,000: 2 of 7 seeds checked)]

21.2.2 Awkwardly scaled variables

Also, brm() has no trouble with the variables on their original scales. The original scales pose a serious hurdle for our metrop() sampler.

# smooth relationship between age and voting--notice the s()
f  <- vote ~ s(age) + educate + income + race
smooth_fit <- brm(f, data = turnout, family = bernoulli, 
           chains = 10, cores = 10, 
           backend = "cmdstanr",
           control = list(adapt_delta = 0.99),  # avoids divergent transitions
           seed = 1234)
# compute and plot expected values (i.e., predicted probabilities)
p <- predictions(smooth_fit, variables = list(age = seq(18, 90, by = 1)), 
            newdata = datagrid(grid_type = "mean_or_mode"))
ggplot(p, aes(x = age, y = estimate, ymin = conf.low, ymax = conf.high)) + 
  geom_ribbon(fill = "grey") + 
  geom_line()