<b>brms</b>: An <i>R</i> Package for Bayesian Multilevel Models Using <i>Stan</i>

The brms package implements Bayesian multilevel models in R using the probabilistic programming language Stan. A wide range of distributions and link functions are supported, allowing users to fit - among others - linear, robust linear, binomial, Poisson, survival, ordinal, zero-inflated, hurdle, and even non-linear models all in a multilevel context. Further modeling options include autocorrelation of the response variable, user defined covariance structures, censored data, as well as meta-analytic standard errors. Prior specifications are flexible and explicitly encourage users to apply prior distributions that actually reflect their beliefs. In addition, model fit can easily be assessed and compared with the Watanabe-Akaike information criterion and leave-one-out cross-validation.

section_content: Multilevel models (MLMs) offer a great flexibility for researchers across sciences (Brown and Prescott 2015; Demidenko 2013; Gelman and Hill 2006; Pinheiro and Bates 2000). They allow the modeling of data measured on different levels at the same time -for instance data of students nested within classes and schools -thus taking complex dependency structures into account. It is not surprising that many packages for R (R Core Team 2017) have been developed to fit MLMs. Possibly the most widely known package in this area is lme4 (Bates, Mächler, Bolker, and Walker 2015), which uses maximum likelihood or restricted maximum likelihood methods for model fitting. Although alternative Bayesian methods have several advantages over frequentist approaches (e. g., the possibility of explicitly incorporating prior knowledge about parameters into the model), their practical use was limited for a long time because the posterior distributions of more complex models (such as MLMs) could not be found analytically. Markov chain Monte Carlo (MCMC) algorithms allowing to draw ran-dom samples from the posterior were not available or too time-consuming. In the last few decades, however, this has changed with the development of new algorithms and the rapid increase of general computing power. Today, several software packages implement these techniques, for instance WinBUGS (Lunn, Thomas, Best, and Spiegelhalter 2000; Spiegelhalter, Thomas, Best, and Lunn 2003), OpenBUGS (Spiegelhalter, Thomas, Best, and Lunn 2007), JAGS (Plummer 2013), MCMCglmm (Hadfield 2010) and Stan (Stan Development Team 2017a; Carpenter, Gelman, Hoffman, Lee, Goodrich, Betancour, Brubaker, Guo, Li, and Ridell 2017) to mention only a few. With the exception of the latter, all of these programs are primarily using combinations of Metropolis-Hastings updates (Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller 1953; Hastings 1970) and Gibbs sampling (Geman and Geman 1984; Gelfand and Smith 1990), sometimes also coupled with slice sampling (Damien, Wakefield, and Walker 1999; Neal 2003). One of the main problems of these algorithms is their rather slow convergence for high-dimensional models with correlated parameters (Neal 2011; Hoffman and Gelman 2014; Gelman, Carlin, Stern, and Rubin 2003). Furthermore, Gibbs sampling requires priors to be conjugate to the likelihood of parameters in order to work efficiently (Gelman et al. 2003), thus reducing the freedom of the researcher in choosing a prior that reflects his or her beliefs. In contrast, Stan implements Hamiltonian Monte Carlo (Duane, Kennedy, Pendleton, and Roweth 1987; Neal 2011 ) and its extension, the No-U-Turn Sampler (NUTS) (Hoffman and Gelman 2014). These algorithms converge much more quickly especially for high-dimensional models regardless of whether the priors are conjugate or not (Hoffman and Gelman 2014). 

Similar to software packages like WinBUGS, Stan comes with its own programming language, allowing for great modeling flexibility (cf., Stan Development Team 2017b; Carpenter et al. 2017). Many researchers may still hesitate to use Stan directly, as every model has to be written, debugged and possibly also optimized. This may be a time-consuming and error prone process even for researchers familiar with Bayesian inference. The package brms (Bürkner 2017), presented in this paper, aims at closing this gap (at least for MLMs) allowing the user to benefit from the merits of Stan only by using simple, lme4-like formula syntax. brms supports a wide range of distributions and link functions, allows for multiple grouping factors each with multiple group-level effects, autocorrelation of the response variable, user defined covariance structures, as well as flexible and explicit prior specifications. 

The purpose of the present article is to provide a general overview of the brms package (version 0. 10. 0). Package brms is available from the Comprehensive R Archive Network (CRAN) at https://CRAN. R-project. org/package=brms. We begin by explaining the underlying structure of MLMs. Next, the software is introduced in detail using recurrence times of infection in kidney patients (McGilchrist and Aisbett 1991) and ratings of inhaler instructions (Ezzet and Whitehead 1991) as examples. We end by comparing brms to other R packages implementing MLMs and describe future plans for extending the package.
section_name: Introduction

section_content: The core of every MLM is the prediction of the response y through the linear combination η of predictors transformed by the inverse link function f assuming a certain distribution D for y. We write

to stress the dependency on the ith data point. In many R packages, D is also called the "family" and we will use this term in the following. The parameter θ describes additional family specific parameters that typically do not vary across data points, such as the standard deviation σ in normal models or the shape α in Gamma or negative binomial models. The linear predictor can generally be written as

In this equation, β and u are the coefficients at population-level and group-level respectively and X, Z are the corresponding design matrices. The response y as well as X and Z make up the data, whereas β, u, and θ are the model parameters being estimated. The coefficients β and u may be more commonly known as fixed and random effects. However, we avoid these terms in the present paper following the recommendations of Gelman and Hill (2006), as they are not used unambiguously in the literature. Also, we want to make explicit that u is a model parameter in the same manner as β so that uncertainty in its estimates can be naturally evaluated. In fact, this is an important advantage of Bayesian MCMC methods as compared to maximum likelihood approaches, which do not treat u as a parameter, but assume that it is part of the error term instead (cf., Fox and Weisberg, 2011). 

Except for linear models, we do not incorporate an additional error term for every observation by default. If desired, such an error term can always be modeled using a grouping factor with as many levels as observations in the data.
section_name: Model description

section_name: Prior distributions

section_content: In brms, population-level parameters are not restricted to have normal priors. Instead, every parameter can have every one-dimensional prior implemented in Stan, for instance uniform, Cauchy or even Gamma priors. As a negative side effect of this flexibility, correlations between them cannot be modeled as parameters. If desired, point estimates of the correlations can be obtained after sampling has been done. By default, population level parameters have an improper flat prior over the reals.
section_name: Regression parameters at population-level

section_content: The group-level parameters u are assumed to come from a multivariate normal distribution with mean zero and unknown covariance matrix Σ:

As is generally the case, covariances between group-level parameters of different grouping factors are assumed to be zero. This implies that Z and u can be split up into several matrices Z k and parameter vectors u k, where k indexes grouping factors, so that the model can be simplified to

Usually, but not always, we can also assume group-level parameters associated with different levels (indexed by j) of the same grouping factor to be independent leading to

The covariance matrices V k are modeled as parameters. In most packages, an Inverse-Wishart distribution is used as a prior for V k. This is mostly because its conjugacy leads to good properties of Gibbs samplers (Gelman et al. 2003). However, there are good arguments against the Inverse-Wishart prior (Natarajan and Kass 2000; Kass and Natarajan 2006). The NUTS sampler implemented in Stan does not require priors to be conjugate. This advantage is utilized in brms: V k is parameterized in terms of a correlation matrix Ω k and a vector of standard deviations σ k through

where D(σ k ) denotes the diagonal matrix with diagonal elements σ k. Priors are then specified for the parameters on the right hand side of the equation. For Ω k, we use the LKJ-correlation prior with parameter ζ > 0 by Lewandowski, Kurowicka, and Joe (2009) 1 :

The expected value of the LKJ-prior is the identity matrix (implying correlations of zero) for any positive value of ζ, which can be interpreted like the shape parameter of a symmetric beta distribution (Stan Development Team 2017b). If ζ = 1 (the default in brms) the density is uniform over correlation matrices of the respective dimension. If ζ > 1, the identity matrix is the mode of the prior, with a sharper peak in the density for larger values of ζ. If 0 < ζ < 1 the prior is U-shaped having a trough at the identity matrix, which leads to higher probabilities for non-zero correlations. For every element of σ k, any prior can be applied that is defined on the non-negative reals only. As default in brms, we use a half Student-t prior with 3 degrees of freedom. This prior often leads to better convergence of the models than a half Cauchy prior, while still being relatively weakly informative. 

Sometimes -for instance when modeling pedigrees -different levels of the same grouping factor cannot be assumed to be independent. In this case, the covariance matrix of u k becomes

where A k is the known covariance matrix between levels and ⊗ is the Kronecker product.
section_name: Regression parameters at group-level

section_content: For some families, additional parameters need to be estimated. In the current section, we only name the most important ones. Normal and Student's distributions need the parameter σ to account for residual error variance. By default, σ has a half Cauchy prior with a scale parameter that depends on the standard deviation of the response variable to remain only weakly informative regardless of response variable's scaling. Furthermore, Student's distribution needs the parameter ν representing the degrees of freedom. By default, ν has a wide Gamma prior as proposed by Juárez and Steel (2010). Gamma, Weibull, and negative binomial distributions need the shape parameter α that also has a wide Gamma prior by default.
section_name: Family specific parameters

section_content: The brms package does not fit models itself but uses Stan on the back-end. Accordingly, all samplers implemented in Stan can be used to fit brms models. Currently, these are the static Hamiltonian Monte Carlo (HMC) sampler sometimes also referred to as hybrid Monte Carlo (Neal 2011 (Neal, 2003;; Duane et al. 1987) and its extension the no-U-turn sampler (NUTS) by Hoffman and Gelman (2014). HMC-like algorithms produce samples that are much less autocorrelated than those of other samplers such as the random-walk Metropolis algorithm (Hoffman and Gelman 2014; Creutz 1988). The main drawback of this increased efficiency is the need to calculate the gradient of the log-posterior, which can be automated using algorithmic differentiation (Griewank and Walther 2008) but is still a time-consuming process for more complex models. Thus, using HMC leads to higher quality samples but takes more time per sample than other algorithms typically applied. Another drawback of HMC is the need to pre-specify at least two parameters, which are both critical for the performance of HMC. The NUTS sampler allows setting these parameters automatically thus eliminating the need for any hand-tuning, while still being at least as efficient as a well tuned HMC (Hoffman and Gelman 2014). For more details on the sampling algorithms applied in Stan, see the Stan user's manual (Stan Development Team 2017b) as well as Hoffman and Gelman (2014). 

In addition to the estimation of model parameters, brms allows drawing samples from the posterior predictive distribution as well as from the pointwise log-likelihood. Both can be used to assess model fit. The former allows a comparison between the actual response y and the response ŷ predicted by the model. The pointwise log-likelihood can be used, among others, to calculate the Watanabe-Akaike information criterion (WAIC) proposed by Watanabe (2010) and leave-one-out cross-validation (LOO; Gelfand, Dey, and Chang 1992; Vehtari, Gelman, and Gabry 2015; see also Ionides 2008) both allowing to compare different models applied to the same data (lower WAICs and LOOs indicate better model fit). The WAIC can be viewed as an improvement of the popular deviance information criterion (DIC), which has been criticized by several authors (Vehtari et al. 2015; Plummer 2008 ; Van der Linde 2005; see also the discussion at the end of the original DIC paper by Spiegelhalter, Best, Carlin, and Van der Linde 2002) in part because of problems arising from fact that the DIC is only a point estimate. In brms, WAIC and LOO are implemented using the loo package (Vehtari, Gelman, and Gabry 2017 ) also following the recommendations of Vehtari et al. (2015).
section_name: Parameter estimation

section_content: The brms package provides functions for fitting MLMs using Stan for full Bayesian inference. 

To install the latest release version of brms from CRAN, type install. packages("brms") within R. The current development version can be downloaded and installed from GitHub via devtools::install_github("paul-buerkner/brms"). 

Additionally, a C++ compiler is required. This is because brms internally creates Stan code, which is translated to C++ and compiled afterwards. The program Rtools (R Core Team 2016) comes with a C++ compiler for Windows. 2 On OS X, one should use Xcode (Apple Inc. 2015) from the App Store. To check whether the compiler can be called within R, run system("g++ -v") when using Rtools or system("clang++ -v") when using Xcode. 

The user passes all model information to brm brm calls make stancode and make standata Model code, data, and additional arguments are passed to rstan

The model is translated to C++, compiled, and fitted in Stan

The fitted model is postprocessed within brms

Results can be investigated using various R methods defined on the fitted model object If no warning occurs and a few lines of difficult to read system code are printed out, the compiler should work correctly. For more detailed instructions on how to get the compilers running, see the prerequisites section on https://github. com/stan-dev/rstan/wiki/ RStan-Getting-Started. 

Models are fitted in brms using the following procedure, which is also summarized in Figure 1. First, the user specifies the model using the brm function in a way typical for most model fitting R functions, that is by defining formula, data, and family, as well as some other optional arguments. Second, this information is processed and the make_stancode and make_standata functions are called. The former generates the model code in Stan language and the latter prepares the data for use in Stan. These two are the mandatory parts of every Stan model and without brms, users would have to specify them themselves. Third, Stan code and data as well as additional arguments (such as the number of iterations and chains) are passed to functions of the rstan package (the R interface of Stan; Stan Development Team, 2017a). Fourth, the model is fitted by Stan after translating and compiling it in C++. Fifth, after the model has been fitted and returned by rstan, the fitted model object is post-processed in brms among others by renaming the model parameters to be understood by the user. Sixth, the results can be investigated in R using various methods such as summary, plot, or predict (for a complete list of methods type methods(class = "brmsfit")).
section_name: Software

section_content: In the following, we use an example about the recurrence time of an infection in kidney patients initially published by McGilchrist and Aisbett (1991). The data set consists of 76 entries of 7 variables:

R> library("brms") R> data("kidney", package = "brms") R> head(kidney, n = 3)
section_name: A worked example

section_content: The core of the brms package is the brm function and we will explain its argument structure using the example above. Suppose we want to predict the (possibly censored) recurrence time using a log-normal model, in which the intercept as well as the effect of age is nested within patients. Then, we may use the following code:

prior = c(set_prior("normal(0,5)", class = "b"), + set_prior("cauchy(0,2)", class = "sd"), + set_prior("lkj(2)", class = "cor")), warmup = 1000, + iter = 2000, chains = 4, control = list(adapt_delta = 0. 95))
section_name: Fitting models with brms

section_content: Without doubt, formula is the most complicated argument, as it contains information on the response variable as well as on predictors at different levels of the model. Everything before the ~sign relates to the response part of formula. In the usual and most simple case, this is just one variable name (e. g., time). However, to incorporate additional information about the response, one can add one or more terms of the form | fun(variable). fun may be one of a few functions defined internally in brms and variable corresponds to a variable in the data set supplied by the user. In this example, cens makes up the internal function that handles censored data, and censored is the variable that contains information on the censoring. Other available functions in this context are weights and disp to allow different sorts of weighting, se to specify known standard errors primarily for meta-analysis, trunc to define truncation boundaries, trials for binomial models 3, and cat to specify the number of categories for ordinal models. 

Everything on the right side of ~specifies predictors. Here, the syntax exactly matches that of lme4. For both, population-level and group-level terms, the + is used to separate different effects from each other. Group-level terms are of the form (coefs | group), where coefs contains one or more variables whose effects are assumed to vary with the levels of the grouping factor given in group. Multiple grouping factors each with multiple group-level coefficients are possible. In the present example, only one group-level term is specified in which 1 + age are the coefficients varying with the grouping factor patient. This implies that the intercept of the model as well as the effect of age is supposed to vary between patients. By default, group-level coefficients within a grouping factor are assumed to be correlated. Correlations can be set to zero by using the (coefs || group) syntax 4. Everything on the right side of formula that is not recognized as part of a group-level term is treated as a population-level effect. In this example, the population-level effects are age, sex, and disease.
section_name: formula: Information on the response and predictors

section_content: Argument family should usually be a family function, a call to a family function or a character string naming the family. If not otherwise specified, default link functions are applied. brms comes with a large variety of families. Linear and robust linear regression can be performed using the gaussian or student family combined with the identity link. For dichotomous and categorical data, families bernoulli, binomial, and categorical combined with the logit link, by default, are perfectly suited. Families poisson, negbinomial, and geometric allow for modeling count data. Families lognormal, Gamma, exponential, and weibull can be used (among others) for survival regression. Ordinal regression can be performed using the families cumulative, cratio, sratio, and acat. Finally, excess zeros in the response can be adequately modeled using families zero_inflated_poisson, zero_inflated_negbinomial, zero_inflated_binomial, zero_inflated_beta, hurdle_poisson, hurdle_negbinomial, and hurdle_gamma. In our example, we use family = lognormal() implying a log-normal "survival" model for the response variable time.
section_name: family: Distribution of the response variable

section_content: Every population-level effect has its corresponding regression parameter. These parameters are named as b_<coef>, where <coef> represents the name of the corresponding populationlevel effect. The default prior is an improper flat prior over the reals. Suppose, for instance, that we want to set a normal prior with mean 0 and standard deviation 10 on the effect of age and a Cauchy prior with location 1 and scale 2 on sexfemale 5. Then, we may write R> prior <-c(set_prior("normal(0,10)", class = "b", coef = "age"), + set_prior("cauchy(1,2)", class = "b", coef = "sexfemale"))

To put the same prior (e. g., a normal prior) on all population-level effects at once, we may write as a shortcut set_prior("normal(0,10)", class = "b"). This also leads to faster sampling, because priors can be vectorized in this case. Note that we could also omit the class argument for population-level effects, as it is the default class in set_prior. 

A special shrinkage prior to be applied on population-level effects is the horseshoe prior (Carvalho, Polson, and Scott 2009, 2010). It is symmetric around zero with fat tails and an infinitely large spike at zero. This makes it ideal for sparse models that have many regression coefficients, although only a minority of them is non-zero. The horseshoe prior can be applied on all population-level effects at once (excluding the intercept) by using set_prior("horseshoe(1)"). The 1 implies that the Student-t prior of the local shrinkage parameters has 1 degree of freedom. In brms it is possible to increase the degrees of freedom (which will often improve convergence), although the prior no longer resembles a horseshoe in this case. 6 For more details see Carvalho et al. (2009 Carvalho et al. (, 2010)). 

Each group-level effect of each grouping factor has a standard deviation parameter, which is restricted to be non-negative and, by default, has a half Student-t prior with 3 degrees of freedom and a scale parameter that is minimally 10. For non-ordinal models, brms tries to evaluate if the scale is large enough to be considered only weakly informative for the model at hand by comparing it with the standard deviation of the response after applying the link function. If this is not the case, it will increase the scale based on the aforementioned standard deviation. 7 Stan implicitly defines a half Student-t prior by using a Student-t prior on a restricted parameter (Stan Development Team 2017b). For other reasonable priors on standard deviations see Gelman (2006). In brms, standard deviation parameters are named as sd_<group>_<coef> so that sd_patient_Intercept and sd_patient_age are the parameter names in the example. If desired, it is possible to set a different prior on each parameter, but statements such as set_prior("student_t(3,0,5)", class = "sd", group = "patient") or even set_prior("student_t(3,0,5)", class = "sd") may also be used and are again faster because of vectorization. 

If there is more than one group-level effect per grouping factor, correlations between grouplevel effects are estimated. As mentioned in Section 2, the LKJ-correlation prior with parameter ζ > 0 (Lewandowski et al. 2009 ) is used for this purpose. In brms, this prior is abbreviated as "lkj(zeta)" and correlation matrix parameters are named as cor_<group>, (e. g., cor_patient), so that set_prior("lkj(2)", class = "cor", group = "patient") is a valid statement. To set the same prior on every correlation matrix in the model, set_prior("lkj(2)", class = "cor") is also allowed, but does not come with any efficiency increases. 

Other model parameters such as the residual standard deviation sigma in normal models or the shape in Gamma models have their priors defined in the same way, where each of them is treated as having its own parameter class. A complete overview on possible prior distributions is given in the Stan user's manual (Stan Development Team 2017b). Note that brms does not thoroughly check if the priors are written in correct Stan language. Instead, Stan will check their syntactical correctness when the model is parsed to C++ and return an error if they are not. This, however, does not imply that priors are always meaningful if they are accepted by Stan. Although brms tries to find common problems (e. g., setting bounded priors on unbounded parameters), there is no guarantee that the defined priors are reasonable for the model.
section_name: prior: Prior distributions of model parameters

section_content: In addition to choosing the number of iterations, warmup samples, and chains, users can control the behavior of the NUTS sampler by using the control argument. The most important reason to use control is to decrease (or eliminate at best) the number of divergent transitions that cause a bias in the obtained posterior samples. Whenever you see the warning "There were x divergent transitions after warmup. ", you should really think about increasing adapt_delta. To do this, write control = list(adapt_delta = <x>), where <x> should usually be a value between 0. 8 (current default) and 1. Increasing adapt_delta will slow down the sampler but will decrease the number of divergent transitions threatening the validity of your posterior samples. 

Another problem arises when the depth of the tree being evaluated in each iteration is exceeded. This is less common than having divergent transitions, but may also bias the posterior samples. When it happens, Stan will throw out a warning suggesting to increase max_treedepth, which can be accomplished by writing control = list(max_treedepth = <x>) with a positive integer <x> that should usually be larger than the current default of 10.
section_name: control: Adjusting the sampling behavior of Stan

section_content: The example model fit1 is fitted using 4 chains, each with 2000 iterations of which the first 1000 are warmup to calibrate the sampler, leading to a total of 4000 posterior samples. 8 For researchers familiar with Gibbs or Metropolis-Hastings sampling, this number may seem far too small to achieve good convergence and reasonable results, especially for multilevel models. However, as brms utilizes the NUTS sampler (Hoffman and Gelman 2014) implemented in Stan, even complex models can often be fitted with not more than a few thousand samples. Of course, every iteration is more computationally intensive and time-consuming than the iterations of other algorithms, but the quality of the samples (i. e., the effective sample size per iteration) is usually higher. 

After the posterior samples have been computed, the brm function returns an R object, containing (among others) the fully commented model code in Stan language, the data to fit the model, and the posterior samples themselves. The model code and data for the present example can be extracted through stancode(fit1) and standata(fit1 Samples were drawn using sampling(NUTS). For each parameter, Eff. Sample is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1). 

On the top of the output, some general information on the model is given, such as family, formula, number of iterations and chains, as well as the WAIC. Next, group-level effects are displayed separately for each grouping factor in terms of standard deviations and correlations between group-level effects. On the bottom of the output, population-level effects are displayed. If incorporated, autocorrelation and family specific parameters (e. g., the residual standard deviation sigma) are also given. 

In general, every parameter is summarized using the mean (Estimate) and the standard deviation (Est. Error) of the posterior distribution as well as two-sided 95% credible intervals (l-95% CI and u-95% CI) based on quantiles. The Eff. Sample value is an estimation of the effective sample size; that is the number of independent samples from the posterior distribution that would be expected to yield the same standard error of the posterior mean as is obtained from the dependent samples returned by the MCMC algorithm. The Rhat value provides information on the convergence of the algorithm (cf., Gelman and Rubin, 1992). If Rhat is considerably greater than 1 (i. e., > 1. 1), the chains have not yet converged and it is necessary to run more iterations and/or set stronger priors. 

To visually investigate the chains as well as the posterior distribution, the plot method can be used (see Figure 2 ). An even more detailed investigation can be achieved by applying the shinystan package (Gabry 2017) through method launch_shiny. With respect to the above summary, sexfemale seems to be the only population-level effect with considerable influence on the response. Because the mean of sexfemale is positive, the model predicts longer periods without an infection for females than for males. Effects of population-level predictors can also be visualized with the marginal_effects method (see Figure 3 ). 

Looking at the group-level effects, the standard deviation parameter of age is suspiciously small. To test whether it is smaller than the standard deviation parameter of Intercept, we apply the hypothesis method:

R> hypothesis(fit1, "Intercept -age > 0", class = "sd", group = "patient") The one-sided 95% credible interval does not contain zero, thus indicating that the standard deviations differ from each other in the expected direction. In accordance with this finding, the Evid. Ratio shows that the hypothesis being tested (i. e., Intercept -age > 0) is about 50 times more likely than the alternative hypothesis Intercept -age < 0. It is important to note that this kind of comparison is not easily possible when applying frequentist methods, because in this case only point estimates are available for group-level standard deviations and correlations. 

When looking at the correlation between both group-level effects, its distribution displayed in Figure 2 and the 95% credible interval in the summary output appear to be rather wide. This indicates that there is not enough evidence in the data to reasonably estimate the correlation. Together, the small standard deviation of age and the uncertainty in the correlation raise the question if age should be modeled as a group-specific term at all. To answer this question, we fit another model without this term:

A good way to compare both models is leave-one-out cross-validation (LOO) 10, which can be called in brms using approximately normal if the number of observations is large so that the standard errors can be very helpful in evaluating differences in the information criteria. However, for small sample sizes, standard errors should be interpreted with care (Vehtari et al. 2015 ). For the present example, it is immediately evident that both models have very similar fit, indicating that there is little benefit in adding group-specific coefficients for age.
section_name: Analyzing the results

section_content: In the following, we want to briefly discuss a second example to demonstrate the capabilities of brms in handling ordinal data. Ezzet and Whitehead (1991) analyze data from a twotreatment, two-period crossover trial to compare 2 inhalation devices for delivering the drug salbutamol in 286 asthma patients. Patients were asked to rate the clarity of leaflet instructions accompanying each device, using a four-point ordinal scale. Ratings are predicted by treat to indicate which of the two inhaler devices was used, period to indicate the time of administration, and carry to model possible carry over effects. R> data("inhaler", package = "brms") R> head(inhaler, n = 1) subject rating treat period carry 1 1 1 0. 5 0. 5 0

Typically, the ordinal response is assumed to originate from the categorization of a latent continuous variable. That is there are K latent thresholds (model intercepts), which partition the continuous scale into the K + 1 observable, ordered categories. Following this approach leads to the cumulative or graded-response model (Samejima 1969) for ordinal data implemented in many R packages. In brms, it is available via family cumulative. Fitting the cumulative model to the inhaler data, also incorporating an intercept varying by subjects, may be accomplished using:

R> fit3 <-brm(formula = rating ~treat + period + carry + (1 | subject), + data = inhaler, family = cumulative)

While the support for ordinal data in many R packages ends here 11, brms allows changes to this basic model in at least three ways. First of all, three additional ordinal families are implemented. Families sratio (stopping ratio) and cratio (continuation ratio) are socalled sequential models (Tutz 1990 ). Both are equivalent to each other for symmetric link functions such as logit but will differ for asymmetric ones such as cloglog. The fourth ordinal family is acat (adjacent category) also known as partial credits model (Masters 1982; Andrich 1978b ). Second, restrictions to the thresholds can be applied. By default, thresholds are ordered for family cumulative or are completely free to vary for the other families. This is indicated by argument threshold = "flexible" (default) in the ordinal family functions. Using threshold = "equidistant" forces the distance between two adjacent thresholds to be the same, that is

for thresholds τ k and distance δ (see also Andrich 1978a; Andrich 1978b; Andersen 1977). Third, the assumption that predictors have constant effects across categories may be relaxed for non-cumulative ordinal models (Van der Ark 2001; Tutz 2000) leading to category-specific effects. For instance, variable treat may only have an impact on the decision between category 3 and 4, but not on the lower categories. Without using category-specific effects, such a pattern would remain invisible. 

To illustrate all three modeling options at once, we fit a (hardly theoretically justified) stopping ratio model with equidistant thresholds and category-specific effects for variable treat on which we apply an informative prior. 

R> fit4 <-brm(formula = + rating ~period + carry + cse(treat) + (1 | subject), + data = inhaler, family = sratio(threshold = "equidistant"), + prior = set_prior("normal(-1,2)", coef = "treat"))

Note that priors are defined on category-specific effects in the same way as for other populationlevel effects. A model summary can be obtained in the same way as before:

R> summary(fit4, waic = TRUE)

Family: sratio(logit) Formula: rating ~period + carry + cse Samples were drawn using sampling(NUTS). For each parameter, Eff. Sample is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence, Rhat = 1). 

Trace and density plots of the model parameters as produced by plot(fit4) can be found in Figure 4. We see that three intercepts (thresholds) and three effects of treat have been estimated, because a four-point scale was used for the ratings. The treatment effect seems to be strongest between category 3 and 4. At the same time, however, the credible interval is also much larger. In fact, the intervals of all three effects of treat are highly overlapping, which indicates that there is not enough evidence in the data to support category-specific effects. 

On the bottom of the output, parameter delta specifies the distance between two adjacent thresholds and indeed the intercepts differ from each other by the magnitude of delta.
section_name: Modeling ordinal data

section_content: Over the years, many R packages have been developed that implement MLMs, each being more or less general in their supported models. Comparing all of them to brms would be too extensive and barely helpful for the purpose of the present paper. Accordingly, we concentrate on a comparison with four packages. These are lme4 (Bates et al. 2015) and MCMCglmm (Hadfield 2010), which are possibly the most general and widely applied R packages for MLMs, as well as rstanarm (Gabry and Goodrich 2017) and rethinking (McElreath 2016), which are both based on Stan. As opposed to the other packages, rethinking was primarily written for teaching purposes and requires the user to specify the full model explicitly using its own simplified BUGS-like syntax thus helping users to better understand the models that are fitted to their data. 

Regarding model families, all five packages support the most common types such as linear and binomial models as well as Poisson models for count data. Currently, brms and MCM-Cglmm provide more flexibility when modeling categorical and ordinal data. In addition, brms supports robust linear regression using Student's distribution. MCMCglmm allows fitting multinomial models that are currently not available in the other packages. 

Generalizing classical MLMs, brms and MCMCglmm allow fitting zero-inflated and hurdle models dealing with excess zeros in the response. Furthermore, brms supports non-linear models similar to the nlme package (Pinheiro, Bates, DebRoy, Sarkar, and R Core Team 2017) providing great flexibility but also requiring more care to produce reasonable results. 

Another flexible model class are generalized additive mixed models (Hastie and Tibshirani 1990; Wood 2004; Zuur 2014), which can be fitted with brms and rstanarm. 

In all five packages, there are quite a few additional modeling options. Variable link functions can be specified in all packages except for MCMCglmm, in which only one link is available per family. MCMCglmm generally supports multivariate responses using data in wide format, whereas brms currently only offers this option for families gaussian and student. It should be noted that it is always possible to transform data from wide to long format for compatibility with the other packages. Autocorrelation of the response can only be fitted in brms, which supports auto-regressive as well as moving-average effects. For ordinal models in brms, effects of predictors may vary across different levels of the response as explained in the inhaler example. A feature currently exclusive to rethinking is the possibility to impute missing values in the predictor variables. 

Information criteria are available in all three packages. The advantage of WAIC and LOO implemented in brms, rstanarm, and rethinking is that their standard errors can be easily estimated to get a better sense of the uncertainty in the criteria. Comparing the prior options of the Bayesian packages, brms and rethinking offer a little more flexibility than MCMCglmm and rstanarm, as virtually any prior distribution can be applied on population-level effects as well as on the standard deviations of group-level effects. In addition, we believe that the way priors are specified in brms and rethinking is more intuitive as it is directly evident what prior is actually applied. A more detailed comparison of the packages can be found in Tables 1 and 2. To facilitate the understanding of the model formulation in brms, Table 3 shows lme4 function calls to fit sample models along with the equivalent brms syntax. 

So far the focus was only on capabilities. Another important topic is speed, especially for more complex models. Of course, lme4 is usually much faster than the other packages as it uses maximum likelihood methods instead of MCMC algorithms, which are slower by design. 

To compare the efficiency of the four Bayesian packages, we fitted multilevel models on real data sets using the minimum effective sample size divided by sampling time as a measure of sampling efficiency. One should always aim at running multiple chains as one cannot be sure that a single chain really explores the whole posterior distribution. However, as MCMCglmm does not come with a built-in option to run multiple chains, we used only a single chain to fit the models after making sure that it leads to the same results as multiple chains. The R code allowing to replicate the results is available as supplementary material. 

The first thing that becomes obvious when fitting the models is that brms and rethinking need to compile the C++ model before actually fitting it, because the Stan code being parsed to C++ is generated on the fly based on the user's input. Compilation takes about a half to one minute depending on the model complexity and computing power of the machine. This is not required by rstanarm and MCMCglmm, although the former is also based on Stan, as compilation takes place only once at installation time. While the latter approach saves the compilation time, the former is more flexible when it comes to model specification. For small and simple models, compilation time dominates the overall computation time, but for larger and more complex models, sampling will take several minutes or hours so that one minute more or less will not really matter, anymore. Accordingly, the following comparisons do not include the compilation time. 

In models containing only group-specific intercepts, MCMCglmm is usually more efficient than the Stan packages. However, when also estimating group-specific slopes, MCMCglmm falls behind the other packages and quite often refuses to sample at all unless one carefully specifies informative priors. Note that these results are obtained by running only a single chain. For all three Stan packages, sampling efficiency can easily be increased by running multiple chains in parallel. Comparing the Stan packages to each other, brms is usually most efficient for models with group-specific terms, whereas rstanarm tends to be roughly 50% to 75% as efficient at least for the analyzed data sets. The efficiency of rethinking is more variable depending on the model formulation and data, sometimes being slightly ahead of the other two packages, but usually being considerably less efficient. Generally, rethinking loses efficiency for models with many population-level effects presumably because one cannot use design matrices and vectorized prior specifications for population-level parameters. Note that it was not possible to specify the exact same priors across packages due to varying Table 3 : Comparison of the model syntax of lme4 and brms using data sets included in lme4. 

Notes: (1) Default links are used so that the link argument may be omitted. 

(2) Fitting this model takes some time. A proper prior on the population-level effects (e. g., prior = set_prior("normal(0,5)")) may help in increasing sampling speed. 

parameterizations. Of course, efficiency depends heavily on the model, chosen priors, and data at hand so that the present results should not be over-interpreted.
section_name: Comparison between packages

section_content: The present paper is meant to provide a general overview on the R package brms implementing MLMs using the probabilistic programming language Stan for full Bayesian inference. 

Although only a small selection of the modeling options available in brms are discussed in detail, I hope that this article can serve as a good starting point to further explore the capabilities of the package. 

For the future, I have several plans on how to improve the functionality of brms. I want to include multivariate models that can handle multiple response variables coming from different distributions as well as new correlation structures for instance for spatial data. Similarly, distributional regression models as well as mixture response distributions appear to be valuable extensions of the package. I am always grateful for any suggestions and ideas regarding new features.
section_name: Conclusion

section_name: Acknowledgments

section_content: First of all, I would like to thank the Stan Development Team for creating the probabilistic programming language Stan, which is an incredibly powerful and flexible tool for performing full Bayesian inference. Without it, brms could not fit a single model. Two anonymous reviewers provided very detailed and thoughtful suggestions to substantially improve both the package and the paper. Furthermore, Prof. Philipp Doebler and Prof. Heinz Holling have given valuable feedback on earlier versions of the paper. Lastly, I want to thank the many users who reported bugs or had ideas for new features, thus helping to continuously improve brms.
section_name: Acknowledgments