There are several reasons why everyone isn’t using Bayesian methods for regression modeling. One reason is that Bayesian modeling requires more thought: you need pesky things like priors, and you can’t assume that if a procedure runs without throwing an error that the answers are valid. A second reason is that MCMC sampling — the bedrock of practical Bayesian modeling — can be slow compared to closed-form or MLE procedures. A third reason is that existing Bayesian solutions have either been highly-specialized (and thus inflexible), or have required knowing how to use a generalized tool like BUGS, JAGS, or Stan. This third reason has recently been shattered in the R world by not one but two packages:
rstanarm. Interestingly, both of these packages are elegant front ends to Stan, via
This article describes
rstanarm, how they help you, and how they differ.
You can install both packages from CRAN, making sure to install dependencies so you get
shinystan as well. [EDIT: Note that you will also need a C++ compiler to make this work, as Charles and Paul Bürkner (@paulbuerkner) have found. On Windows that means you’ll need to install
Rtools, and on the Mac you may have to install Xcode (which is also free).
brms‘s help refers to the RStan Getting Started, which is very helpful.] If you like having the latest development versions — which may have a few bug fixes that the CRAN versions don’t yet have — you can use
devtools to install them following instructions at the
brms github site or the
Let’s start with a quick multinomial logistic regression with the famous Iris dataset, using
brms. You may want to skip the actual
brm call, below, because it’s so slow (we’ll fix that in the next step):
library (brms) rstan_options (auto_write=TRUE) options (mc.cores=parallel::detectCores ()) # Run on multiple cores set.seed (3875) ir <- data.frame (scale (iris[, -5]), Species=iris[, 5]) ### With improper prior it takes about 12 minutes, with about 40% CPU utilization and fans running, ### so you probably don't want to casually run the next line... system.time (b1 <- brm (Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width, data=ir, family="categorical", n.chains=3, n.iter=3000, n.warmup=600))
First, note that the
brm call looks like
glm or other standard regression functions. Second, I advised you not to run the
brm because on my couple-of-year-old Macbook Pro, it takes about 12 minutes to run. Why so long? Let’s look at some of the results of running it:
b1 # ===> Result, below Family: categorical (logit) Formula: Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width Data: ir (Number of observations: 150) Samples: 3 chains, each with n.iter = 3000; n.warmup = 600; n.thin = 1; total post-warmup samples = 7200 WAIC: NaN Fixed Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 1811.49 1282.51 -669.27 4171.99 16 1.26 Intercept 1773.48 1282.87 -707.92 4129.22 16 1.26 Petal.Length 3814.69 6080.50 -8398.92 15011.29 2 2.32 Petal.Length 3848.02 6080.70 -8353.52 15032.85 2 2.32 Petal.Width 14769.65 18021.35 -2921.08 54798.11 2 3.36 Petal.Width 14794.32 18021.10 -2902.81 54829.05 2 3.36 Sepal.Length 1519.97 1897.12 -2270.30 5334.05 7 1.43 Sepal.Length 1515.83 1897.17 -2274.31 5332.95 7 1.43 Sepal.Width -7371.98 5370.24 -18512.35 -935.85 2 2.51 Sepal.Width -7377.22 5370.22 -18515.78 -941.65 2 2.51
A multinomial logistic regression involves multiple pair-wise logistic regressions, and the default is a baseline level versus the other levels. In this case, the last level (
virginica) is the baseline, so we see results for 1)
virginica, and 2)
brms provides three other options for ordinal regressions, too.)
The first “diagnostic” we might notice is that it took way longer to run than we might’ve expected (12 minutes) for such a small dataset. Turning to the formal results above, we see huge estimated coefficients, huge error margins, a tiny effective sample size (2-16 effective samples out of 7200 actual samples), and an
Rhatsignificantly different from 1. So we can officially say something (everything, actually) is very wrong.
If we were coding in
Stan ourselves, we’d have to think about bugs we might’ve introduced, but with
brms, we can assume for now that the code is correct. So the first thing that comes to mind is that the default flat (improper) priors are so broad that the sampler is wandering aimlessly, which gives poor results and takes a long time because of many rejections. The first graph in this posting was generated by
plot (b1), and it clearly shows non-convergence of
virginica). This is a good reason to run multiple chains, since you can see how poor the mixing is and how different the densities are. (If we want nice interactive interface to all of the results, we could
Let’s try again, with more reasonable priors. In the case of a logistic regression, the exponentiated coefficients reflect the increase in probability for a unit increase of the variable, so let’s try using a
normal (0, 8) prior, (95% CI is , which easily covers reasonable odds):
system.time (b2 <- brm (Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width, data=ir, family="categorical", n.chains=3, n.iter=3000, n.warmup=600, prior=c(set_prior ("normal (0, 8)"))))
This only takes about a minute to run — about half of which involves compiling the model in C++ — which is a more reasonable time, and the results are much better:
Family: categorical (logit) Formula: Species ~ Petal.Length + Petal.Width + Sepal.Length + Sepal.Width Data: ir (Number of observations: 150) Samples: 3 chains, each with iter = 3000; warmup = 600; thin = 1; total post-warmup samples = 7200 WAIC: 21.87 Fixed Effects: Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat Intercept 6.90 3.21 1.52 14.06 2869 1 Intercept -5.85 4.18 -13.97 2.50 2865 1 Petal.Length 4.15 4.45 -4.26 12.76 2835 1 Petal.Length 15.48 5.06 5.88 25.62 3399 1 Petal.Width 4.32 4.50 -4.50 13.11 2760 1 Petal.Width 13.29 4.91 3.76 22.90 2800 1 Sepal.Length 4.92 4.11 -2.85 13.29 2360 1 Sepal.Length 3.19 4.16 -4.52 11.69 2546 1 Sepal.Width -4.00 2.37 -9.13 0.19 2187 1 Sepal.Width -5.69 2.57 -11.16 -1.06 2262 1
Much more reasonable estimates, errors, effective samples, and
Rhat. And let’s compare a plot for
Ahhh, that looks a lot better. The chains mix well and seem to randomly explore the space, and the densities closely agree. Go back to the top graph and compare, to see how much of an influence poor priors can have.
We have a wide choice of priors, including Normal, Student’s t, Cauchy, Laplace (double exponential), and many others. The formula argument is a lot like
lmer‘s (from package
lme4, an excellent place to start a hierarchical/mixed model to check things before moving to Bayesian solutions) with an addition:
response | addition ~ fixed + (random | group)
addition can be replaced with function calls
trunc, to specify SE of the observations (for meta-analysis), weighted regression, to specify the number of trials underlying each observation, the number of categories, and censoring or truncation, respectively. (See details of
brm for which families these apply to, and how they are used.) You can do zero-inflated and hurdle models, specify multiplicative effects, and of course do the usual hierarchical/mixed effects (
group) as well.
acat options provide different options than the default baseline model for (ordinal) categorical models.)
All in one function! Oh, did I mention that you can specify AR, MA, and ARMA correlation structures?
The only downside to
brms could be that it generates Stan code on the fly and passes it to Stan via
rstan, which will result in it being compiled. For larger models, the 40-60 seconds for compilation won’t matter much, but for small models like this Iris model it dominates the run time. This technique provides great flexibility, as listed above, so I think it’s worth it.
Another example of
brms‘ emphasis on flexibility is that the Stan code it generates is simple, and could be straightforwardly used to learn Stan coding or as a basis for a more-complex Stan model that you’d modify and then run directly via
rstan. In keeping with this,
brms provides the
make_standata functions to allow direct access to this functionality.
rstanarm package takes a similar, but different approach. Three differences stand out: First, instead of a single main function call,
rstanarm has several calls that are meant to be similar to pre-existing functions you are probably already using. Just prepend with “stan_”:
stan_polr (Ordinal Logistic). Oh, and you’ll probably want to provide some priors, too.
rstanarm pre-compiles the models it supports when it’s installed, so it skips the compilation step when you use it. You’ll notice that it immediately jumps to running the sampler rather than having a “Compiling C++” step. The code it generates is considerably more extensive and involved than the code
brms generates, which seems to allow it to sample faster but which would also make it more difficult to take and adapt it by hand to go beyond what the
rstanarm/brms approach offers explicitly. In keeping with this philosophy, there is no explicit function for seeing the generated code, though you can always look into the fitted model to see it. For example, if the model were
br2, you could look at
Third, as its excellent vignettes emphasize, Bayesian modeling is a series of steps that include posterior checks, and
rstanarm provides a couple of functions to help you, including
pp_check. (This difference is not as hard-wired as the first two, and
brmscould/should someday include similar functions, but it’s a statement that’s consistent with the Stan team’s emphasis.)
As a quick example of
rstanarm use, let’s build a (poor, toy) model on the
mm <- stan_glm (mpg ~ ., data=mtcars, prior=normal (0, 8)) mm #===> Results stan_glm(formula = mpg ~ ., data = mtcars, prior = normal(0, 8)) Estimates: Median MAD_SD (Intercept) 11.7 19.1 cyl -0.1 1.1 disp 0.0 0.0 hp 0.0 0.0 drat 0.8 1.7 wt -3.7 2.0 qsec 0.8 0.8 vs 0.3 2.1 am 2.5 2.2 gear 0.7 1.5 carb -0.2 0.9 sigma 2.7 0.4 Sample avg. posterior predictive distribution of y (X = xbar): Median MAD_SD mean_PPD 20.1 0.7
Note the more sparse output, which Gelman promotes. You can get more detail with
summary (br), and you can also use
shinystan to look at most everything that a Bayesian regression can give you. We can look at the values and CIs of the coefficients with
plot (mm), and we can compare posterior sample distributions with the actual distribution with:
pp_check (mm, "dist", nreps=30):
I could go into more detail, but this is getting a bit long, and
rstanarm is a very nice package, too, so let’s wrap things up with a comparison of the two packages and some tips.
Commonalities and other differences
Both packages support a wide variety of regression models — pretty much everything you’ll ever need. Both packages use Stan, via
shinystan, which means you can also use
rstan capabilities as well, and you get parallel execution support — mainly useful for multiple chains, which you should always do. Both packages support sparse solutions,
brms via Laplace or Horseshoe priors, and
rstanarm via Hierarchical Shrinkage Family priors. Both packages support Stan 2.9’s new Variational Bayes methods, which are much faster then MCMC sampling (an order of magnitude or more), but approximate and only valid for initial explorations, not final results.
Because of its pre-compiled-model approach,
rstanarm is faster in starting to sample for small models, and is slightly faster overall, though a bit less flexible with things like priors.
brms supports (non-ordinal) multinomial logistic regression, several ordinal logistic regression types, and time-series correlation structures.
rstanarm supports GAMMs (via
rstanarm is done by the Stan/
make_stancode makes Stan less of a black box and allows you to go beyond pre-packaged capabilities, while
pp_check provides a useful tool for the important step of posterior checking.
Bayesian modeling is a general machine that can model any kind of regression you can think of. Until recently, if you wanted to take advantage of this general machinery, you’d have to learn a general tool and its language. If you simply wanted to use Bayesian methods, you were often forced to use very-specialized functions that weren’t flexible. With the advent of
rstanarm, R users can now use extremely flexible functions from within the familiar and powerful R framework. Perhaps we won’t all become Bayesians now, but we now have significantly fewer excuses for not doing so. This is very exciting!
First, Stan’s HMC/NUTS sampler is slower per sample, but better explores the probability space, so you should be able to use fewer samples than you might’ve come to expect with other samplers. (Probably an order of magnitude fewer.) Second, Stan transforms code to C++ and then compiles the C++, which introduces an initial delay at the start of sampling. (This is bypassed in
rstanarm.) Third, don’t forget the
options statements I started with: you really need to run multiple chains and the fastest way to do that is by having Stan run multiple processes/threads.
Remember that the results of the
stan_ plots, such as
stan_dens or the results of
plot (mod, "dens") are
ggplot2 objects and can be modified with additional
geoms. For example, if you want to zoom in on a density plot:
stan_plot (b2$fit, show_density=TRUE) + coord_cartesian (xlim=c(-15, 15))
(Note: you want to use
coord_cartesian rather than
xlim which eliminates points and screws up your plot.) If you want to jitter and adjust the opacity of
pp_check (mm, check="scatter", nreps=12) + geom_jitter (width=0.1, height=0.1, color=rgb (0, 0, 0, 0.2))