brms intercept only model runs very slow

Question

I am trying to learn brms package for multilevel modeling. A reproducible code is as below:

library(MCMCglmm)
library(nlme)
library(lme4)
library(brms)

data(MathAchieve,package='nlme')
data(MathAchSchool,package='nlme')
dat=merge(MathAchSchool,MathAchieve,by='School')
str(dat)

set.seed(23429989)
m1=brm(MathAch~1+(1|School),dat)
m2=lmer(MathAch~1+(1|School),dat)

system.time(lmer(MathAch~1+(1|School),dat))
   user  system elapsed 
  0.054   0.002   0.056 

system.time(brm(MathAch~1+(1|School),dat))
   user  system elapsed 
124.646   5.787 143.627 

While the model seems run properly, it runs much slower than a frequentist MLM model (e.g., using lmer()) even for an intercept-only model. What did I miss? Thank you very much.

Answer 1

Under the hood, the brms package builds a Stan model. There are two things that happen that take some time. First, Stan compiles some C++ code. After that, Stan runs an MCMC (Markov Chain Monte Carlo) algorithm that draws samples from the posterior distribution. (The actual details are way more complicated than that, but I've tried to capture the essence in a nutshell.)

Both steps can take quite a bit of time, depending on the number of parameters in the model and the number of data points.

If you end up running your model more than once, you can save the compiled Stan program using

rstan_options(auto_write = TRUE)

Or you can use the update method in brms to accomplish something similar.

As for drawing samples from the posterior, you can run multiple chains at the same time using parallel processes. But each chain can still take a long time when the model and/or the data is large. (In the posted example, the model is very simple, but the dataset is large-ish with around 7000 observations.)

There are lots of little tricks for re-parameterizing models and sometimes they speed things up quite a bit. For the simple model posed in the question, there's not really much more to do, but you could read up on non-centered parameterization for one of the best and easiest-to-implement tricks.

The benefit of using Stan over other fitting methods (IMO) is in fitting complex models. Functions like lmer don't have a syntax that can accommodate all the possible complexity in certain models, whereas in Stan, there is no theoretical limit to the complexity. (Of course, brm has some of the same limitations as lmer due to the concise formula syntax.) Even when frequentist algorithms can process a model, their mode of inference does not match the way a Bayesian model propagates uncertainty through all the parameter space. I won't get into all the benefits of being Bayesian here, but if you care about inference, it's worth thinking about using tools that give you results that can be interpreted and explored in helpful ways.