Kiran Gauthier writes:
After attending your talk at the University of Minnesota, I wanted to ask a follow up regarding the structure of hierarchical / multilevel models but we ran out of time. Do you have any insight on the thought that probabilistic programming languages are so flexible, and the Bayesian inference algorithms so fast, that there is a balance to be struck between “simple” hierarchical models and more “complex” hierarchical models that augment the simple frameworks with more modeled interactions when analyzing real data?
I think that a real benefit of the Bayesian paradigm is that (in theory) if the data doesn’t converge my uncertainty in a parameter, then the inference engine should return my prior (or something close to it). Does this happen in reality? I know you’ve written about canary variables before as an indication of model misspecification which I think is an awesome idea, I’m just wondering how to strike that balance between a simple / approximate model, and a more complicated model given that the true generative process is unknown, and noisy data with bad models can lead good inference engines astray.
My reply: I think complex models are better. As Radford Neal put it so memorably, nearly thirty years ago,
Sometimes a simple model will outperform a more complex model . . . Nevertheless, I believe that deliberately limiting the complexity of the model is not fruitful when the problem is evidently complex. Instead, if a simple model is found that outperforms some particular complex model, the appropriate response is to define a different complex model that captures whatever aspect of the problem led to the simple model performing well.
That said, I don’t recommend fitting the complex model on its own. Rather, I recommend building up to it from something simpler. This building-up occurs on two time scales:
1. When working on your particular problem, start with simple comparisons and then fit more and more complicated models until you have what you want.
2. Taking the long view, as our understanding of statistics progresses, we can understand more complicated models and fit them routinely. This is kind of the converse of the idea that statistical analysis recapitulates the development of statistical methods.
Hi Andrew,
Thank you very much for your reply. I think this all fits with the larger paradigm of the iterative Bayesian model building / model checking workflow and your recent post on breaking models to grok them.
Something that I’ve thought a lot about but never formalized is the idea of approximate inference for less important quantities to focus the compute on the latent variables of interest.
An example to make it more tangible, say I have k observations for j individuals and I want to summarize the individual parameters using a multivariate Gaussian instead of doing full Bayesian inference. My inference of the population level parameters (i) should be pretty ok as long as my Laplace approximation is a good one and I’ve probably saved a lot of compute.
Has this been implemented? I think I’ve seen a (two-stage / one-stage?) estimation for fixed- and random-effects in multilevel models but I haven’t dug into the details yet.
If not, I can try to take a crack at it so that we can try to get some answers on how useful it’ll be relative to other tricks for hierarchical models (narrowing priors, centered parametrization). I did something similar during my PhD.
Maybe a bit far afield, but John Kruschke (author of the highly recommended “Doing Bayesian Data Analysis”, or as I call it, the “puppy book”) actually proposed a modeling framework similar to what you describe. Only it was not intended as a statistical model, but as a cognitive model of associative learning! Kruschke called it “Locally Bayesian Learning”. The idea is that there are multiple layers arranged in a hierarchy. Each individual layer does full Bayesian inference based on its immediate inputs, but then only passes the marginal means as input to the next layer.
You could imagine that, instead of sensory inputs, the different “layers” of the model could correspond to different levels of grouping (e.g., individuals, counties, etc.). And I suppose the “output” of each layer need not be the marginal—it could be the MAP estimate too.
Anyway, here’s a link to Kruschke’s 2006 paper:
https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=55cf38dd71586e3ec567dd1dac5ba7ff291be3d5
gec, this is awesome, thank you.
@Kiran: What you’re describing sounds a lot like nested Laplace approximations, which are the basis of INLA. Or more generally, empirical Bayes. You see the same thing in frequentist settings with something like lme4 (the “lme” being for linear mixed effects).
Efficient and robust computation for black-box versions of this kind of algorithm are ongoing. Andrew even threw his hat into he ring at one point with an algorithm. Charles Margossian and Steve Bronder are scheming to put something like this into Stan.
@gec: It sounds like you’re talking about what is called a “cut” in graphical model inference. It also arises from multiple imputation, where you impute and uncertainty flows forward but not back into your imputation. In something like lme4, what you’re calling “output” is done through optimization. First it marginalizes, then optimizes the marginalized density—that’s why it’s called “maximum marginal likelihood”.
Ironically, neither empirical Bayes nor cut are full Bayesian inference methods.
Hi Bob,
Thanks for letting me know, I’ve been aware of INLA but I’ll look more deeply into it. I’ve always been a bit suspicious of empirical Bayes methods because I think they’re using optimization methods instead of samplers under the hood, and that unconstrained optimization leads to weird minima in the hierarchical case from the examples I’ve worked on.
I guess with tight-ish priors and running a MAP-finding optimization, and then approximating the Hessian as a multivariate Normal approximation for the individual parameters you could get the best of both worlds. Hopefully INLA over the individual (random-effects) is what I’m looking for, and I can adapt that to do full Bayesian inference on the population parameters of interest.
Empirical Bayes typically runs optimization on the marginalized likelihood. The maximum (marginal) likelihood is usually well defined in these settings. It’s not trying to run optimization on the hierarchical model. An example is lme4 in R. In fact, the resulting marginalized density usually has better geometry—that was the upshot of (part of) Charles Margossian’s Ph.D. thesis and the basis for INLA.
You can use zero-avoiding priors to make MAP estimates well defined in hierarchical models, but Andrew moved away from that approach about 10 years ago in favor of zero-embracing priors like half-normal because the zero-avoiding estimate is typically not what you want.
Could you recommend any readings that would help me understand the difference between a multivariate Gaussian model and a “full” Bayesian model?
“Full Bayes” usually means integrating over the posterior as opposed to, for example, taking an optimization step like empirical Bayes.
A multivariate Gaussian is a density, which I would not call a “model”. A Bayesian model typically has parameters and a data-generating process. If the model is literally just y ~ normal(mu, Sigma), then full Bayesian estimation for mu is hat(mu) = INTEGRAL_{mu, Sigma} mu * p(mu, Sigma | y) d(mu, Sigma).
I’m really surprised Andrew didn’t go where I was thinking he would with this: the unfolding flower. In an earlier post, he says:
This is what Andrew originally hired Matt Hoffman and me to work on—a language for specifying regressions that was more flexible than lme4 (we gave up and generalized even further, but not as usefully for regression, with Stan).
The issue you run into is that if you don’t have very much data, it’s hard to fit a complicated model. We’d like to have priors so that more complicated models reduce to simpler models when there is not enough data to support them. This is the basis for my co-favorite idea in Bayesian priors along with weakly informative priors: penalized complexity. It’s perhaps not surprising it was introduced by the INLA developers—they deal with a lot of complex spatio-temporal data with natural hierarchical structure (though space and time together are tricky as there are two ways to go from heterogeneous to spatial + temporal).
That reduction to simpler models is also something I’ve been super interested in, I tried to do it using spike-and-slab or horseshoe priors in PyMC a couple of years ago. I remember that the posterior geometry got pretty ugly and that my sampler had a hard time navigating it but I loved the conceptual idea of a sampler shutting off a parameter with some probability if the data doesn’t support it’s inclusion. I think Aki wrote about this recently as well, but this gives me even more hope that I can benefit from the INLA community.
Thanks for the references, Bob.
One useful way to think about the unfolding from simple to complex models is that the simple model, typically a variant of naive OLS, is a baseline against which increasingly complex models can be evaluated.