Bayesian model averaging without reversible-jump MCMC on highly-dimensional models

[Ben Lambert]

Hi,

I have currently fit about 7 different quite highly-dimensional hierarchical models for a particular problem using Stan; all the fits, as well as convergence diagnostics (Rhat, n divergent, and all the other Shiny Stan diagnostics) are ok. The individual models take about an hour to run individually.

The models each have different parameterisations, with differing numbers of parameters. I would ideally like to do some sort of mixture modelling, where I assign priors on the probability of each model a priori, but then let the data decide the ultimate weights given to each model. (I know that BDA section 7.5 suggests that it is preferable to do more of a continuous model-type expansion, but as far as I can tell this is not possible in this paper, since its point is to compare the results obtained across existent models in the literature, where there is a considerable lack of consensus on model choice).

Since each model has different numbers of parameters, I would normally need to do reversible jump MCMC here. Firstly, I want to check - is this possible in Stan? (I think from this discussion that it is not: https://groups.google.com/forum/?fromgroups#!searchin/stan-users/reversible$20jump$20mcmc/stan-users/CXAT_XUTfwE/Nxgqjm_Zx1AJ0).

Secondly, does anyone know of any approximate results that could be used to average across all the models? (I have kept the log-likelihood in order to calculate LOO-CV, and WAIC, and was wondering whether this could be used.) I suspect that there isn't really anything I can do here, other than coding up the models myself in another language, but I wanted to tap the community brain before I begin this (arduous) process.

Best,

Ben

[Andrew Gelman]

My advice:

Forget reversible jump.

Instead, fit the 7 models separately—actually, you’ve already done that, so you’ve done all the hard work.  For each of the 7 fitted models you get a predicted value for each data point.  Call these y_hat_{nk} for data points n=1,…,N and models k=1,…,7.   Now just create a new model that is a mixture of the 7, with mixture fractions lambda_1,…,lambda_7 which are simplex-constrained (i.e., nonnegative and sum to 1).  You can estimate the lambdas that minimize the sum of the squared prediction errors (that’s just constrained least squares, but it’s easiest to do via a simple Stan model.)

You should be good to go.

Please let us know how this works.

A

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

If we had "cut", you could do that model mixture optimization directly
in Stan in one go. But it's not a Bayesian procedure.

If you literally built up a mixture model in Stan, you'd get a different
result, because each mixture component could fit the data its responsible for
rather than making each component fit all the data equally weighted.

You want to check out the marginalizing discrete parameters section of
the manual to figure out how to write that new model Andrew's suggesting.
Normal with fixed scale gives you a density proportional to the sum of square
errors because we work on the log scale.

- Bob

[Andrew Gelman]

I’d just suggest doing directly what I suggested below:
- First fit the 7 models separately. Use these to get point predictions for each data point.
- Write a new Stan program with data vector[N] y; matrix[N,y] y_pred; and parameter simplex[7] lambda; to minimize the sum of squared prediction errors.
- Then feed in the data and the Nx7 matrix of predicted values and run Stan-optimize on this new program to estimate lambda.
- Voilà, your new model: a weighted avg of your original 7 models with weights lambda_1,…lambda_7.

No marginalizing discrete parameters is necessary here.

I agree that this is not the same as fitting a mixture model, but I think the above makes sense as a form of model averaging. As I said, I’m curious how it works out.

[Bob Carpenter]

Andrew's approach already marginalized the discrete
responsibility parameter into a simplex for you. :-)

To fit this as a mixture model in Stan, we'd need cut.
Otherwise, you won't get the same fit for the components as
you would by fitting them separately because the data's
implicitly weighted by responsibility in the mixture
formulation.

- Bob

[Roy]

Hi Andrew

That sounds like an interesting idea as it could be applied to any set
of fitted models, even Frequentist ones, right? Are there any references
in the literature of model averaging you are basing your suggestion on?
I'm just curious :)

Cheers

Roy

[Michael Weylandt]

Andrew's approach sounds like (model) stacking from the machine
learning literature.

See, e.g., Leo Brieman, "Stacked Regression". Machine Learning 24(1),
pp. 49-64 (July 1996) or Section 8.8 of the Elements of Statistical
Learning by Hastie, Tibshirani, and Friedman (free PDF available on
Tibshirani's website).

[Bob Carpenter]

I think there's a big gap between doing the proper mixture modeling,
Andrew's half-way solution that cuts the inference on each model from
the mixing, and the typical Bayesian model averaging that weights by
model posterior (an aggregated cut-based solution):

http://www.stat.colostate.edu/~jah/papers/statsci.pdf

Often in machine learning (as in the 4th author of the above
paper's work winning the Netflix Prize!) you optimize on held-out
or cross-validated data rather than in-training-sample data, often
running a logistic regression with further predictors for items in
addition to an intercept:

http://www.netflixprize.com/assets/ProgressPrize2007_KorBell.pdf

- Bob

P.S. By "cut" I just mean running a model and gathering its results to
use downstream rather than allowing "feedback" in the model. It's
a "feature" of BUGS :-)

[Andrew Gelman]

Bob: I don’t know why you keep bringing up “cut”—this is just a distraction. It would be no problem to just fit a mixture in Stan directly but I’m addressing a specific question which comes up a lot, which is that a person fits separate models A, B, C, …, and then wants weights for them.

[Andrew Gelman]

Aki and I have discussed this approach. I’ve talked about it a lot but have never done it, hence my curiouity about how well it might work on this particular example. I don’t recall seeing the method in the literature but it could well be there. It’s also related to a recent paper by Christian Robert and Judith Rousseau.

> On Feb 7, 2016, at 4:32 PM, Roy <royco...@gmail.com> wrote:
>

[Andrew Gelman]

I took a look at the references that Bob sent, and my mixture-modeling optimization does seem pretty much the same as the “stacking” idea from the statistics literature. I did some quick googling and also found this paper from 2003 that seems relevant to this discussion:
http://delivery.acm.org/10.1145/970000/964281/4-683-clarke.pdf
This idea seems so natural, I don’t know why I’d not thought of it before, or why it doesn’t seem to come up much in the Bayesian model averaging literature.
A

[Roy Costilla]

Hi Andrew, Bob and all
Thanks for your response on a spin-off question!
Cheers
Roy

[Dustin Tran]

@Andrew, the link you wrote is broken.

Why not just use the evaluation criteria (normalized) as weights? Since Ben already has all the numbers to calculate LOO-CV, it’d be easy to weight based on these estimates of the log pointwise predictive density. I have no idea if this works, maybe the scaling is off. However it already sounds less sensitive than using the marginal likelihoods/Bayes factors as the loss function to weight things on, which is the more traditional thing to do.

Dustin

[Andrew Gelman]

Here’s the reference:
http://dl.acm.org/citation.cfm?id=964281
What we’re talking about _does_ appear to be stacking, and I have no idea why this is not a standard approach in Bayesian statistics. And why is it not in chapter 7 of BDA3? I’ve seen so much wasted effort going into attempts to compute marginal likelihoods for Bayesian model averaging, and this idea seems much better.

I guess it didn’t used to be so feasible. Now that we have pointwise loo at our fingertips, this stacking thing seems like the right way to go.

I prefer the term “mixture model” to “stacking,” though.

[Michael Betancourt]

We need to stop having these discussions over the weekend! Mike West has been
doing some brilliant work that generalizes all of this as a hierarchical model over
models. It goes beyond the typical machine learning approaches and Bayesian
model comparison approaches into one coherent framework where you can choose
the mixture weights based on your own self-consistent criterion.

[Bob Carpenter]

I think overloading "mixture model" here's going to be confusing
(not that "stacking" really helps matters). A standard mixture
model puts a prior on the simplex and estimates it along with the
individual model's parameters. You'll get a different marginal
posterior for the model parameters that way than if you fit the
models separately and then combine them after the fact. At least
I'm pretty sure you will.

How about "pseudomixture"? Or "postmixture"?

- Bob

[Bob Carpenter]

> On Feb 7, 2016, at 11:14 PM, Michael Betancourt <betan...@gmail.com> wrote:
>
> We need to stop having these discussions over the weekend! Mike West has been
> doing some brilliant work that generalizes all of this as a hierarchical model over
> models.

Is it related to what Andrew and Sebastian have been doing
for PK/PD models? This is why I brought up the term Andrew
wants me to drop in the first place!

- Bob

[Dustin Tran]

Following Christian Robert and Judith Rousseau’s language, it’s a mixture model conditional on the hypotheses. The hypotheses are the estimated predictive posteriors. Mixture modeling makes sense in that it’s all about the decision to condition on the models or to infer them jointly with the mixture weights.

Dustin

[Andrew Gelman]

Bob:

No, it is a mixture model. It’s just not _full Bayes inference_ for the mixture model.

A

[Michael Betancourt]

Right, the difference is whether the mixture is on the likelihood (which can
be written immediately in Stan) or on the component posteriors (which
require approximations or extra machinery to get marginal likelihood estimates).

[Andrew Gelman]

That would be great. The whole “Bayesian model averaging” thing has been a useless time sink for decades.

> On Feb 7, 2016, at 11:14 PM, Michael Betancourt <betan...@gmail.com> wrote:
>

[Bob Carpenter]

Is the motivation for not doing full Bayes purely computational?
Or is it about fitting the models independently before mixing?

- Bob

[Andrew Gelman]

1. Computational

2. Robustness to specification of priors within the separate models. For problems with pure Bayes model averaging, see my 1995 paper with Rubin, also this is discussed a bit in chap 7 of BDA3.

[Bob Carpenter]

Thanks. 2 makes sense. The mixture weights are like data weighting --- that's
why the mixture component won't get same parameter estimates as full model.
And now that I think about priors, of course they'll be affected --- there's
less overall data weight, and the weight there is comes from where the model's
the best at fitting, so it's not even a random subsample.

- Bob

[Michael Betancourt]

> That would be great. The whole “Bayesian model averaging” thing has been a useless time sink for decades.

I wouldn’t go that far — I agree that it’s a mess in the social sciences, especially
when models are hierarchical, but it has its use in applications with strong
generative/causal models. Although even in those applications where BMA
is useful the predictive performance should still be good!

The way I like to think about it is that _conditioned on the assumptions_ BMA
is the mathematically self-consistent and conceptually elegant right answer.
But in practice it’s almost impossible to build an analysis with the assumptions
held constant and BMA becomes highly vulnerable to overfitting and we need
some predictive diagnostics.

Anyways, there’s a nice write up of Mike West’s approach in http://arxiv.org/abs/1601.07463.
Basically if you assume that your individual models are exchangeable then
de Finetti’s theorem holds and any model weighting has to take the form of
a hierarchical model*. The individual weights could be marginal likelihoods,
in which case you recover BMA, or they could be sensitive to predictive
structure and the whole process can be iterated with new measurements.
Andrew’s two-stage strategy would be an approximation of this. In the end
you get this lovely self-consistent theoretical approach that unifies BMA and
predictive weighting.

* Caveat: Mike West didn’t agree with the appeal to exchangeability here.
But Dan Simpson agreed with me — in the end I believe de Finetti is an
if and only if, so the mixture form that is used has to imply exchangeability
of models.

[Ben Lambert]

Thanks all for your replies. So I am now attempting to code up Andrew's suggestion however, I had a query. I have K=4000 posterior predictions for each data point y, for each model.  Was the thinking that I would essentially average over all predictions for each data point, thus taking the mean, or median for example?

If so, was the model below what you were thinking? (Where yPred[i,j] is the average prediction for point i, in model j). Also, would you run this in RStan using: fit <- optimising(belowModel, data)? I have not used this functionality of Stan before.

data {

  int<lower=0> N; ## Total number of observations

  int<lower=0> K; ## Number of posterior predictive samples from model

  int yTrue[N];

  matrix[N,7] yPred;

  } 

parameters {

  simplex lambda[7]; ## Weightings parameters on each model

} 

transformed parameters {

  real<lower=0> squared_error;

  squared_error <- dot_self(yTrue - yPred * lambda);

} 

model {

  increment_log_prob(-squared_error);

}

generated quantities {

  real<lower=0> sigma_squared;

  sigma_squared <- squared_error / N;

}

[Dustin Tran]

This is cool. The language they use is a bit foreign, and I would never have guessed the problem they address based on their title. I’ll dig deeper. Thanks Michael

Dustin

[Aki Vehtari]

First summary

1) variations of stacking make separate models and then weight the predictions

2) BMA is just plain Bayesian approach with model uncertainty

3) in usual mixture model, the mixture components affect each other

4) in Xian's mixture model there are additional restrictions and so far they've been able to get practical computation only for two component mixture

Note that the usual mixture model is more flexible. Consider that we

have data from t_4 and two candiate models N and t_20 (which is almost

same as N). 

1) separate models will be almost identical and stacking will not be

   much different from either one

2) not much difference to 1 in this case

3) in a mixture N and t_20 can have different scales and the mixture

   can fit long tailed t_4 much better than either single model

4) ..?

So obviously in 1 and 2 we need to have models which behave

differently to get anything better than a single model. 

Andrew:

> What we’re talking about _does_ appear to be stacking, and I have no

> idea why this is not a standard approach in Bayesian statistics.

I think like Mike:

> The way I like to think about it is that _conditioned on the

> assumptions_ BMA is the mathematically self-consistent and

> conceptually elegant right answer.

Naturally the problem is "conditioned on the assumptions", so that if

we make bad assumptions then something which is not self-consistent

and elegant may produce better results. But if we can find a way to

produce better results than BMA, why we don't use that same approach

for continuous parameters if it's so good? It seems that it's possible

to find examples where something else works better than BMA, but not

something which would better often enough. 

> And why is it not in chapter 7 of BDA3?

Because there hasn't been convincing theory nor convincing

experiments. Maybe we can have it in the next edition?

I've been aware of some versions of stacking. I don't remember this

paper by Clarke, but can be also because it's not very convincing.

Although I advocate cross-validation I haven't used cross-validation

weighting expect for a few experiments because I haven't been aware of

theoretical justification similar to BMA. Clarke's paper provides some

experimental evidence but not good theoretical justification (and in

later paper he switch mean to median). There are also some papers

showing that point estimates of the parameters using cross-validation

are much noisier than by optimising the marginal likelihood. Which

raise the question why I'm advocating cross-validation for model

selection since parameter estimation is one kind of model selection,

too. In ... we show why we recommend cross-validation for model

selection only if there is a small number of models. It's possible

that when doing model averaging the noisiness is not as big problem

when doing model averaging.

BMA can be sensitive to priors if models have a large difference in

the number of parameters, in which case "pseudo BMA" (refering to term

pseudo Bayes factor) might be more stable, although I've seen also

case where uniform weighting performed better than cross-validation

weighting.

BMA can also difficult to compute (although often it's possible to

switch continuous model, like instead of spike-and-slab you can use

hierarchical shrinkage prior), and then pseudo BMA might be quick and

accurate enough approximation.

Now I have to go to sleep, but I'll read Mike's paper which Mike

recommended...

[Andrew Gelman]

This is all very interesting. I want to work it out in a simple example and then maybe I can come back to you with more questions. In particular, I agree that cross-validation optimization methods (including stacking) will not be optimal—and, given that stacking has been around for awhile and is _not_ so popular, it probably has problems—but asymptotically it should give the right answer, no? Regarding Aki’s other point, yes, this is a point that I have raised too, i think in my paper with Shalizi: when we do Bayesian continuous-parameter models, we’re doing “Bayesian model averaging” and I have no problem with it (usually), so why is BMA so bad for some (discrete) models. I think we can understand this based on some continuous-parameter models where the interpretation of one parameter changes when another parameter changes. Anyway, we can discuss more, I think this is important, also I agree with Bob that it can make sense to think of the mixture model, even if only as a comparison point.
A

[Ben Lambert]

Hello all, thanks again for your comments. Andrew, I wanted to ask what exactly you were thinking, in regards to the case where there are multiple predictions per data point. In particular, I have K=4000 predictions per data point. Were you thinking I would just 'stack' my observations on top of one another, and minimise the sum of squares, or just take the mean/median for each point, and use that instead? 

If the latter, was the model below what you were thinking? (Where yPred[i,j] is the average prediction for point i, in model j). Also, would you run this in RStan using: fit <- optimising(belowModel, data)? I have not used this functionality of Stan before.

data {

  int<lower=0> N; ## Total number of observations

[Andrew Gelman]

Hi, when you say K=4000 predictions, do you mean that you have 4000 posterior simulation draws? I guess I’m thinking you just take the mean of the 4000 sims so you get 1 prediction per data point per model, and that’s what you’d use as yPred.

And then, yes, what you have below is what I was thinking of. I assume it would be very fast.

Also, I’d never heard of dot_self()—that’s cool.

A

[Ben Lambert]

Great - thanks for confirming Andrew. I will let you know as soon as I have the results. Best, Ben

[Ben Lambert]

Ok, so I've run the above simplex model using the mean posterior predictions across 4000 samples. 

The trouble is that, the predictions across the 7 models are very highly autocorrelated - upwards of 0.97. This degree of correlation in the predictors means that the regression model is very unstable (sensitive to inclusion/dropping certain variables, i.e. the predictions from each model). 

The sum of this is that I end up with weights that are overly-biased towards the model which produces only marginally better predictions than its nearest competitor. For example, the correlation of each of the 7 model's predictions with the outcome is approximately: 0.81, 0.849, 0.82, 0.82, 0.847, 0.83 0.81. The weights the model suggests giving to each model are: 0 1 0 0 0 0 0. Note here that model 2 essentially gets all the weight even though its nearest competitor has a correlation of only 0.002 less with the outcome. I know that in my case, since all the models are very similar in terms of prediction, that essentially the weights I should want to find should be close to 1/7. It's just that I would like to do so in a less hand-wavy, more explicit and replicable way.

Basically the issue here is one of high degrees of autocorrelation amongst regressors. A way of dealing with it might be to use a penalty for weights, although what strength to give that penalty is unclear. Perhaps I could partition the dataset into training and test sets, and choose the weight which minimises the MSE?

Does anyone else have any other ideas?

Best,

Ben

[Andrew Gelman]

Interesting! If they all give very similar predictions, does it really matter much what you do?
A

[Ben Lambert]

Hi Andrew,

Whilst each model does yield similar predictions, in terms of parameters of interest, they differ. In particular, I am looking at determinants of organisms' mean lifetimes, through mark-release-recapture studies. 

Different assumptions can be made about how mortality varies with time, and dependent on the particular assumption made, this results in different predictions on mean lifetime. However, the crucial problem we have is that most studies are too short, or have too-low sample sizes, for us to differentiate between the different models in terms of their LOO-CV/WAIC; within error bounds of these criteria, all models perform relatively similarly. 

As such, I want to take this epistemic uncertainty in model choice into account when making estimating mean lifetime, by using some sort of weighting, where the weights are determined by just how well each model predicts the data. Intuitively, since each of the models does a similarly-good job of predicting the data, we would like assign roughly equal weights to each model's estimates. It's just whether there is a more data-driven way to find these weights. 

One particularly crude, and non-theory-derived way to do this would be to use something like:

weight of model i = 1 - (LOO-CV_i /(sum of LOO-CV across all models)

This would mean model i would be assigned a high weight if it has a good predictive accuracy, manifested through a low value of LOO-CV. Note: I only bring this up as a point for discussion, and am not actually suggesting that this is a sensible thing to do. I am only thinking that the above rule would get me approximately what I intuition suggests might be sensible.

Best,

Ben

On Fri, Feb 12, 2016 at 3:41 AM, Andrew Gelman <gel...@stat.columbia.edu> wrote:

Interesting!  If they all give very similar predictions, does it really matter much what you do?
A

[Andrew Gelman]

Hi, yes, this is a great example and worth thinking about further.  Yet another approach that has been suggested is to weight each model proportional to exp(-LOO).  In the asymptotic limit, this will pick the single best-fit model.  To the extent that different models fit well to different data, a mixture could makes sense.  This is somethign worth discussing.

A

[Avraham Adler]

On Thursday, February 11, 2016 at 11:05:39 PM UTC-5, Andrew Gelman wrote:

Hi, yes, this is a great example and worth thinking about further.  Yet another approach that has been suggested is to weight each model proportional to exp(-LOO).  In the asymptotic limit, this will pick the single best-fit model.  To the extent that different models fit well to different data, a mixture could makes sense.  This is somethign worth discussing.

A

Is that based on Burnham and Anderson's Akaike weighting? They weight to -AICc/2 because AIC is measured on the deviance scale.

Avi
 

[andy...@cubist-asia.com]

Ben,

I have also noticed the kind of stability issue you pointed out where I tried it on something that's very similar to what you have.

Instead of using the optimizer to find the weights, for me, going down the route of matrix algebra min variance seems to provide much more stable results.

Essentially, if we have a covariance matrix Sigma, it had been shown that the linear weight for min variance combination is inv(Sigma) * 1 / (1' * inv(Sigma) * 1). This is popular in finance, but I have also come across it from belief aggregation literature, e.g. PDF:http://faculty.washington.edu/ezivot/econ424/portfolioTheoryMatrix.pdf

Stan code to do this:
  vector compute_min_var_weights (matrix Sigma, int J) {
    matrix[J, J]  Sigma_inv;
    vector[J]     vec_1;
    row_vector[J] vec_1_t;
   
    vec_1   <- rep_vector(1, J);
    vec_1_t <- rep_row_vector(1, J);
    Sigma_inv <- inverse(Sigma);

    return (Sigma_inv * vec_1)/(vec_1_t * Sigma_inv * vec_1);
  }

I have looked at the following approach:
- 1: (from Andrew): First fit the 7 models separately.  Use these to get point predictions for each data point.
- 2: fit the historical error (actual - forecast) as a multi-normal (0, Sigma)
- 3: compute weights with above function.

Vs optimization, for my problem at least, this way seems to be more stable weights (mean of weights not going out of whack across a different seeds).
 
There is something tricky in step 2 with the covar/correlation estimation. I think it had been discussed elsewhere on the list, if your data size is not big and you use a LJK prior, the posterior of your correlation matrix can "under-estimate" the level of positive correlation and thus affects the min var computation. I have tried a few ways to get around this issue, but haven't given it much thought. Is there a "best practice" on how to fit a MV whose correlation (prior) is highly correlated?

Andy

[Bob Carpenter]

> On Feb 16, 2016, at 3:12 AM, andy...@cubist-asia.com wrote:
>
> ...

> There is something tricky in step 2 with the covar/correlation estimation. I think it had been discussed elsewhere on the list, if your data size is not big and you use a LJK prior, the posterior of your correlation matrix can "under-estimate" the level of positive correlation and thus affects the min var computation.

It works just like any other prior. The larger the parameter (1 is uniform),
the more the prior concentrates around the unit correlation matrix and the more
it pulls the posterior to be a unit. You also have the issue of the posterior
mean not being the same as the posterior mode, which might make it look like
an underestimate compared to a point estimate.

You can get an intuitive idea of how this works by using simulated data
where you know the answer and seeing what the fits look like.

> I have tried a few ways to get around this issue, but haven't given it much thought. Is there a "best practice" on how to fit a MV whose correlation (prior) is highly correlated?

Other than using the Cholesky factors for stability and some
regularization with LKJ, I don't think we have any other recommendations.
At least not written out.

- Bob