Stochastic Search Variable Selection (SSVS)

Python · R · Part 1 of the variable-selection arc — spike-and-slab

Overview

Part 1 of the Bayesian variable-selection arc: spike-and-slab, from scratch. Given a response and a pile of candidate predictors, which belong in the model? The Bayesian answer is not a single "best" model but a posterior distribution over the $2^p$ possible models. George & McCulloch (1993) made this computable with a trick now central to Bayesian model selection: attach to each coefficient a latent inclusion indicator $\gamma_j\in\{0,1\}$ and a spike-and-slab prior, then let a Gibbs sampler wander through model space — the stochastic search. This notebook builds SSVS by hand and reproduces Congdon's (2005) selection example — the classic Draper–Smith turnip-greens experiment (vitamin-B2 / riboflavin content against three environmental drivers, $n=27$ plots) — where the search should discover that essentially one predictor, soil-moisture tension, carries the signal.

\beta_j\mid\gamma_j=0 \sim \mathcal N(0,\tau^2)\ \text{(spike)}, \quad \beta_j\mid\gamma_j=1 \sim \mathcal N(0,c^2\tau^2)\ \text{(slab)}, \quad \gamma_j\sim\text{Bernoulli}(w)

The spike-and-slab Gibbs sampler — the search

Each coefficient gets a spike-and-slab prior indexed by $\gamma_j$ : a spike $\mathcal N(0,\tau^2)$ pinning a switched-off coefficient near zero, and a slab $\mathcal N(0,c^2\tau^2)$ allowing a real effect, with $\gamma_j\sim\text{Bernoulli}(w)$ . The key to George–McCulloch is that the spike is a narrow but proper normal, not a point mass — so with two conjugate normals every full conditional is closed form and the sampler mixes far better than a hard point-mass prior. The three-block Gibbs sampler cycles: (1) the coefficients $\beta$ as one joint Gaussian draw (a single Cholesky of the posterior precision); (2) the indicators $\gamma_j$ , an independent Bernoulli from the ratio of the slab and spike densities at the current $\beta_j$ — this single line is the "search", letting the chain toggle predictors in and out; and (3) the error variance $\sigma^2$ from its Inverse-Gamma conditional. Because the indicators update jointly with continuously-varying coefficients, the chain visits many models and returns a full posterior over model space — the quantity of interest being the posterior inclusion probability $P(\gamma_j=1\mid y)$ . The engine lives in a self-contained ssvs.py; the same problem is fit a second way in PyMC (a gradient sampler for $\beta,\sigma$ married to a discrete sampler for $\gamma$ ).

Results

One predictor survives. Soil-moisture tension has a posterior inclusion probability of ~0.69, every other term sits near or below its prior value of 0.5, and the single most-visited model is "moisture only" (~18% of the posterior) — exactly Congdon's Model B, rediscovered by the search (the driest plots yield the least riboflavin). The spike-and-slab is visible in the posteriors: the moisture coefficient is broad and centred away from zero (it lives in the slab), while a dropped predictor like temperature is squeezed into a tight spike at zero — the sampler is not shrinking it a little, it is repeatedly setting its indicator to 0. The only near-rival is the moisture×temperature interaction (~0.56), which picks up correlated signal but rarely enters without moisture. The selection is a property of the model and data, not the algorithm: the from-scratch conjugate Gibbs and PyMC's NUTS-plus-binary-Gibbs give inclusion probabilities within a few points everywhere (moisture ≈ 0.69 in both) — and the exercise shows why George–McCulloch derived a fully-conjugate Gibbs scheme, since the discrete sampler is the natural fit for the sharp spike/slab geometry. What SSVS delivers over a single fit is not a point estimate but posterior inclusion probabilities and a ranked model space — a coherent statement of which predictors matter and how sure we are.

Notebooks

Python SSVS from scratch — the George–McCulloch spike-and-slab Gibbs sampler, plus a PyMC compound-sampler cross-check, on the Draper–Smith / Congdon riboflavin example (Python) Open → R R cross-check — BAS exact enumeration of all 2⁶ = 64 models, confirming the inclusion probabilities (R) Open →

Downloads

ssvs.py The from-scratch George–McCulloch spike-and-slab Gibbs engine — joint coefficient draw, indicator Bernoulli update, and variance conditional (new self-contained module) ssvs_data.csv The Draper–Smith (1966) turnip-greens vitamin-B2 experiment (Congdon Example 3.10): riboflavin vs radiation, soil-moisture tension, and temperature, n = 27

References

George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. Journal of the American Statistical Association 88(423), 881–889. — the SSVS spike-and-slab Gibbs sampler reproduced here
George, E. I. & McCulloch, R. E. (1997). Approaches for Bayesian variable selection. Statistica Sinica 7(2), 339–373. — spike-and-slab variants and the conjugate formulation
Mitchell, T. J. & Beauchamp, J. J. (1988). Bayesian variable selection in linear regression. Journal of the American Statistical Association 83(404), 1023–1032. — the original spike-and-slab prior
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 exact-enumeration cross-check
O'Hara, R. B. & Sillanpää, M. J. (2009). A review of Bayesian variable selection methods: what, how and which. Bayesian Analysis 4(1), 85–117. — a survey placing SSVS among the alternatives
Congdon, P. (2005). Bayesian Models for Categorical Data. Wiley. — the worked riboflavin selection example (Model B)
Draper, N. R. & Smith, H. (1966). Applied Regression Analysis. Wiley. — the turnip-greens vitamin-B2 experiment