Jonah Gabry asks:
What are your latest thoughts on how to use a posterior from a previous fit as the prior for a subsequent fit after collecting new data? This came up recently when I was talking on the phone with some scientists in Utah who are devoted Stan users.
My reply:
The basic model is p(theta | y1, y2) is proportional to p(theta | y1) * p(y2 | theta), so that you use the posterior given y1 as the prior for the analysis of y2. In that case, it’s just a computational challenge, two which there are two obvious solutions: (a) Run Stan on p(theta | y1), get your 4000 simulations or whatever, and approximate them as a mixture of multivariate normals (there are various automated ways of doing this, and it shouldn’t really matter so much exactly how you do it), and then use that as your prior at the next step; or (b) Just rerun your entire analysis on (y1, y2) and cut out the middleman. Why not just do (b)? Maybe it’s a computational issue that you don’t want to re-fit the whole model. In that case you could do importance reweighting. Another trick is that if the model has hierarchical structure, you can do some factorization to get conditional posteriors; that’s what we did in our sports consulting project when the team wanted updated estimates after every game or every week and didn’t want to be re-fitting each time.
More generally, though, I’m wondering if maybe you really want a model with p(y1 | theta1) and p(y2 | theta2)? The point here is that once you’re talking about breaking up your analysis in the middle, maybe you should think about your parameters changing. In that case, you need a joint distribution for theta1, theta2, or, to put it another way, some model for how theta1 and theta2 might differ.
Good luck!
Suppose you are following a Bayesian sequential experimental design, you learn a posterior with your data, and then use this posterior to simulate unseen data using your learnt generative model. Select the experiment that maximizes a given objective, perform the experiment and then use your posterior as the prior for performing the above again.
Of course, you can just input the result of your experiment to get your new posterior, but for how many iterations is it feasible to relearn the whole model? I see it as preferable to distribution estimation, but still, there should be more efficient ways to go about it, right?
> really want a model with p(y1 | theta1) and p(y2 | theta2)?
You need to choose wisely what you accept as common. And there is much more discussion of that here – https://www.stat.columbia.edu/~gelman/research/unpublished/Amalgamating6.pdf
Briefly, “Statistical science has evolved from the growing awareness, extraction, and assessment of commonness in the midst of diversity. … Awareness of commonness can lead to an increase in evidence regarding the target; disregarding commonness wastes evidence; and mistaken acceptance of commonness destroys otherwise available evidence. It is the tension between these last two processes that drives many of the theoretical and
practical controversies within statistics.”
could you point us to a mwe that illustrates these options? in my experience, in hierarchical data, if we fit y1, y2 together as opposed to separately, the posteriors are quite different. reason: the variance components’ posteriors can change quite a bit depending on how one does the analysis.
The math – p(theta | y1, y2) is proportional to p(theta | y1) * p(y2 | theta) – presumes it is the same data generating model and prior – so if you are changing these all bets are off but maybe you _should_ change the them.
Would you be willing to provide more details/code (or a pointer to a paper/textbook) for how to do something like you suggest for the sports consulting example? That sounds incredibly useful.
David:
Unfortunately, I never wrote this up! What I remember is that we tried several different strategies, and I think the one that worked the best took advantage of the hierarchical structure of the problem.
Why did the team want to avoid re-fitting the model each week?
Dan:
The model was kinda complicated so it took awhile to fit. Also, we’re always playing with the model and improving it in various ways, so re-fitting each version could take a long time.
As a philosopher of science, it’s really notable to me that the question is, roughly, “how do you implement Bayesian updating?” And the answer is “why would you bother with updating?”
Dan:
Just to be clear, my final proposal is not to avoid updating; it’s just to update using a hierarchical model. The simple updating rule assumes a constant theta; the hierarchical model allows theta to vary, but it’s still Bayesian updating. See this paper with Weber et al. for an example.
Yes I am confused also, and must be missing something or perhaps it all depends on the situation. Isn’t this what is done all the time in Kalman filters and more general state-space models?
Roy:
Yes, indeed, a time-series model for theta could be fine; it all just depends on the context. In my linked article, there’s no time series, there are just two different experiments. If there is a time-series structure, that’s cool too.
Might be wrong, but there is something folks seem to miss about likelihood and prior multiplication – the basic mechanics.
For instance, if you have a simple regression through the origin with known variance, there only is only one unknown parameter and the prior and individual observation likelihoods are all easily plotted one at a time.
The combined posterior ~ prior * like1 * like2 * … * likn so is invariant to the order of multiplication – do it anyway you want.
Now the same is true for with multiple parameters but the parameters do not become identified until a certain number of likelihoods are multiplied – so plots of individual observation likelihoods are not (easily) plottable. (There is a fudge here, in that you can use the full sample to “estimate” all the parameters but one to get plottable views of almost individual observation likelihood for that one parameter – see https://statmodeling.stat.columbia.edu/wp-content/uploads/2011/05/plot13.pdf )
In this blog post, I think Andrew is presuming everyone understands the mechanics and is pointing to the value of questioning the all too often assumption that the parameters are common across all observations (e.g. perhaps some vary by study).
p.s. I had learned this in my first course in statistics from Don Fraser and when clinicians asked me to analyze multiple studies in 1985 – it seemed almost obvious to me.
Keith:
This is related to our paper on expectation propagation, in particular figure 2 of that paper.
Thanks, I missed that paper.
There indeed — within the closed system of the experiment at-hand — is the recursive updating of a Gaussian density, one observation at a time; subject to a dynamical model of the linear (or linearized) “device” through which some hidden process travels just before it generates an observed datum.
But open the system up and ask, with what do I prime the pump to start everything going? Engineering specification on the noise in a channel? The result of somebody’s published experiment? The result of my last computer simulation run out to time T?
I think it makes more sense to make use of the whole procedure in this way: as a system for studying parametrically the dependence of the (evolution of) the output on the assumptions one makes about the input.
In other words, the “Prior” that one puts in to prime the pump represents a “knob” and one should vary the knob. What is the *range* of variation of that knob? Now that is where the link to the world comes in: Lab A found such and such and Lab B found such and such. So one somehow or other smoothly varies one’s input (PDF) to include A’s and B’s and everything else “in between”.