brms

Overview

The brms package provides an interface to fit Bayesian generalized
(non-)linear multivariate multilevel models using Stan, which is a C++
package for performing full Bayesian inference (see
http://mc-stan.org/). The formula syntax is very similar to that of
the package lme4 to provide a familiar and simple interface for
performing regression analyses. A wide range of response distributions
are supported, allowing users to fit – among others – linear, robust
linear, count data, survival, response times, ordinal, zero-inflated,
and even self-defined mixture models all in a multilevel context.
Further modeling options include non-linear and smooth terms,
auto-correlation structures, censored data, missing value imputation,
and quite a few more. In addition, all parameters of the response
distribution can be predicted in order to perform distributional
regression. Multivariate models (i.e., models with multiple response
variables) can be fit, as well. Prior specifications are flexible and
explicitly encourage users to apply prior distributions that actually
reflect their beliefs. Model fit can easily be assessed and compared
with posterior predictive checks, cross-validation, and Bayes factors.

Resources

• Introduction to
  brms (Journal of
  Statistical Software)
• Advanced multilevel modeling with
  brms
  (The R Journal)
• Website (Website of brms
  with documentation and vignettes)
• Blog posts
  (List of blog posts about brms)
• Ask a question (Stan Forums on
  Discourse)
• Open an issue
  (GitHub issues for bug reports and feature requests)

How to use brms

library(brms)

As a simple example, we use poisson regression to model the seizure
counts in epileptic patients to investigate whether the treatment
(represented by variable Trt) can reduce the seizure counts and
whether the effect of the treatment varies with the (standardized)
baseline number of seizures a person had before treatment (variable
zBase). As we have multiple observations per person, a group-level
intercept is incorporated to account for the resulting dependency in the
data.

fit1 <- brm(count ~ zAge + zBase * Trt + (1