Stochastic Search Variable Selection (SSVS)

Bayesian variable selection, Part 1 — spike-and-slab from scratch

Given a response and a pile of candidate predictors, which predictors 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 that is now the workhorse of 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 first notebook builds SSVS from scratch and reproduces the selection example of Congdon (2005), Bayesian Models for Categorical Data, Example 3.10, whose data is the classic turnip-greens experiment of Draper & Smith (1966), *Applied Regression Analysis***: the **vitamin B2 (riboflavin) content of turnip greens measured against three environmental drivers, over $n=27$ plots. The search should discover that essentially one predictor — soil-moisture tension — carries the signal. Two engines — a hand-derived conjugate Gibbs sampler and a PyMC compound sampler — plus an R cross-check (BAS, next notebook) must all agree.

The data (ssvs_data.csv; predictors are the book's pre-scaled values):

column	variable	definition
Y	riboflavin	vitamin B2 content of the turnip greens
X1	radiation	relative gram-calories per minute over the preceding half-day, ÷100
X2	soil-moisture tension	average soil-moisture tension, ÷1000 — set at three controlled levels (0.02, 0.07, 0.474), nine plots each
X3	temperature	Fahrenheit air temperature, ÷10 (values 6.2–8.9 ⇒ 62–89 °F)

From these three mains we also form the three pairwise interactions, giving $p=6$ candidate terms. Soil-moisture tension is a designed 3-level factor and, on its own, correlates $-0.86$ with riboflavin (the driest plots yield markedly less B2); radiation and temperature correlate only $+0.10$ and $+0.12$. So this is a clean test case: the search should recover moisture alone — good for trusting the machinery before harder problems.

The model

A linear regression with $p$ candidate predictors (here $p=6$: the mains radiation, moisture, temperature and their pairwise interactions rad:moist, rad:temp, moist:temp), all columns standardized so a single pair of prior scales governs every coefficient:

$$ y_i = \sum_{j=1}^p \beta_j x_{ij} + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal N(0,\sigma^2). $$

Each coefficient gets a spike-and-slab prior indexed by a latent 0/1 inclusion indicator $\gamma_j$:

$$ \beta_j \mid \gamma_j=0 \;\sim\; \mathcal N(0,\;\tau^2)\quad\text{(spike: pinned near zero)},\qquad \beta_j \mid \gamma_j=1 \;\sim\; \mathcal N(0,\;c^2\tau^2)\quad\text{(slab: a real effect is allowed)}, $$ $$ \gamma_j \sim \text{Bernoulli}(w), \qquad \sigma^2 \sim \text{Inverse-Gamma}(a_0,b_0). $$

The spike is a narrow but proper normal (not a point mass at 0), which is the key to George-McCulloch: with two conjugate normals every full conditional is closed form, and the sampler mixes far better than a hard point-mass prior. The quantity we care about is the posterior inclusion probability $P(\gamma_j=1\mid y)$, estimated by the fraction of Gibbs sweeps in which predictor $j$ is switched on.

Hyperparameters (all fixed explicitly):

symbol	value	meaning
$\tau$	$0.05$	spike standard deviation — a "switched-off" coefficient is confined to $[-0.1,0.1]$ (standardized units)
$c$	$20$	slab-to-spike ratio; slab sd $=c\tau = 1.0$, wide enough for any standardized effect
$w$	$0.5$	prior inclusion probability — each predictor is a coin flip a priori, so the prior is neutral over model size
$a_0,b_0$	$10^{-3},10^{-3}$	vague Inverse-Gamma prior on the error variance $\sigma^2$

The threshold at which the spike and slab densities cross, $\beta^\star = \tau\sqrt{2c^2\log c/(c^2-1)}\approx 0.12$, is the implicit "is this coefficient big enough to matter?" boundary the search uses.

The Gibbs sampler

SSVS cycles three closed-form full conditionals (with $y$ and $X$ standardized). Let $D^2=\text{diag}(d_j^2)$ hold the current prior variances, $d_j^2 = \tau^2$ if $\gamma_j=0$ and $c^2\tau^2$ if $\gamma_j=1$.

1. Coefficients $\;\beta \mid \gamma,\sigma^2,y \sim \mathcal N(\mu,\Sigma)$, a joint Gaussian with $$ \Sigma=\Big(\tfrac{X'X}{\sigma^2}+D^{-2}\Big)^{-1},\qquad \mu=\Sigma\,\tfrac{X'y}{\sigma^2}. $$ Drawn with one Cholesky factor of the posterior precision.

2. Indicators $\;\gamma_j \mid \beta_j$ — an independent Bernoulli for each $j$, from the ratio of the two prior densities evaluated at the current $\beta_j$: $$ P(\gamma_j=1\mid\beta_j)=\frac{w\,\mathcal N(\beta_j;0,c^2\tau^2)}{w\,\mathcal N(\beta_j;0,c^2\tau^2)+(1-w)\,\mathcal N(\beta_j;0,\tau^2)}. $$ A large $|\beta_j|$ is far more likely under the wide slab, so $\gamma_j\to1$; a $\beta_j$ near zero is better explained by the spike, so $\gamma_j\to0$. This single line is the "search" — it is what lets the chain toggle predictors in and out.

3. Error variance $\;\sigma^2 \mid \beta,y \sim \text{Inverse-Gamma}\big(a_0+\tfrac n2,\; b_0+\tfrac12\lVert y-X\beta\rVert^2\big).$

Because the indicators are updated jointly with continuously-varying coefficients, the chain visits many models and returns a full posterior over the model space, not just a point estimate. The engine lives in ssvs.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import ssvs                                                          # the from-scratch George-McCulloch sampler

d = pd.read_csv("ssvs_data.csv"); X = d[["X1","X2","X3"]].values; Y = d["Y"].values  # radiation, moisture, temp
names = ["radiation","moisture","temp","rad:moist","rad:temp","moist:temp"]
D = np.column_stack([X[:,0], X[:,1], X[:,2],                        # 3 main effects ...
                     X[:,0]*X[:,1], X[:,0]*X[:,2], X[:,1]*X[:,2]])  # ... + 3 pairwise interactions

res = ssvs.ssvs_gibbs(Y, D, ndraw=12000, burn=6000, tau=0.05, c=20.0, w=0.5, seed=1)
print("marginal inclusion probabilities  P(gamma_j = 1 | data):")
for nm, pi in zip(names, res["incl"]): print("  %-11s  %.3f" % (nm, pi))
print("\ntop models by posterior probability:")
for inc, pr in ssvs.top_models(res, names, 5): print("  %.3f  %s" % (pr, " + ".join(inc)))

# --- Fig: the visited model space, and the spike-and-slab mechanism ---
fig, (axL, axR) = plt.subplots(1, 2, figsize=(12.5, 4.3))

# left: posterior probability of the most-visited models (what the stochastic search explored)
tm     = ssvs.top_models(res, names, 8)
labels = [" + ".join(inc) for inc, pr in tm][::-1]
probs  = [pr for inc, pr in tm][::-1]
axL.barh(range(len(probs)), probs, color="firebrick")
axL.set_yticks(range(len(probs))); axL.set_yticklabels(labels, fontsize=9)
axL.set_xlabel("posterior model probability")
axL.set_title("visited model space (top 8 models)")
axL.grid(alpha=.25, axis="x")

# right: spike-and-slab in action -- a selected coef lives in the broad slab, a dropped one in the tight spike
axR.hist(res["beta"][:, 1], bins=60, density=True, color="firebrick", alpha=.6, label="moisture (selected -> slab)")
axR.hist(res["beta"][:, 2], bins=60, density=True, color="steelblue", alpha=.6, label="temp (dropped -> spike)")
axR.axvline(0, color="k", lw=.5, ls=":")
axR.set_xlabel("standardized coefficient"); axR.set_ylabel("posterior density")
axR.set_title("spike-and-slab mechanism")
axR.legend(fontsize=9); axR.grid(alpha=.25)

fig.tight_layout(); plt.show()

marginal inclusion probabilities  P(gamma_j = 1 | data):
  radiation    0.202
  moisture     0.694
  temp         0.163
  rad:moist    0.271
  rad:temp     0.236
  moist:temp   0.558

top models by posterior probability:
  0.176  moisture
  0.119  moisture + moist:temp
  0.105  moist:temp
  0.052  moisture + rad:moist
  0.039  moisture + rad:temp

Reading the from-scratch result

• 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). This is exactly Congdon's Model B — the search rediscovers that moisture alone carries the signal (the driest plots yield the least riboflavin).
• The spike-and-slab in action (right panel). The moisture coefficient's posterior is broad and centred away from zero (it lives in the slab); a dropped predictor like temperature is squeezed into a tight spike at zero. The sampler is not shrinking temperature a little — it is repeatedly setting its indicator to 0 and pinning it.
• moist:temp is the only near-rival (~0.56): because moisture matters so much, its interaction with temperature picks up correlated signal. But it rarely enters without moisture, and no model beats moisture-alone.

import pymc as pm                                                    # second engine: same model, gradient + discrete sampler
Dz = (D - D.mean(0)) / D.std(0); yz = (Y - Y.mean()) / Y.std()      # PyMC works on the standardized design
with pm.Model() as m:
    tau, c, w = 0.05, 20.0, 0.5
    g    = pm.Bernoulli("g", w, shape=6)                            # inclusion indicators (discrete)
    b    = pm.Normal("b", 0, 1, shape=6)
    beta = pm.Deterministic("beta", b * pm.math.switch(g, c*tau, tau))   # spike/slab scale set by g
    sig  = pm.HalfNormal("sig", 1.0)
    pm.Normal("y", pm.math.dot(Dz, beta), sig, observed=yz)
    idata = pm.sample(3000, tune=2000, chains=4,                    # NUTS for (b, sig), binary Gibbs for g
                      step=[pm.BinaryGibbsMetropolis([g]), pm.NUTS([b, sig], target_accept=0.95)],
                      random_seed=1, progressbar=False)
pmincl = idata.posterior["g"].mean(("chain","draw")).values
print("pymc inclusion: ", "  ".join("%s=%.3f" % (n, p) for n, p in zip(names, pmincl)))
print("-> from-scratch and PyMC agree: soil-moisture tension is the one predictor that survives; moist:temp a distant second.")

# --- Fig: posterior inclusion probabilities, from-scratch Gibbs vs PyMC ---
fig, ax = plt.subplots(figsize=(8.2, 4.2))
xpos = np.arange(len(names)); wdt = 0.4
ax.bar(xpos - wdt/2, res["incl"], wdt, color="firebrick", label="from-scratch Gibbs")
ax.bar(xpos + wdt/2, pmincl,      wdt, color="#4a4a4a",   label="PyMC")
ax.axhline(0.5, color="k", lw=.5, ls=":")
ax.set_xticks(xpos); ax.set_xticklabels(names, rotation=30, ha="right")
ax.set_ylabel("P(gamma_j = 1 | data)"); ax.set_ylim(0, 1)
ax.set_title("posterior inclusion probabilities: from-scratch Gibbs vs PyMC")
ax.legend(); ax.grid(alpha=.25, axis="y")
fig.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 3_000 draw iterations (8_000 + 12_000 draws total) took 9 seconds.

pymc inclusion:  radiation=0.191  moisture=0.689  temp=0.166  rad:moist=0.273  rad:temp=0.233  moist:temp=0.582
-> from-scratch and PyMC agree: soil-moisture tension is the one predictor that survives; moist:temp a distant second.

Results — the search agrees across two independent engines

• Same answer, two samplers. 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). The selection is a property of the model and the data, not the algorithm.
• Discrete Gibbs is the natural sampler for spike-and-slab. PyMC has to marry a gradient sampler (for $\beta,\sigma$) to a discrete one (for $\gamma$) and needs a high target_accept to tame the sharp spike/slab geometry — precisely why George-McCulloch derived a fully-conjugate Gibbs scheme in the first place. Both work here; the from-scratch sampler is cleaner and faster on this problem.
• What SSVS delivers over a single fit. 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. The next notebook checks the whole thing against BAS, which enumerates all $2^6=64$ models exactly.

The method's family (for the wider selection arc)

This is the spike-and-slab / SSVS node. Congdon 3.10 actually codes the Kuo-Mallick (1998) product form ($\beta_j\!\cdot\!\gamma_j$ with a genuine point mass); we use George-McCulloch's proper-normal spike, which mixes better and gives the same selection. Later notebooks in this arc will pick up the conjugate / model-enumeration branch (George-McCulloch 1997; Congdon 3.3, cystic-fibrosis), VAR restriction search (George-Sun-Ni 2008), and unknown dimension via reversible-jump MCMC (Congdon 3.8/3.9).

References

• George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. JASA 88, 881–889.
• George, E. I. & McCulloch, R. E. (1997). Approaches for Bayesian variable selection. Statistica Sinica 7, 339–373.
• Kuo, L. & Mallick, B. (1998). Variable selection for regression models. Sankhyā B 60, 65–81.
• Draper, N. R. & Smith, H. (1966). Applied Regression Analysis. Wiley — source of the turnip-greens riboflavin data.
• Congdon, P. (2005). Bayesian Models for Categorical Data, Example 3.10. Wiley.

Data: ssvs_data.csv — turnip-greens vitamin-B2 experiment (Draper & Smith 1966; Congdon Example 3.10): riboflavin vs radiation (X1), soil-moisture tension (X2), temperature (X3); n = 27. Next: ssvs_R.ipynb — the R cross-check with BAS (exact enumeration of all 64 models).