brms in R: Bayesian Regression Without Writing a Single Line of Stan

In the previous post you fit a Stan model the long way: write the model in the Stan language, compile it, sample, inspect. That is the right path when you need a custom likelihood. For most regressions you do not. The brms package takes the lme4-style formula syntax you already know (y ~ x + (1 | group)), generates the Stan code for you, compiles it, and gives you a fitted model with all the same diagnostics. This post walks through fitting a Bayesian regression with brms end to end, setting priors, checking the fit, and choosing between brms, rstanarm, and writing Stan yourself.

What can brms do that lm() and glmer() can't?

Frequentist tools like lm() and glmer() fit one set of point estimates plus standard errors. They do not give you a posterior over the unknowns, and they do not let you bring prior knowledge into the fit. brms does both. You write the same formula you would hand to glmer(), and brms returns a posterior over every coefficient, every variance component, and every prediction.

The win is most obvious in three places. The first is uncertainty quantification. Frequentist confidence intervals are statements about repeated sampling under the null; Bayesian credible intervals are direct probability statements over the parameter ("90% probability that the slope is between 0.4 and 0.6"). Stakeholders read the Bayesian version more easily.

The second is partial pooling in hierarchical models. glmer() does a one-shot REML fit; brms gives you the joint posterior over all group-level intercepts and slopes plus the population-level parameters, with proper uncertainty propagation between levels.

The third is flexibility. brms supports custom families (student, negbinomial, bernoulli, categorical, ordinal, zero_inflated_*, hurdle_*, and dozens more), distributional regression (where the variance can also be modelled as a function of predictors), and non-linear models, all in formula syntax.

Below is a complete brms fit on the built-in mtcars dataset. We model fuel efficiency (mpg) as a function of weight and horsepower. The brms call is one line. The output is a posterior over both coefficients with built-in convergence diagnostics.

Walk through what just happened. Line 1 loaded brms; line 2 told brms to use cmdstanr as the Stan backend (faster than rstan and tracks the latest Stan releases). The brm() call took the formula mpg ~ wt + hp, the dataset, and standard MCMC settings. Behind the scenes brms generated Stan code, compiled it (30 to 60 seconds the first time), and ran four chains.

The summary printed the posterior estimates. Estimate is the posterior mean for each coefficient. Est.Error is the posterior standard deviation. l-95% CI and u-95% CI are the bounds of the 95% credible interval.

Now read the numbers. The intercept posterior mean is 37.2 mpg, meaning a hypothetical car with zero weight and zero horsepower would average about 37.2 mpg. The wt posterior mean is -3.85 mpg per 1000 lbs of weight, with 95% credible interval [-5.16, -2.44]. The hp posterior mean is -0.03 mpg per horsepower with 95% credible interval [-0.05, -0.01].

Both intervals exclude zero, so we are confident weight and horsepower both reduce fuel efficiency. The Rhat of 1.00 and Bulk_ESS near 4000 confirm the chains converged.

That is the whole pipeline. One formula, one dataset, four chains, full Bayesian inference. Compare to writing the equivalent Stan model from scratch (data block, parameters block, model block with priors and likelihood, plus the cmdstanr scaffolding from the previous post): that took 30 lines for one model, this took five lines and is more flexible.

Figure 1: brms takes an R formula, generates Stan code, hands it to cmdstanr or rstan to compile and sample, and returns a brmsfit object that supports the standard summary, plot, and predict methods.

How do I install brms and connect it to cmdstanr?

brms is on CRAN. The Stan backend is not. You need both, and brms needs to know which backend to call.

Walk through what each step does. The first installs brms itself from CRAN. brms is a thin layer of R code that does the formula-to-Stan translation; it does not contain a sampler.

The second installs cmdstanr from the Stan r-universe (not CRAN, because the package updates faster than CRAN allows). The third compiles the actual Stan engine if you do not already have it; you only run this once per machine.

The fourth tells brms to use cmdstanr as its Stan backend for the rest of this R session. Without that line, brms defaults to rstan, which is slower to compile and lags Stan releases. To make cmdstanr the default for every session, add options(brms.backend = "cmdstanr") to your ~/.Rprofile.

What does a brms model formula actually compile to?

brms reads your formula and emits Stan code. You rarely need to look at it, but seeing the translation once demystifies what is happening. The function make_stancode() returns the generated Stan program without fitting anything.

Walk through the generated code. The data block declares N (number of observations), Y (response vector), K (number of population-level effects), and X (the design matrix). brms built X for you by parsing ~ wt + hp and pulling the relevant columns from mtcars.

The parameters block declares b (the coefficient vector for wt and hp), the Intercept, and sigma (residual standard deviation).

The model block has three lines. The first sets a Student-t(3, 19.2, 7.4) prior on the intercept, with the location and scale derived from the data's mean and standard deviation. The second sets a half-Student-t(3, 0, 7.4) prior on sigma (truncated at zero because sigma must be positive). The third is the likelihood: each Y[i] is Normal-distributed around Intercept + X[i] %*% b with standard deviation sigma. brms uses the optimised normal_id_glm_lpdf form for speed.

The data-driven defaults come from default_prior(), which inspects your data and picks sensible defaults. Look at them with prior_summary(fit_mtcars):

Walk through the columns. class is the parameter type. coef is the coefficient name. prior is the assigned prior. The (flat) lines for b and b (hp, wt) say the regression coefficients use an implicit improper uniform prior, which is the brms default for population-level effects.

The student_t(3, 19.2, 7.4) is the Student-t prior on the intercept, with 19.2 (close to mean(mpg) = 20.1) as the centre and 7.4 (close to sd(mpg) = 6.0 plus a 2.5 floor) as the scale. The student_t(3, 0, 7.4) on sigma is a half-Student-t (the lb column shows 0), the standard weakly informative choice for residual standard deviation.

How do I set priors in brms?

Default priors are weakly informative. When you have real prior knowledge, override them with prior() and pass to brm(). The prior() function takes a Stan-language distribution string and a parameter class.

Walk through what we set. The first prior, normal(0, 5) on class b, applies to every population-level coefficient (here wt and hp). The 0 mean says we have no prior reason to believe the effects are positive or negative; the 5 standard deviation says we expect the effects to be within roughly [-10, 10] mpg per unit predictor.

The second prior, normal(35, 5) on Intercept, says we expect a typical car's mpg with predictors at zero to be near 35, with reasonable spread. (Note: brms internally centres the predictors, so the "intercept" prior applies to a centred intercept by default. Read ?brmsformula for the gory detail.)

The third prior, student_t(3, 0, 5) on sigma, is a tighter half-Student-t than the default, expressing more confidence that residuals will be small.

prior_summary() confirms what was actually used. The source column shows user for everything we set explicitly.

How do I check the fit?

Three diagnostics matter for every brms fit. Convergence diagnostics, posterior predictive checks, and out-of-sample predictive accuracy.

The first you already saw: Rhat, Bulk_ESS, Tail_ESS in the summary output. brms also reports divergent transitions if any occur. For our fit_mtcars model, all those values were healthy.

The second is a posterior predictive check. It generates predicted values from the fitted model and compares them to the observed data. If the model captures the data's structure, simulated draws should look like the real data.

Walk through what pp_check() did. It drew 50 sets of fake mpg values from the fitted model (using posterior_predict() internally) and overlaid their density curves on the observed data's density. If the model is well-specified, the two should be visually indistinguishable.

For our model, the curves overlap closely in the bulk of the distribution but slightly miss the long right tail. That hint suggests the Normal likelihood may be too restrictive for the highest-mpg cars; we could try family = student("identity") for heavier tails. The next post in this section drills deeper into posterior predictive checks.

The third diagnostic is cross-validated predictive accuracy via loo(), which estimates expected log predictive density on held-out data without refitting.

Walk through what loo() reports. elpd_loo is the expected log predictive density (higher is better); p_loo is the effective number of parameters; looic is -2 * elpd_loo for comparing models on the AIC scale (lower is better). The "All Pareto k estimates are good" line confirms the LOO approximation is reliable.

To compare two models, fit both and call loo_compare(). The model with the higher elpd_loo is preferred; if the difference is at least 4 standard errors, the preference is definitive.

When should I reach for brms vs rstanarm vs cmdstanr?

Three R interfaces cover almost all Bayesian modelling. The choice is mostly about flexibility versus speed of iteration.

Figure 2: brms is the right default when you want formula syntax. rstanarm is faster for standard GLM/multilevel designs. cmdstanr gives you full Stan-language control for custom models.

Walk through the decision flow. If your model is a standard regression (linear, GLM, multilevel) and you want it to fit fast, use rstanarm. Its models are pre-compiled into the package, so there is no compile step. The trade is that rstanarm only supports the model classes its authors built in.

If your model is almost standard but needs something rstanarm does not have (a custom family, a specific link, a non-linear term, distributional regression on the variance), use brms. brms generates Stan code for a much wider range of models, at the cost of a 30 to 60 second compile per model.

If your model needs a custom likelihood that brms does not support (a mixture with non-standard components, a survival model with complex censoring, a state-space model), drop down to cmdstanr and write the Stan code yourself.

For everything from this post through the rest of the Bayesian section, our default tool is brms.

Practice Exercises

Exercise 1: Logistic regression with brms

Fit a Bayesian logistic regression on the mtcars dataset, predicting whether a car has 8 cylinders (vs == 0) from mpg and wt. Use family = bernoulli("logit") and the brms defaults.

Exercise 2: Hierarchical model with random intercepts

Fit a Bayesian linear model on mtcars with mpg as the outcome, wt as the fixed effect, and a random intercept for cyl (treating each cyl group as having its own baseline mpg).

Exercise 3: Compare three models with LOO

Fit mpg ~ wt, mpg ~ wt + hp, and mpg ~ wt * hp (with interaction). Compare them with loo_compare() and report which one wins.

Summary

brms is the standard R interface for fitting Bayesian regressions when you want formula syntax instead of writing Stan code. It generates Stan, calls cmdstanr or rstan to sample, and returns a brmsfit object that supports the full Bayesian workflow.

Use brms by default for Bayesian regression. Drop down to cmdstanr for custom likelihoods, and lift up to rstanarm when your model is a standard GLM and you want zero compile time.

References

1. Bürkner, P. "brms: An R Package for Bayesian Multilevel Models Using Stan." Journal of Statistical Software 80 (2017). The brms paper.
2. Bürkner, P. brms documentation. paulbuerkner.com/brms. The full reference manual and vignettes.
3. Bürkner, P. Applied Bayesian Regression Modelling with brms. paulbuerkner.com/software/brms-book. Book-length treatment with worked examples.
4. McElreath, R. Statistical Rethinking, 2nd ed. CRC Press, 2020. Uses brms in many examples; the best Bayesian textbook for working data analysts.
5. Vehtari, A., Gelman, A., Gabry, J. "Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC." Statistics and Computing 27 (2017). The methodology behind loo().

Continue Learning

• Stan in R, the previous post in this curriculum. Read it to understand what brms is generating under the hood when you call brm().
• Hamiltonian Monte Carlo in R, the algorithm Stan and brms run beneath the formula syntax.
• Bayesian Statistics in R, the section opener for the prior-likelihood-posterior intuition every brms model relies on.