Carlos Parada writes:
If I do a spike-and-slab regression in Stan (or Turing.jl, I guess; I don’t believe Stan allows for spike-and-slab priors), will that give equivalent results to doing Bayesian Model Averaging of two models, one of which includes a predictor variable and one of which doesn’t? (Assuming equivalent priors.) Mostly I’m somewhat concerned I picked the wrong side after I learned that BMA isn’t asymptotically efficient. Or, less jokingly, I’ve learned to expect that any time a Bayesian procedure does worse than a frequentist procedure, it’s because you’re using the wrong Bayesian procedure. I’m not sure what the alternative Bayesian procedure is. (Leave-one-out cross-validation doesn’t count!)
I know that BMA will successfully maximize the log-probability assigned to the correct model, assuming the model is in the candidate set. I guess there’s no reason this has to produce equivalent outcomes to the ones I’d get by minimizing the expected K-L divergence. But the idea of averaging models using anything other than Bayesian updating still rubs me the wrong way. As the kind of person who came to Bayes for mostly philosophical reasons, it feels ugly — what happened to Savage’s Dutch Book arguments, or Cox’s theorem?
I have two responses. First, if you want to select a subset of predictors, I recommend the horseshoe model. Second, I recommend stacking rather than Bayesian model averaging; see this paper. Regarding Parada’s question: realistically, if you’re thinking about model averaging at all, you’re going to be in the M-open setting in which you don’t think that any of the candidate models is a good description of the underlying process. As we discuss in the stacking paper, so-called BMA is designed for the M-closed setting in which one of the candidate models is correct, and it doesn’t always work so well in M-open.
As I mentioned in the Yuling Yao thread, you have to be clear on whether your goal is parameter estimation or prediction.
We know the equation for volume of a box: V = L * W * H
But lets say we didn’t. We could collect all sorts of data about boxes. The color, the weight, the material, what it contains, etc. You can then use a model (or ensemble of them) based on these variables to get decent skill on predicting the volume, but none of the coefficients actually mean anything.
> realistically, if you’re thinking about model averaging at all, you’re going to be in the M-open setting in which you don’t think that any of the candidate models is a good description of the underlying process.
I don’t think this is true at all.
Suppose you’re looking at hyperspectral images of agricultural land. You think a plot of land can be either, corn, soybeans, wheat, barley, rye, or a sample from a broad-spectrum of stuff you’ve seen before that isn’t any of those things.
There’s lots of cases like this. Radar reflection from airplanes, bacterial metabolites, tomography of subterranean rock formations… whatever.
Bayesian model averaging works well to both classify data into which type of thing generated it, and give uncertainty among the possibilities.
Daniel:
I think of your example as a “classification” problem, not a “model averaging” problem. That is, I don’t think of corn, soybeans, etc., as different “models.” I understand that a model averaging problem is mathematically the same as a classification problem, but in practice they are different, for reasons discussed in our paper.
Thanks, it’s good to understand if there’s a terminology difference. I agree that this is similar to a classification problem but it’s not hard to turn it into something like “we have a blurry picture of a few acres of multi-crop fields and we want to infer the proportion of corn, soybeans, and squash.”
In this context “Corn” is in fact a model, specifically because you’re getting absorption and reflectance spectrum over a few hundred frequency bands across a few thousand nanometers, you have a mechanistic model of what the output from “corn” should be which is different from what “soybeans” is because there are different pigments in the leaves and such.
Daniel, perhaps you are thinking about some naive Bayes classifier, which mathematically resembles BMA?
Well just imagine f(nu) is a measured spectrum of brightness in a variety of wavelengths. And we want to model it as
f(nu) = fc(nu) *pc+fs(nu)*ps+fo(nu)*po+e(nu)
And pc+ps+po=1 and fc is a corn spectrum, fs is a soybean spectrum, and fo is an “other” spectrum, and we have some data with ground truth to learn those spectra simultaneously with learning the pc,ps,po weights for each photograph of hundreds of locations.
If there is no extra error term e(nu), then the model you have is the naive Bayes classifier with x=nu and y=category of corn/soybean/other. We could also say it is BMA with three modes over three models of x.
For people who wish they were able to use spike-and-slab priors in Stan, I don’t think the horseshoe priors are the best substitute. Horseshoe priors can specialize to a spike or a slab but those are each limiting cases as the local scale approaches zero or infinity. So, Stan is not going to be able to reach either of those in the unconstrained space and you will get divergent transitions and / or hit the maximum treedepth if you try. For finite values of the local scale, the horseshoe priors have a mode with finite density at zero, rather than a spike.
The closest thing to a spike-and-slab prior that is viable in Stan is a product normal distribution
https://en.wikipedia.org/wiki/Distribution_of_the_product_of_two_random_variables#Independent_central-normal_distributions
whose PDF is proportional to a Bessel K function in the two component case, which is infinite at zero. This singularity would be a disaster for Stan if you tried to use the log-kernel directly, but multiplying two (or more) parameters with mean-zero normal priors together to define the coefficient in the transformed parameter block works “well” (there is a label-switching problem with the primitive components so different chains will go to different modes, but the effective sample size of their product should be acceptable). With n components, the asymptotic (as the product tends toward zero or infinity) densities are given in equation 1.5 of
https://arxiv.org/pdf/1309.4344.pdf
Apparently, the product-normal distribution is an extreme special case of some other distribution that was derived by some Bayesian whose last name I have momentarily forgotten (it starts with P I think) that also has infinite density at zero.
In Julia you can do Gibbs sampling of discrete parameters while doing NUTS on continuous ones. Not that this will sample super well, but you could construct a model where it’s something like:
a ~ Bernoulli(0.5)
b ~ Normal(0,1.0)
c = a * b
Now c is either 0, or a normal distribution.
To make this sample better it might be useful to do something else. In particular it might make sense to use a continuous function to map something like Normal(0,1) to a different variable where f(x) = 0 for x between -a,a such that the probability to have f(x) = 0 is a known value, but it can continuously vary away from zero.
Thinking about this, suppose you take uniform(-1,1) and shove it through a hinged 3 section function. Lets say its like three lines connected by smooth hinges. Between (-1,-a) it goes from say -1e6 to 0, then its zero from (-a,a) then goes from(0,1e6) between (a,1)
Now a is the probability of being 0 and (1-a) is the probability of being anything between (-1e6,1e6) other than 0.
See Andrews “continuous hinge” post from a few years back to figure out how to make this work
Some controversial views to throw in:
1) We could replace the taxonomy of “M-open, M-closed and M-complete” by the degree of model misspecification, which has a probability 0 to be exact 0.
2) The distinction between “our goal is prediction” and “our goal is estimation” is misleading. When we are interested in parameter estimation, we are interested in causal inference, which is some sort of prediction. But then a prediction without causal inference is not robust, so in the end there is only one goal.
For the question in OP: “If I do a spike-and-slab regression in Stan, will that give equivalent results to doing (a) Bayesian Model Averaging, (b) stacking of two models (c) mixture of expert, one of which includes a predictor variable and one of which doesn’t?”
The answer is NO in all three cases. Model averaging is not the same as model extension. We recommend to use model averaging when you have already exhausted your energy in individual modeling, rather than using model averaging as a short cut for model building.
The only time to care about parameter estimation is if you derived the model from a set of well-defined assumptions that someone thinks is plausible.
Eg, the number of cells after D divisions could be modeled as:
N = N_0 * 2^D
where: D = r*t
So number of cells at time t would be
N(t) = N_0 * 2^(r*t)
This is overly simplified and so wrong, but still the parameters (N_0 = initial number of cells, r = rate of division) have well defined meaning.
Instead you could do:
N(t) = a*N_0 + b*t
The a and b coefficients do not mean anything because this model was not derived from assumptions about what is going on.
I wouldn’t be willing to say that this is the only time one might care about the specific values/ranges of parameters – mostly out of ignorance though I cannot conjure an counterexample at present – but in general seems to me to be the more correct way of thinking.
I forget who Meehl was quoting but in his 1989 lecture series he quotes something to the following effect: “the existence of twilight doesn’t preclude noon and midnight from being tolerably distinct”
Of course there is a distinction between prediction only models and those trying to derive parameters with some semblance of causal meaning – even if the latter case can lead quite nicely to the former.
“If I do a spike-and-slab regression… will that give equivalent results to doing Bayesian Model Averaging of two models, one of which includes a predictor variable and one of which doesn’t?”
The answer is YES and we explain why this is with a simple coin flipping example here:
https://arxiv.org/pdf/2204.06054.pdf
It would be really great to apply stacking to a similarly simple toy example to help readers better understand the mechanics (and advantages) of this approach.