The problem:
When fitting multilevel models, the group-level variance parameters can be difficult to estimate. Posterior distributions are wide, and point estimates are noisy. The maximum marginal likelihood estimate of the variance parameter can often be zero, which is a problem for computational algorithms such as lme4 which are based on this marginal mode. For models with multiple varying coefficients (varying-intercept, varying-slope models), the bigger the group-level covariance matrix, the more likely it is that its max marginal likelihood estimate will be degenerate. This leads to computational problems as well as problems with the estimated coefficients, as they get not just partially pooled but completely pooled toward the fitted model.
The solution:
Priors. Zero-avoiding or boundary-avoiding priors to avoid zero or degenerate estimates of group-level variances, also informative priors to get more reasonable estimates when the number of groups is small.
The research papers:
[2013] A nondegenerate estimator for hierarchical variance parameters via penalized likelihood estimation. {\em Psychometrika} {\bf 78}, 685–709. (Yeojin Chung, Sophia Rabe-Hesketh, Andrew Gelman, Jingchen Liu, and Vincent Dorie)
[2014] Weakly informative prior for point estimation of covariance matrices in hierarchical models. {\em Journal of Educational and Behavioral Statistics} {\bf 40}, 136–157. (Yeojin Chung, Andrew Gelman, Sophia Rabe-Hesketh, Jingchen Liu, and Vincent Dorie)
The R package:
blme: Bayesian Linear Mixed-Effects Models, by Vince Dorie, Doug Bates, Martin Maechler, Ben Bolker, and Steven Walker
Going forward:
blme is great but we’d also like to have full Bayes. Stan does full Bayes but can be slow if you have a lot of data and a lot of groups. Just for example, suppose you have longitudinal data with 5 observations on each of 100,000 people. Then a hierarchical model will have hundreds of thousands of parameters—that’s a lot! On the other hand, the Bayesian central limit theorem should be working in your favor (see appendix B of BDA, for example). So some combination of approximate and full Bayesian inference should work.
Also, lme4, and even blme, can have trouble when you have lots of variance parameters running around, and lme4 has its own issues, which unfortunately blme inherits, regarding computation with empty groups and various issues like that which should not really be a problem with Bayesian inference with informative priors.
Right now, though, we don’t have this best-of-both-worlds Bayesian solution that does full Bayes when computationally feasible and uses appropriate approximations otherwise. So blme is part of our toolbox. Thanks to Vince!
Perhaps a redundancy or error: Under “The research papers”, the 2014 Chung et al paper is listed twice
Fixed; thanks.
blme looks like blame when you glance over it fast.
Andrew, do you think that blme fits a different niche than brms, and if so, would you mind highlighting where it may be preferred?
Kiran,
I’d rather use rstanarm or brms because they do full Bayesian inference, or code the model straight up in Stan which gives me unlimited ways to expand the model (for example, you can take a logistic regression and have the probabilities range from 0.01 to 0.99 instead of from 0 to 1, so as to allow for occasional outliers).
The only trouble with rstanarm, brms, or full Stan is that it can be slow when you have a lot of data, and blme can be faster. So that’s when I’d use blme, is when I want an answer faster than I can get with full Bayes.
Someone said this: if it’s worth doing, it’s worth doing slowly. But more seriously, when working exploratievly, I sometimes start with lmer for this reason, to get an idea of what the MLEs will be like and which variance components will need regularization, before I go into full Bayesian modeling. I think I got that workflow from reading Gelman and Hill 2007.
Shravan:
We call it that trail of breadcrumbs or Bayesian workflow.
I still have issues with MLE in hierarchical models with thousands of groups, even with zero avoiding priors. I’ve tried going beyond zero avoiding priors like the suggested log gamma and gone all the way to strongly biased ones like log normal. In the end, I ended up hard coding a loose variance as data, then inferring variance empirically, like a single-iteration EM type algorithm. Not pretty, but it worked well enough.
Maybe I needed tighter control over my residuals, so the hierarchical parameters didn’t have an escape route by blowing the predictions up?
I am aware that Integrated Nested Laplace Approximation is a valid alternative to estimate group-level variance parameters when you have a lot of data (https://www.r-bloggers.com/2018/11/lmer-vs-inla-for-variance-components/). Is there any theoretical/statistical reason why you do not suggest using it?
In addition to the distinction between (a) MAP (blme) vs (b) full-Bayes (stan), I think there is an orthonormal companion: there are generally two ways to tune group variance parameter: either by (i) Bayes which effectively weights the parameter by the marginal likelihood, or (ii) some sort of cross validation which seeks the optimal amount of shrinkage. Can we use (ii) to construct better default prior for (i)? Can we use (i) to better stabilize (ii)? There are many open questions even in the simplest 8 school model.
Actually I often find that blme + arm::sim() produce posteriors that are very similar to full stan sampling, at least for linear models.
It would be great if these two technologies could be used as a possible backend for fast inference with stan and, indirectly, for brms.
I would already use it as drop-in solution for brms when exploring model specifications, but unfortunately the prior syntax and choices are very different…