brms is a free R package that lets you fit Bayesian statistical models using simple, readable code. Developed by Paul-Christian Bürkner, it acts as a user-friendly front end to Stan, a powerful but complex probabilistic programming language. Instead of writing Stan code from scratch, you describe your model with a formula that looks a lot like standard R syntax, and brms translates it into Stan behind the scenes, runs the analysis, and hands you back the results.
The name stands for “Bayesian Regression Models using Stan,” and the package has become one of the most widely used tools for Bayesian analysis in research. It’s flexible enough to handle everything from basic linear regression to complex hierarchical models, non-linear relationships, and dozens of outcome types.
What Problem brms Solves
Bayesian statistics offers real advantages over traditional methods: you get full probability distributions for your estimates rather than single point values, you can incorporate prior knowledge, and you can make direct probability statements about your results. But the software side has historically been a barrier. Languages like Stan give you enormous modeling flexibility, but every model has to be written, debugged, and optimized by hand. That’s time-consuming and error-prone, even for experienced statisticians.
brms closes that gap. If you’ve used R’s lme4 package for mixed-effects models, the formula syntax will feel immediately familiar. You write something like brm(y ~ x1 + x2 + (1 | group), data = mydata), and brms generates the Stan code, compiles it, runs four independent sampling chains, and gives you a full posterior distribution for every parameter. What might take hours of Stan programming takes a single function call.
Supported Model Types
The range of models brms can fit is genuinely broad. It supports over 40 response distribution families, covering most situations a researcher would encounter:
- Continuous outcomes: Gaussian (standard linear regression), Student-t (robust regression that handles outliers), log-normal, skew-normal, and others
- Count data: Poisson, negative binomial, and geometric distributions for things like event counts or disease incidence
- Binary and categorical outcomes: Bernoulli and binomial families for yes/no data (logistic regression), plus categorical and ordinal models for multi-category outcomes
- Zero-inflated and hurdle models: For data with excess zeros, like healthcare utilization counts where many patients have zero visits
- Survival and time-to-event data: Weibull, exponential, and other distributions used in clinical and epidemiological research
- Specialized distributions: Beta for proportions, Dirichlet for compositional data, von Mises for circular data, and ex-Gaussian and Wiener distributions used in cognitive science reaction-time studies
Each family supports multiple link functions. For binary outcomes, for example, you can choose logit, probit, or complementary log-log links depending on what makes sense for your data.
Multilevel and Hierarchical Models
Where brms really shines is in multilevel modeling, sometimes called hierarchical or mixed-effects modeling. These models handle data with natural grouping structures: patients within hospitals, students within schools, repeated measurements within individuals.
The syntax extends lme4’s conventions in useful ways. You can specify multiple grouping factors, each with their own random effects. Nested structures work intuitively: writing (1 | school/classroom) automatically expands to separate terms for schools and classrooms within schools. If modeling correlations between random effects causes convergence problems, switching | to || drops those correlations to simplify the model.
brms also supports multi-membership structures, where an observation belongs to more than one group simultaneously. A student who transferred between two schools during a study, for instance, can be linked to both. This is handled through a special mm() function inside the grouping term.
Non-Linear Models
Beyond the standard linear predictor framework, brms can estimate genuinely non-linear relationships. You write the non-linear formula directly, flag it with nl = TRUE, and specify which symbols are parameters to be estimated versus predictors in the data.
For example, a growth curve model like cumulative ~ ultimate * (1 - exp(-(time/theta)^omega)) can be fit with each parameter getting its own sub-model. The parameter “ultimate” could vary by group as a random effect, while “theta” and “omega” stay fixed across groups. Each non-linear parameter is actually a placeholder for its own linear predictor, so you can build multilevel structure into any part of the non-linear equation. This makes brms unusually flexible for pharmacokinetic models, learning curves, dose-response relationships, and item-response models in psychometrics.
Setting Priors
Bayesian analysis requires prior distributions, which represent what you know (or assume) about parameters before seeing the data. brms handles this through a set_prior() function where you specify distributions using Stan’s syntax.
You can set priors on individual coefficients, on entire classes of parameters at once, or on specific model components like the standard deviation of random effects. A prior like set_prior("normal(0, 5)", class = "b", coef = "age") puts a normal distribution centered at zero with a standard deviation of 5 on the age coefficient. For correlation matrices in random effects, brms uses the LKJ prior, which with its default setting treats all correlation structures as equally plausible.
If you don’t specify priors, brms uses flat (non-informative) priors for population-level effects by default. This makes it easy to get started, but you’ll get better results, and faster sampling, by choosing weakly informative priors that reflect reasonable parameter ranges for your problem.
Checking and Comparing Models
Once a model is fit, brms provides tools to evaluate whether it captures the data well and to compare competing models. Posterior predictive checks let you simulate new data from the fitted model and visually compare it against the observed data. These plots, generated through the bayesplot package, make it straightforward to spot problems like a model that predicts too few zeros or misses the spread of the data.
For formal model comparison, brms integrates with the loo package to compute approximate leave-one-out cross-validation. Rather than refitting the model hundreds of times with one observation removed each time, it uses a technique called Pareto-smoothed importance sampling to efficiently estimate how well the model would predict each data point if that point had been left out of the fitting process. The result is a measure called ELPD (expected log pointwise predictive density), and models with higher ELPD values predict new data better. The method also provides diagnostics that warn you when the approximation is unreliable for specific observations.
How brms Compares to Alternatives
The most common alternatives in R are rstanarm (another Stan-based package) and lme4 (frequentist mixed models). brms occupies a middle ground that leans toward maximum flexibility. rstanarm uses pre-compiled Stan models, which means faster startup but a narrower range of model types. brms generates and compiles new Stan code for each model, adding some overhead but allowing it to handle non-linear models, distributional models, and a much wider set of response families.
Compared to older Bayesian tools that rely on Gibbs sampling or Metropolis-Hastings algorithms, brms benefits from Stan’s use of the No-U-Turn Sampler (NUTS), a form of Hamiltonian Monte Carlo. This sampling method converges more reliably for high-dimensional models with correlated parameters, which is exactly the situation you encounter in complex multilevel models.
Use in Regulatory and Applied Research
brms has moved well beyond academic statistics departments. The U.S. Food and Drug Administration has published example statistical analysis plans that use brms as the primary analysis tool. One such plan for a parallel-group randomized clinical trial of a hypertension treatment specified a Bayesian t-test implemented through brms, using informative priors on treatment differences and allowing for unequal variances and heavier-tailed distributions than the normal. The code called brm() directly, with explicit prior specifications for each parameter class, and ran four independent chains of 4,000 iterations each.
This kind of adoption signals that brms has reached a level of reliability and transparency that satisfies rigorous review standards. The fact that it generates readable Stan code means reviewers can inspect exactly what the model is doing, which matters in high-stakes applications.
The Other BRMS
If you arrived here searching for a medical term, BRMS also stands for the Bech-Rafaelsen Melancholia Scale, a clinician-rated tool for assessing the severity of depression. It’s a psychiatric rating scale with high discriminant validity for distinguishing mild, moderate, and severe depressive episodes as defined by international diagnostic criteria. The scale correlates well with other clinician-rated depression measures (around 0.70) but less strongly with self-report questionnaires (around 0.32 to 0.53), which is a common pattern reflecting the difference between how clinicians and patients perceive symptom severity.