1
$\begingroup$

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.

$\endgroup$
8
  • $\begingroup$ Is there a specific reason why you have the impression that brms should run faster than functions from other packages? There is quite a bit happening under the hood of a brms model because it has to compile a Stan model and sample from it. (BTW, if that last sentence solves the problem, I'm happy to turn it into an official answer.) $\endgroup$ – Sean Raleigh Jun 29 '20 at 4:27
  • $\begingroup$ @SeanRaleigh Thank you very much for the quick response. I didnt expect brms to run significantly faster than the frequentist models, but when I see it run a simple intercept-only model much more slower than the frequentist intercept-only model, I was thinking if I do anything wrong. Usually, the application of MLM involves level-1 or level-2 covariates, sometimes even with cross level interactions. Does the specification of prior make it faster? or is any other way to make it faster (on a same machine)? $\endgroup$ – user11806155 Jun 29 '20 at 5:20
  • $\begingroup$ There are always little things you can do to speed things up, including some tricks with priors. But what you're seeing is a sort of floor effect...it just takes Stan some time to spin up, even for simple models (at least the first time you run the code). The payoff comes for much larger and more complex models when the frequentist packages struggle to get any kind of good estimate, but the Bayesian models have less trouble. Actually, for me, the best payoff is the philosophical one that comes from doing inference sensibly, but that's a topic for a different thread. ;-) $\endgroup$ – Sean Raleigh Jun 29 '20 at 5:43
  • $\begingroup$ I'll also add that the data set is somewhat large-ish (over 7000 observations). $\endgroup$ – Sean Raleigh Jun 29 '20 at 5:46
  • $\begingroup$ @SeanRaleigh Dear Prof. Raleigh Thank you very much for the suggestions. In this case I think i need to think about other options (e.g., upgrade my hardware...), because this dataset is only for reproducible purpose. The real analysis will be performed on an industry-grade dataset, which is larger... $\endgroup$ – user11806155 Jun 29 '20 at 5:50
4
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ just a follow-up, is it possible to have a frequentist-equivalent (to some extent) bayesian model, for instance, to set chain=1, iter=1, warmup=0, or prior= something that is essentially not a distribution but a number? My goal is to do multilevel multinomial logistic regression, but I could not find a frequentist R package to do it, thats why I turn to brms. I dont mind not getting bayesian inference, I just need a solution to run the multilevel multinomial logistic model in R. $\endgroup$ – user11806155 Jul 2 '20 at 14:28
  • $\begingroup$ I mean, to trick the brm() function to think that I am running a bayesian model, when I am actually running a frequentist model. $\endgroup$ – user11806155 Jul 2 '20 at 14:43
  • $\begingroup$ I mean, to trick the brm() function to think that I am running a bayesian model, when I am actually running a frequentist model. $\endgroup$ – user11806155 Jul 2 '20 at 14:43
  • $\begingroup$ To my knowledge, you cannot trick brms into giving you frequentist inference. The number of chains and all that is irrelevant because MCMC algorithms are not inherently Bayesian or frequentist; they just sample from a probability distribution. The "guts" of a brm call sample from a posterior distribution and there's no way around that. In some cases, flat uniform priors give answers similar to frequentist answers, but I don't know that it's the case for multilevel models. And you're still not getting a frequentist answer...just a Bayesian one that kind of agrees with a frequentist one. $\endgroup$ – Sean Raleigh Jul 3 '20 at 15:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.