Bayesian Model Comparison: Bayes Factors, PSIS-LOO and WAIC

Two regression models differ only in the prior on one coefficient: normal with standard deviation 1 in one, 10 in the other. The estimate sits well inside both priors and the data are informative, so the two posteriors are almost indistinguishable and predictions agree to the second decimal. Yet the Bayes factor favours the tighter-prior model by roughly a factor of ten. Both facts are correct, because the two procedures answer different questions. One asks which prior predicted the data better before seeing it. The other asks which fitted model will predict new data better.

Marginal likelihood and Bayes factors

The marginal likelihood of model \(M\) integrates the likelihood over the prior:

and the Bayes factor is the ratio \(\text{BF}_{12} = p(y \mid M_1) / p(y \mid M_2)\). Multiplied by prior model odds it gives posterior model odds, which makes it the coherent answer when one of the candidate models is true (Kass and Raftery, 1995, Bayes Factors, JASA 90(430)).

The sensitivity in the opening example follows from a Laplace-style approximation. When the likelihood is concentrated relative to the prior, \(p(y \mid M) \approx p(y \mid \hat\theta)\, p(\hat\theta \mid M)\, w_{\text{post}}\), where \(w_{\text{post}}\) is the posterior width. Widening a roughly flat prior tenfold divides the prior density at \(\hat\theta\) by about ten and leaves the other terms nearly unchanged. The marginal likelihood falls tenfold while the posterior barely moves. Taken to the limit, an improper flat prior makes the Bayes factor undefined, and a very diffuse prior on an extra parameter makes the simpler model win regardless of the data, the Jeffreys-Lindley paradox.

Computation is its own problem. The marginal likelihood is an integral over the prior, where most mass sits in regions of negligible likelihood, and the harmonic mean of likelihoods over posterior draws, the tempting shortcut, has notoriously unstable variance. Bridge sampling and nested sampling work, at substantially higher cost than fitting the model.

Expected log predictive density

The predictive alternative targets out-of-sample fit. For \(n\) observations, the expected log pointwise predictive density is

where \(\tilde y_i\) is a new draw from the true data-generating process. Leave-one-out cross-validation estimates it with \(\text{elpd}_{\text{loo}} = \sum_i \log p(y_i \mid y_{-i})\). Refitting the model \(n\) times is usually infeasible, but with \(S\) posterior draws \(\theta^{s}\) each leave-one-out predictive can be estimated by importance sampling, using ratios \(r_i^s \propto 1/p(y_i \mid \theta^s)\) to reweight the full posterior toward the posterior without \(y_i\).

Raw importance ratios have heavy tails exactly when observation \(i\) is influential, and the estimate becomes dominated by a few draws. Pareto-smoothed importance sampling fits a generalised Pareto distribution to the largest ratios and replaces them with smoothed quantiles (Vehtari, Gelman and Gabry, 2017, Practical Bayesian model evaluation using leave-one-out cross-validation and WAIC, Statistics and Computing 27). The fitted shape parameter \(\hat{k}\) doubles as a diagnostic. The current guidance treats the estimate as reliable when \(\hat{k} < \min(1 - 1/\log_{10} S,\ 0.7)\) (Vehtari et al., 2024, Pareto Smoothed Importance Sampling, JMLR 25). With \(S = 1{,}000\) draws that threshold is 0.67; with \(S = 4{,}000\) it is 0.7. Observations above it should be refitted exactly or handled with \(K\)-fold cross-validation.

WAIC, and why LOO is usually preferred

WAIC estimates the same quantity analytically from pointwise posterior log-likelihoods:

The first term is the in-sample log predictive density; the second, the effective number of parameters, corrects its optimism. Watanabe proved WAIC is asymptotically equivalent to Bayesian leave-one-out cross-validation even for singular models such as mixtures and neural networks, where AIC's parameter count is meaningless (Watanabe, 2010, JMLR 11).

In finite samples with weak priors or influential observations WAIC can be noticeably biased, and it offers no signal when that is happening. PSIS-LOO has the \(\hat k\) diagnostic, which is the main reason Vehtari, Gelman and Gabry recommend it over WAIC.

When it breaks

The two camps answer different questions, and they disagree about which matters. Bayes factors assume the true model is among the candidates (the \(\mathcal{M}\)-closed view). Predictive criteria assume none is, and ask which approximation predicts best. For combining models the disagreement is concrete: Bayesian model averaging weights by posterior model probability and, as data grow, puts all weight on one model, while stacking optimises weights for predictive performance and keeps a mixture when models are complementary (Yao, Vehtari, Simpson and Gelman, 2018, Using Stacking to Average Bayesian Predictive Distributions, Bayesian Analysis 13(3)).

Differences need standard errors. The elpd difference between models has a standard error computed from the pointwise differences, roughly \(\sqrt{n \operatorname{Var}_i(\text{elpd}_i^{A} - \text{elpd}_i^{B})}\). A difference of 3 with a standard error of 4 is not a preference.

Leave-one-out is the wrong unit for structured data. For time series, leaving out one point lets the model interpolate from its neighbours; the relevant quantity is leave-future-out. For hierarchical models predicting new groups, leave-one-group-out is the target, and PSIS approximations for it tend to fail more often.

Models must predict the same thing. Comparing a model of \(y\) against a model of \(\log y\) without the Jacobian adjustment compares densities on different scales and produces meaningless numbers.