Discrete model comparision: Bayes factor in generated parameters?

[Kevin Miller]

Hi Stan Users,

I'd like to do some comparison of discrete models based on their model evidence p(M|y), and I'm wondering what the best way is to do this in stan. It seems to me that a good approach is to use Bayes factors, for which I need p(y|M) for each model. I'm imagining using stan to approximate the integral p(y|theta)*p(theta)*dTheta by sampling p(theta) in the Model block, computing the likelihood at each sample in the generated quantities block, and taking the average of these outside of stan. Code illustrating this idea for a simple hierarchical gaussian model is below the fold.

Alternatives on my mind include BIC, which I understand to constitute a minor statistical sin, and a technique called "thermodynamic integration", which I don't totally understand. I guess my questions are as much about statistics as they are about stan itself:

1) Does the approach below make sense as a way of getting estimates of p(y|M) for use in a Bayes Factor?

2) More generally: What are best practices in stan for comparing models based on model evidence, rather than on predictive performance?

Thanks!

Kevin

// Stan code for generating samples from p(theta), computing at each theta p(y|theta)

data {
     int num_subjects;
     int num_observations;
     real y[num_subjects, num_observations];

     real prior_mean;
     real prior_var;
     real prior_location;
     real prior_scale;
}

parameters {

     real pop_mean_mean;
     real<lower=0> pop_mean_var;
     real<lower=0> shared_var;
     real subj_mean[num_subjects];

}

model {

     pop_mean_mean ~ normal(prior_mean,prior_var);
     pop_mean_var ~ cauchy(prior_location, prior_scale);

     shared_var ~ cauchy(prior_location, prior_scale);
     subj_mean ~ normal(pop_mean_mean, pop_mean_var);

}

generated quantities {

     real log_lik;
     log_lik <- 0;

     for (subj_i in 1:num_subjects) {
          for (obs_i in 1:num_observations){
               log_lik <- log_lik + normal_log(y[subj_i,obs_i],subj_mean[subj_i],shared_var);
          }
     }
}

// Then back in matlab:

log_p_D_given_M = log(mean(exp(model_output.extract.log_lik)))

[Michael Betancourt]

I'm imagining using stan to approximate the integral p(y|theta)*p(theta)*dTheta by sampling p(theta) in the Model block, computing the likelihood at each sample in the generated quantities block, and taking the average of these outside of stan. Code illustrating this idea for a simple hierarchical gaussian model is below the fold.

Do not do this!  This estimator has infinite Monte Carlo variance

and is essentially worthless in practice.  Computing the evidence

is really hard and requires methods that transcend Markov chain

Monte Carlo.

Alternatives on my mind include BIC, which I understand to constitute a minor statistical sin, and a technique called "thermodynamic integration", which I don't totally understand. I guess my questions are as much about statistics as they are about stan itself:

1) Does the approach below make sense as a way of getting estimates of p(y|M) for use in a Bayes Factor?

No; see above.  BIC is also of almost no practical utility, and

thermodynamic integration is incredibly hard to tune and ensure

accurate results.  

Evidence calculations are in the next generation of problems

for Stan, and robust algorithms are currently being developed

theoretically.

2) More generally: What are best practices in stan for comparing models based on model evidence, rather than on predictive performance?

I’m sure Aki and Andrew will have more to say, but in general

you have to be very, very careful with Bayesian Model Comparison.

The evidence is incredibly sensitive to the choice of prior, which

makes it easy to “game” if you’re not careful.  There are some

circumstances where BMC can be useful (say, identifying the

number of modes in a mixture) but in general you have to

remember that Bayes is all about reporting uncertainty and

not making any hard choices as BMC does.  One strategy

that is becoming popular is to mix the models — i.e. don’t

try to select between models but rather make one big model

that encompasses both.  It’s easier algorithmically and much

more robust statistically (no coincidence that those two go

together).

[Andrew Gelman]

Hi, for reasons described in BDA (chapter 7 of 3rd edition) and in my 1995 paper with Rubin, I don’t generally recommend Bayes factors or BIC.  I think this stuff is a many-decades-long dead end.

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.
For more options, visit https://groups.google.com/d/optout.

[Bob Carpenter]

So as not to be entirely negative, here's the latest from
Andrew, Aki, and Jonah on model comparison:

http://andrewgelman.com/2015/07/16/new-papers-on-loowaic-and-stan/

I believe their thinking is still evolving given the chain of
papers over the last few years. It's a much trickier topic than
short intros to Bayes lead people to believe.

- Bob

[Kevin Miller]

Thanks, everyone, for the quick and helpful replies!  I'm seeing that model comparison is indeed a trickier topic than I believed!  I think that what I need might be advice at a higher level of abstraction. Let me explain what I'm trying to do:

The common situation in my field (psychology/neuroscience) is to have several models which instantiate hypotheses about the generative process that gave rise to a dataset, and our usual goal for model comparison is to determine the relative plausibility of the hypotheses.  It seems to many of us (though I'd love to find out that we're wrong about this) that predictive-performance-based methods like cross validation and WAIC are not ideal for this purpose, since they often favor models that are too complex to be plausible hypotheses about how the brain works, especially when datasets are large and overfitting is minimal. The typical solution is to use BIC, but I'd like to do better, especially since I'd like to use hierarchical models, where BIC seems particularly terrible.

Do you have any recommendations for hypothesis testing in situations like these? Perhaps build a mixture model (as Michael suggests), and ask whether the posterior over the mixing weight is mostly to one side or the other of 0.5? Or perhaps our reasons for avoiding cross validation, etc., are wrongheaded, and there is a way to make these methods work?
I guess it's also possible that there really is no good general quantitative solution, and we should do something like report a measure of predictive performance and report model structure, and do our "plausibility" evaluations qualitatively.

I realize I'm now outside the realm of stan and asking for general statistics help - any advice will be very much appreciated!

Thanks a ton!

Kevin

[Aki Vehtari]

On Wednesday, August 5, 2015 at 9:41:15 AM UTC+3, Kevin Miller wrote:

The common situation in my field (psychology/neuroscience) is to have several models which instantiate hypotheses about the generative process that gave rise to a dataset, and our usual goal for model comparison is to determine the relative plausibility of the hypotheses.  It seems to many of us (though I'd love to find out that we're wrong about this) that predictive-performance-based methods like cross validation and WAIC are not ideal for this purpose, since they often favor models that are too complex to be plausible hypotheses about how the brain works, especially when datasets are large and overfitting is minimal. The typical solution is to use BIC, but I'd like to do better, especially since I'd like to use hierarchical models, where BIC seems particularly terrible.

In psychology and neuroscience the predictive approach is very natural as you should be able to predict the behavior of the not yet observed person or brain. If you have several observations per subject, you could use leave-one-subject-out cross-validation (which is commonly used by neuroscientist I know). If you a large number of subjects, you could use k-fold-CV such that data division is made by subjects.

Aki

[Kevin Miller]

Hi Aki,

Thanks for your response! I'm enjoying your papers with Andrew on WAIC and on PSIS-LOO, and I'm excited about using these tools in some of my work.

I'm not sure, though, that predictive performance measures are ideal for all of our applications. Often in neuroscience we consider out models not as tools to make predictions about the future (as one might in e.g. weather prediction or stock picking), but rather as competing hypotheses about the generative process that gave rise to the dataset (e.g. the neural algorithm a rat is using to solve a particular task). In this case, it seems to me that the Bayes factor, with it's "automatic Occam's razor", is the natural model comparison metric.

To take an extreme example, a model which is just a lookup table giving the probability of each outcome occurring in each condition will have excellent predictive performance if constrained by enough data, but this should not make us expect to find a lookup table in the brain. A Bayes factor analysis would correctly give this model a very low score, due to the huge number of free parameters.  I guess I've learned here that Bayes factors are problematic to reliably estimate, so I'm not sure where that leaves us if we want to compare models of this type.

Thanks!

Kevin

[Andrew Gelman]

Hi, Kevin.  I agree that predictive performance is not everything; this is a point we discuss here:

http://www.stat.columbia.edu/~gelman/research/published/final_sub.pdf

But marginal likelihood is not in general a good solution because it depends on irrelevant aspects of the model (as mentioned in an earlier email, change a prior from N(0,100) to N(0,1000) and you change the marginal likelihood by a factor of 10).  BIC is what it is, it can be useful, but don’t kid yourself that it is Bayesian in any real sense!

A

[Sebastian Weber]

Hi!

To my knowledge you can calculate Bayes factors using INLA. Of course, if you can write your model in INLA and this also doesn't mean it is useful.

... but I have never done this myself as INLA was always too restrictive for my models.

Best,

Sebastian

BTW, very interesting discussion!

[Andrew Gelman]

As long as you have proper priors, you can calculate Bayes factors.  One way to do so is using path sampling, which we could probably implement as an R function calling Stan.  Maybe this is worth doing!  But I think that in most applications it would be a bad idea.

A

[Aki Vehtari]

Often in neuroscience we consider out models not as tools to make predictions about the future (as one might in e.g. weather prediction or stock picking), but rather as competing hypotheses about the generative process that gave rise to the dataset (e.g. the neural algorithm a rat is using to solve a particular task).

I assume that you would like that your model could predict how the other rats will perform in task solving? For me the first step is to check that model generalizes out of the observed data. If there are several models which give similar predictions then it's possible to continue examining which models are plausible given other information known or measurable in the future.

I have worked with neuroscientist, and in several cases Bayes factor has been supersensitive to priors on parameters and thus completely useless for model selection (see also Andrew's comment). So it's not just the difficulty of estimating the Bayes factors, it's also that Bayes factor can be super sensitive to priors.

Aki