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 $\sigma^2$ integrate out analytically, leaving a closed-form marginal likelihood $p(y\mid\text{model})$ . With $p$ predictors there are $2^p$ models; for $p\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 $2^9=512$ models rather than working through them in forward stages.

p(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 $\gamma\in\{0,1\}^p$ , the conjugate choice is Zellner's $g$ -prior on the coefficients, $\beta_\gamma\mid\sigma^2 \sim \mathcal N(0, g\sigma^2(X_\gamma^\top X_\gamma)^{-1})$ , with a reference prior $p(\alpha,\sigma^2)\propto 1/\sigma^2$ . Integrating out $\alpha,\beta_\gamma,\sigma^2$ gives a marginal likelihood that depends on the data only through the model's $R^2_\gamma$ (Liang et al. 2008) — the numerator penalises size, the denominator rewards fit, an automatic fully-Bayesian Occam balance. The hyperparameter is $g=n=25$ , the unit-information prior that 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 $y$ and $X$ ; for each of the $2^p$ subsets regress by least squares, read off $R^2_\gamma$ , and evaluate the closed-form $\log p(y\mid\gamma)$ ; normalise over all models for exact posterior model probabilities; then sum over models containing predictor $j$ for its inclusion probability, and weight the $g$ -shrunk estimates for model-averaged coefficients. Because it is exact, there is nothing to converge or diagnose — the only limit is that $2^p$ must be enumerable (beyond $p\approx 25$ one falls back to Part 1's search). The engine is a self-contained bvsenum.py.

Results

Three predictors carry the signal — weight (0.64), FEV1 (0.60) and BMP (0.49) have the highest inclusion probabilities, the other six at or below the neutral 0.5 line — clinically sensible, since respiratory muscle strength tracks body size/nutrition and lung function. 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 $g$ -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 leading predictors up front (weight and FEV1 clearly, with BMP still among the leaders); the magnitudes and exact ordering differ a little (inclusion probabilities are genuinely prior-dependent), but the set of predictors that matters is robust to both prior and algorithm. 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

George, E. I. & McCulloch, R. E. (1997). Approaches for Bayesian variable selection. Statistica Sinica 7(2), 339–373. — the conjugate model-enumeration framework reproduced here
Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g-prior distributions. In Bayesian Inference and Decision Techniques, 233–243. — the g-prior on the coefficients
Liang, F., Paulo, R., Molina, G., Clyde, M. A. & Berger, J. O. (2008). Mixtures of g priors for Bayesian variable selection. Journal of the American Statistical Association 103(481), 410–423. — the closed-form marginal likelihood in terms of R²_γ
Kass, R. E. & Wasserman, L. (1995). A reference Bayesian test for nested hypotheses and its relationship to the Schwarz criterion. Journal of the American Statistical Association 90(431), 928–934. — the unit-information prior (g = n) and the BIC connection
Raftery, A. E., Madigan, D. & Hoeting, J. A. (1997). Bayesian model averaging for linear regression models. Journal of the American Statistical Association 92(437), 179–191. — Bayesian model averaging
Hoeting, J. A., Madigan, D., Raftery, A. E. & Volinsky, C. T. (1999). Bayesian model averaging: a tutorial. Statistical Science 14(4), 382–401. — the BMA tutorial
Clyde, M. A., Ghosh, J. & Littman, M. L. (2011). Bayesian adaptive sampling for variable selection and model averaging. Journal of Computational and Graphical Statistics 20(1), 80–101. — the R BAS package used in the cross-check