Exact model enumeration & Bayesian model averaging

Bayesian variable selection, Part 2 — 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 that $\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\le 20$-ish we can score all of them and get exact posterior model probabilities — no MCMC, no Monte-Carlo error.

This reproduces Congdon (2005), Example 3.3, on the famous cystic-fibrosis PEmax data (O'Neill et al. 1983; Altman 1991; ISwR::cystfibr): lung function in 25 patients against 9 clinical predictors. Congdon works through the models in forward stages; we do the stronger thing and enumerate all $2^9=512$. The from-scratch enumerator is checked two ways: against R's BAS (exact, same prior — must match to the decimal) and against a PyMC spike-and-slab search (shows the sampler of Part 1 recovers the same answer when enumeration would be infeasible).

The data (cystfibr_data.csv, $n=25$):

column	meaning
pemax	maximal static expiratory pressure (response) — a measure of respiratory muscle strength
age, sex	age (yr), sex (0 = male, 1 = female)
height, weight	height (cm), weight (kg)
bmp	body mass as % of normal for age/sex
fev1	forced expiratory volume in 1 s
rv, frc, tlc	residual volume, functional residual capacity, total lung capacity

The model and the closed-form marginal likelihood

For a model indexed by an inclusion vector $\gamma\in\{0,1\}^p$, keep the intercept and the columns $X_\gamma$ with $\gamma_j=1$: $$ y = \alpha\mathbf 1 + X_\gamma\beta_\gamma + \varepsilon, \qquad \varepsilon\sim\mathcal N(0,\sigma^2 I). $$

The conjugate choice is Zellner's $g$-prior on the (centered) coefficients, with a reference prior on $(\alpha,\sigma^2)$: $$ \beta_\gamma\mid\sigma^2,\gamma \sim \mathcal N\!\big(0,\; g\,\sigma^2 (X_\gamma^\top X_\gamma)^{-1}\big),\qquad 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 ordinary $R^2_\gamma$ (Liang et al. 2008): $$ p(y\mid\gamma) \;\propto\; \frac{(1+g)^{(n-1-p_\gamma)/2}}{\big[\,1+g(1-R^2_\gamma)\,\big]^{(n-1)/2}}, $$ with the intercept-only model ($p_\gamma=0,\;R^2=0$) as the reference (log marginal likelihood $0$). The numerator penalises size ($p_\gamma$), the denominator rewards fit ($R^2_\gamma$) — an automatic, fully-Bayesian Occam balance.

Hyperparameter (defined explicitly): $g=n=25$, the unit-information prior. It puts the weight of one observation on the prior, makes model comparison behave like BIC, and is exactly the setting used by R's BAS (alpha = n) — so the two must agree. A uniform prior is placed over the $2^p$ models. From the model scores we read three things: posterior model probabilities $p(\gamma\mid y)\propto p(y\mid\gamma)$, marginal inclusion probabilities $P(\gamma_j=1\mid y)=\sum_{\gamma:\gamma_j=1}p(\gamma\mid y)$, and model-averaged coefficients $\bar\beta_j=\sum_\gamma p(\gamma\mid y)\,\mathbb E[\beta_j\mid\gamma]$.

The algorithm — enumerate, don't search

No Markov chain here; the computation is a deterministic sweep:

1. Center $y$ and $X$ (the intercept is then implicit and always included).
2. For each of the $2^p$ subsets $\gamma$: regress $y$ on $X_\gamma$ by least squares, read off $R^2_\gamma$, and evaluate the closed-form $\log p(y\mid\gamma)$ above.
3. Normalise $\exp\{\log p(y\mid\gamma)\}$ over all models → exact posterior model probabilities.
4. Sum those probabilities over the models containing predictor $j$ → its inclusion probability; sum $p(\gamma\mid y)\cdot\tfrac{g}{1+g}\hat\beta_{\gamma}$ → the model-averaged coefficients (the $\tfrac{g}{1+g}$ factor is the $g$-prior's posterior shrinkage).

Because it is exact, there is nothing to converge or diagnose — the only limit is that $2^p$ must be enumerable ($512$ here is trivial; beyond $p\approx25$ one falls back to the stochastic search of Part 1). The engine is bvsenum.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import bvsenum as B                                                 # the from-scratch g-prior enumerator

names = ["age","sex","height","weight","bmp","fev1","rv","frc","tlc"]
d = pd.read_csv("cystfibr_data.csv"); X = d[names].values; y = d["pemax"].values

res = B.enumerate_bvs(y, X, g=len(y), names=names)                  # g = n = 25; enumerate all 2^9 = 512 models
print("exact enumeration of 2^9 = 512 models  (g-prior, g = n = 25)\n")
print("marginal inclusion probabilities:")
for nm, pi in sorted(zip(names, res["incl"]), key=lambda t: -t[1]): print("  %-7s  %.3f" % (nm, pi))
print("\ntop models by posterior probability:")
for inc, pr in B.top_models(res, 5): print("  %.3f  %s" % (pr, " + ".join(inc)))
size = np.arange(res["p"]+1); psize = np.array([res["w"][res["size"]==k].sum() for k in size])
print("\nE[model size] = %.1f" % (size*psize).sum())

# --- Fig: model uncertainty -- the top models, and the posterior over model size ---   -> bvs_2.png
tm     = B.top_models(res, 8)
tlabel = [" + ".join(inc) for inc, pr in tm]
tprob  = np.array([pr for inc, pr in tm])
ypos   = np.arange(len(tm))[::-1]
Emsize = (size * psize).sum()
fig, ax = plt.subplots(1, 2, figsize=(12.6, 4.4))
ax[0].barh(ypos, tprob, color="firebrick")
ax[0].set_yticks(ypos); ax[0].set_yticklabels(tlabel, fontsize=9)
ax[0].set_xlabel("posterior model probability")
ax[0].set_title("top models (no single winner)")
ax[0].grid(alpha=.25, axis="x")
ax[1].bar(size, psize, color="firebrick")
ax[1].axvline(Emsize, color="navy", lw=1.2, ls="--", label="E[size] = %.1f" % Emsize)
ax[1].set_xlabel("model size"); ax[1].set_ylabel("posterior probability")
ax[1].set_title("posterior over model size")
ax[1].set_xticks(size); ax[1].legend(); ax[1].grid(alpha=.25, axis="y")
plt.tight_layout(); plt.show()

exact enumeration of 2^9 = 512 models  (g-prior, g = n = 25)

marginal inclusion probabilities:
  weight   0.641
  fev1     0.596
  bmp      0.492
  age      0.376
  rv       0.303
  height   0.296
  tlc      0.235
  frc      0.225
  sex      0.213

top models by posterior probability:
  0.051  weight + bmp + fev1
  0.035  weight
  0.033  weight + bmp + fev1 + rv
  0.029  weight + bmp
  0.022  weight + bmp + fev1 + frc

E[model size] = 3.4

Reading the enumeration

• Three predictors carry the signal. Weight (0.64), FEV1 (0.60) and BMP (0.49) have the highest inclusion probabilities; the remaining six sit at or below the neutral 0.5 line. Clinically sensible — respiratory muscle strength tracks body size/nutrition (weight, BMP) and lung function (FEV1).
• But no single model dominates (left panel). The best model, weight + bmp + fev1, holds only ~5% of the posterior; dozens of models are plausible. With 25 patients and 9 correlated predictors, which model is genuinely uncertain — a point a single stepwise "winner" would hide entirely.
• The data prefers small models (right panel). The posterior over model size peaks at 2-3 predictors, averaging ~3.4 — the $g$-prior's built-in Occam penalty resisting the temptation 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 everything on one fit chosen under uncertainty.

import pymc as pm                                                    # cross-check: the Part-1 spike-and-slab SEARCH on the same data
Xz = (X - X.mean(0)) / X.std(0); yz = (y - y.mean()) / y.std(); p = 9
with pm.Model() as m:
    tau, c = 0.05, 20.0
    g    = pm.Bernoulli("g", 0.5, shape=p)
    b    = pm.Normal("b", 0, 1, shape=p)
    beta = pm.Deterministic("beta", b * pm.math.switch(g, c*tau, tau))
    sig  = pm.HalfNormal("sig", 1.0)
    pm.Normal("y", pm.math.dot(Xz, beta), sig, observed=yz)
    idata = pm.sample(4000, tune=2000, chains=4,
                      step=[pm.BinaryGibbsMetropolis([g]), pm.NUTS([b, sig], target_accept=0.95)],
                      random_seed=1, progressbar=False)
ssvs_ip = idata.posterior["g"].mean(("chain","draw")).values
print("pymc SSVS inclusion: ", "  ".join("%s=%.3f" % (n, pi)
      for n, pi in sorted(zip(names, ssvs_ip), key=lambda t: -t[1])))
print("-> stochastic search and exact enumeration agree on the leading predictors (weight and FEV1 clearly; BMP stays high, though the sampled run nudges age just above it).")

# --- Fig: inclusion probabilities, exact enumeration vs the spike-and-slab search ---   -> bvs_1.png
order = np.argsort(-res["incl"])
xlab  = [names[j] for j in order]
exact = res["incl"][order]
ssvs  = np.asarray(ssvs_ip)[order]
xp = np.arange(p); wbar = 0.4
fig, ax = plt.subplots(figsize=(8.4, 4.6))
ax.bar(xp - wbar/2, exact, wbar, color="firebrick", label="exact enumeration (g-prior)")
ax.bar(xp + wbar/2, ssvs,  wbar, color="navy",      label="spike-and-slab search (PyMC)")
ax.axhline(0.5, color="k", lw=.5, ls=":")
ax.set_xticks(xp); ax.set_xticklabels(xlab)
ax.set_ylabel("marginal inclusion probability")
ax.set_title("inclusion probabilities: exact enumeration vs spike-and-slab search")
ax.legend(); ax.grid(alpha=.25, axis="y")
plt.tight_layout(); plt.show()

Multiprocess sampling (4 chains in 4 jobs)

CompoundStep

>BinaryGibbsMetropolis: [g]

>NUTS: [b, sig]

Sampling 4 chains for 2_000 tune and 4_000 draw iterations (8_000 + 16_000 draws total) took 10 seconds.

pymc SSVS inclusion:  weight=0.651  fev1=0.549  age=0.442  bmp=0.417  height=0.387  rv=0.318  frc=0.272  tlc=0.229  sex=0.189
-> stochastic search and exact enumeration agree on the leading predictors (weight and FEV1 clearly; BMP stays high, though the sampled run nudges age just above it).

Results — exact evidence, and a search that finds it

• Enumeration and search agree on the substance. The 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 (here the sampled run nudges age just past BMP into third). 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 the prior and the algorithm.
• Why have both. Enumeration is exact and assumption-transparent but only feasible for small $p$; the Part-1 search scales to $p$ far beyond $2^p$-enumerable. Seeing them coincide here (where we can enumerate) is what licenses trusting the search where we cannot.
• The honest headline is uncertainty. This dataset is the textbook cautionary tale: classical stepwise selection reports a single confident model, but the full posterior shows the best model is only ~5% probable. Bayesian model averaging reports that uncertainty instead of hiding it.

The method's family

This is the conjugate / model-enumeration node (George-McCulloch 1997), complementing Part 1's stochastic-search node. The exact marginal likelihood is the same object BAS computes (next notebook), and the same quantity that, when $2^p$ is too large, is sampled rather than summed. Still ahead in the arc: VAR restriction search (George-Sun-Ni 2008), log-linear/contingency tables (Albert 1996), and unknown dimension via reversible-jump MCMC.

References

• George, E. I. & McCulloch, R. E. (1997). Approaches for Bayesian variable selection. Statistica Sinica 7, 339–373.
• Liang, F., Paulo, R., Molina, G., Clyde, M. A. & Berger, J. O. (2008). Mixtures of g-priors for Bayesian variable selection. JASA 103, 410–423.
• O'Neill, S. et al. (1983). The effects of chronic hyperinflation, nutritional status, and posture on respiratory muscle strength in cystic fibrosis. Am. Rev. Respir. Dis. 128, 1051–1054.
• Congdon, P. (2005). Bayesian Models for Categorical Data, Example 3.3. Wiley.

Data: cystfibr_data.csv — cystic-fibrosis PEmax data, 25 patients, 9 predictors (ISwR::cystfibr; Congdon Example 3.3). Next: bvs_R.ipynb — the R cross-check with BAS (exact enumeration, same g-prior).