Exact Model Enumeration & Bayesian Model Averaging

Python · R · Part 2 of the variable-selection arc — the conjugate route

Overview

Part 2 of the variable-selection arc: the conjugate route. Part 1 searched the model space with a spike-and-slab Gibbs sampler. But if the candidate set is small, we don't have to search at all — we can enumerate every model exactly. George & McCulloch (1997) showed how: put a conjugate prior on each model's coefficients so β\beta and σ2\sigma^2 integrate out analytically, leaving a closed-form marginal likelihood p(y∣model)p(y\mid\text{model}). With pp predictors there are 2p2^p models; for p≲20p\lesssim 20 we score all of them and get exact posterior model probabilities — no MCMC, no Monte-Carlo error. This reproduces Congdon's (2005) Example 3.3 on the cystic-fibrosis PEmax data (respiratory muscle strength in 25 patients against 9 clinical predictors), enumerating all 29=5122^9=512 models rather than working through them in forward stages.

p(y∣γ)  ∝  (1+g)(n−1−pγ)/2[ 1+g(1−Rγ2) ](n−1)/2p(y\mid\gamma) \;\propto\; \frac{(1+g)^{(n-1-p_\gamma)/2}}{\big[\,1+g(1-R^2_\gamma)\,\big]^{(n-1)/2}}

The g-prior and the closed-form marginal likelihood

For a model indexed by an inclusion vector γ∈{0,1}p\gamma\in\{0,1\}^p, the conjugate choice is Zellner's gg-prior on the coefficients, βγ∣σ2∼N(0,gσ2(Xγ⊤Xγ)−1)\beta_\gamma\mid\sigma^2 \sim \mathcal N(0, g\sigma^2(X_\gamma^\top X_\gamma)^{-1}), with a reference prior p(α,σ2)∝1/σ2p(\alpha,\sigma^2)\propto 1/\sigma^2. Integrating out α,βγ,σ2\alpha,\beta_\gamma,\sigma^2 gives a marginal likelihood that depends on the data only through the model's Rγ2R^2_\gamma (Liang et al. 2008) — the numerator penalises size, the denominator rewards fit, an automatic fully-Bayesian Occam balance. The single hyperparameter gg sets how wide the slab is — how large an effect the prior is prepared to entertain — and because the penalty for size grows with it, gg is precisely the dial controlling how hard the Occam balance bites. Here g=n=25g=n=25, the unit-information prior: the prior is given the weight of a single observation, the least informative choice that still lets the marginal likelihood be computed, which is what makes model comparison behave like BIC and matches R's BAS (alpha = n) exactly. The algorithm is a deterministic sweep, not a Markov chain: center yy and XX; for each of the 2p2^p subsets regress by least squares, read off Rγ2R^2_\gamma, and evaluate the closed-form log⁡p(y∣γ)\log p(y\mid\gamma); normalise over all models for exact posterior model probabilities; then sum over models containing predictor jj for its inclusion probability, and weight the gg-shrunk estimates for model-averaged coefficients. Because it is exact, there is nothing to converge or diagnose — the only limit is that 2p2^p must be enumerable (beyond p≈25p\approx 25 one falls back to Part 1's search). The engine is a self-contained bvsenum.py.

Results

Two predictors clear the line, and a third sits on it — weight (0.6410.641) and FEV1 (0.5960.596) are the only inclusion probabilities above the neutral 0.5 that a predictor starts with, with BMP third at 0.4920.492, fractionally below it, and the remaining six trailing from 0.376 down to 0.213. Clinically that is sensible, since respiratory muscle strength tracks body size, nutrition and lung function — but it is worth reading the numbers exactly rather than rounding them into a top three. With 25 patients and nine correlated predictors, only two predictors have had their prior odds moved upward at all. But no single model dominates: the best model, weight + bmp + fev1, holds only ~5% of the posterior; dozens are plausible, and with 25 patients and 9 correlated predictors which model is genuinely uncertain — a point a stepwise "winner" hides entirely. The posterior over model size peaks at 3–4 predictors (averaging ~3.4) — the gg-prior's built-in Occam penalty resisting the urge to use all nine. This is why we average: inclusion probabilities and model-averaged coefficients summarise all 512 models weighted by evidence rather than betting on one fit chosen under uncertainty. Enumeration and search agree on the substance: the Part-1 spike-and-slab MCMC — a different prior reached by sampling — puts the same two predictors up front — weight (0.651) and FEV1 (0.549) — while the ordering behind them shifts: the sampled run nudges age (0.442) just above BMP (0.417), reversing third and fourth place. The leading pair is robust to both prior and algorithm; the tail ordering is not, which is the honest reading of inclusion probabilities that are genuinely prior-dependent and separated by less than the spacing between them. Seeing them coincide here, where we can enumerate, is what licenses trusting the search where we cannot. The exact marginal likelihood is the same object R's BAS computes.

Notebooks

Python Exact model enumeration from scratch — the g-prior marginal likelihood over all 2⁹ = 512 models, posterior model probabilities, inclusion probabilities and Bayesian model averaging, with a PyMC spike-and-slab cross-check (Python) Open → R R cross-check — BAS exact g-prior enumeration, matching the posterior model and inclusion probabilities to the decimal (R) Open →

Downloads

bvsenum.py The from-scratch g-prior model enumerator — closed-form marginal likelihoods, posterior model probabilities, inclusion probabilities, and model-averaged coefficients (new self-contained module) cystfibr_data.csv Cystic-fibrosis PEmax data (O'Neill et al. 1983; ISwR::cystfibr; Congdon Example 3.3) — respiratory muscle strength vs 9 clinical predictors, n = 25

References