Variable Selection
Which predictors belong in the model? The Bayesian answer is not a single winner but a posterior over the whole model space — inclusion probabilities and a ranked list of models, with honest uncertainty. This arc builds that machinery from scratch: stochastic search variable selection (SSVS) with George–McCulloch spike-and-slab priors; exact model enumeration with conjugate g-priors and Bayesian model averaging; VAR restriction search that selects a whole system's lag structure and error covariance at once; log-linear model selection that chooses which associations belong in a contingency table; and, when even the model's size is unknown, a trilogy of trans-dimensional samplers — a reversible-jump change-point model, a Dirichlet-process mixture, and free-knot splines — inferring the number of regimes, of components, and of knots. Further nodes extend the reach across outcome types — a logistic selector (Pólya–Gamma augmentation carries the whole Gaussian machinery over to a binary outcome) and a survival selector that handles censored time-to-event data — and across scale: a p ≫ n genomics horseshoe, where continuous global-local shrinkage replaces discrete indicators and a specialised sampler makes selection over thousands of predictors routine. Each is cross-checked across Python and R.
Stochastic Search Variable Selection (SSVS)
Part 1 of the variable-selection arc: spike-and-slab, from scratch. Given a pile of candidate predictors, George & McCulloch (1993) attach to each coefficient a latent inclusion indicator and a spike-and-slab prior — a narrow spike pinning switched-off coefficients near zero, a wide slab allowing a real effect — and let a Gibbs sampler wander through the model space (the stochastic search) — a posterior over 2p models. Because the spike is a proper normal, every full conditional is closed form; the indicator update — a Bernoulli from the ratio of slab to spike density at the current coefficient — is the search. Reproducing Congdon's turnip-greens riboflavin example (Draper–Smith), the sampler recovers that soil-moisture tension alone carries the signal (posterior inclusion ~0.69, the "moisture only" model most visited), delivering not a point estimate but posterior inclusion probabilities and a ranked model space. The from-scratch Gibbs and a PyMC compound sampler agree to within a few points — showing why George–McCulloch derived a fully-conjugate scheme for the sharp spike/slab geometry.
View example →Exact Model Enumeration & Bayesian Model Averaging
Part 2 of the variable-selection arc: the conjugate route. Where
Part 1 searched the model space, a small candidate set
can be enumerated exactly. George & McCulloch (1997) put a
conjugate Zellner g-prior on each model so the coefficients integrate out, leaving a
closed-form marginal likelihood that depends on the data only through the model's
R² — an automatic Occam balance of fit against size. Scoring all 2⁹ = 512 models of
the cystic-fibrosis PEmax data gives exact posterior model probabilities, inclusion
probabilities, and model-averaged coefficients — no MCMC. Weight, FEV1 and BMP carry
the signal, yet the best single model holds only ~5% of the posterior: the honest
headline is uncertainty, which stepwise selection hides. The exact enumeration and Part 1's
spike-and-slab search agree on which predictors matter — licensing the search where enumeration is
infeasible. Cross-checked against PyMC and R's BAS.
SSVS for VARs — Bayesian Restriction Search
Part 3 of the variable-selection arc: restriction search on a whole system. An unrestricted VAR is famously over-parameterised, and George, Sun & Ni (2008) turn Part 1's spike-and-slab loose on it — selecting on two fronts at once: which lag coefficients are real, and which contemporaneous links are real (via spike-and-slab on the off-diagonals of the error-precision Cholesky factor). On a US monetary VAR the search discovers, rather than imposes, textbook dynamics — a Taylor-rule-like interest-rate equation — switches off about half the lag terms, and prunes the covariance too (keeping ffr↔inflation and ffr↔output, dropping the direct output↔inflation link). Unlike a shrinkage prior it returns genuine 0/1 restrictions and their probabilities — a posterior over VAR specifications — and it is the only method in the series that selects the error covariance as well as the mean. The selection counterpart to the Minnesota and horseshoe shrinkage priors used elsewhere. Cross-checked against PyMC and base R.
View example →Bayesian Selection of Log-Linear Models
Part 4 of the variable-selection arc: selection for categorical data. For a contingency table, "which predictors matter" becomes "which associations are present". Albert (1996) models the cell counts as Poisson with a log-linear mean, puts robust Cauchy priors on the interaction terms, and compares models by marginal likelihood — so the data, not the prior's tails, decide. The showcase is the Berkeley admissions table, where aggregate figures suggest gender bias: the selection asks whether a direct admission–gender association survives once department is accounted for. It does not — the conditional-independence model carries 90% of the posterior (Admit:Dept and Gender:Dept certain, Admit:Gender only ~0.10) — resolving Simpson's paradox as model selection. The same spike-and-slab / marginal-likelihood logic as Parts 1–3, now over the interactions of a log-linear model. Cross-checked against PyMC and R.
View example →Reversible-Jump MCMC — Inference over Model Dimension
Part 5 of the variable-selection arc — the closing node. Every earlier method chose among models of fixed size; the last question is harder: how many components does the model have? When the parameter's dimension is itself unknown, an ordinary sampler cannot move — the vector changes length. Green's (1995) reversible-jump MCMC constructs birth / death / shift moves that jump between dimensions while preserving detailed balance. On Green's flagship coal-mining disaster series it infers the number of change-points jointly with their locations and rates: a change is certain, the posterior favours two regime shifts (a major drop ~1890 and a smaller one ~1947) — the celebrated "single 1890 change-point" is just the k = 1 slice. Collapsing the conjugate rates lets the jumps move over discrete partitions with no Jacobian. It carries the arc's idea — let the data weigh the models — from selecting coefficients to selecting the model's very size. Cross-checked against PyMC and base R.
View example →Dirichlet-Process Mixtures — How Many Components?
Part 6 of the variable-selection arc: the mixture face of unknown dimension. Where
Part 5 asked "how many change-points?", its sibling asks
"how many components generated this data?" — a mixture whose number of components
K is itself a parameter. On the benchmark galaxy-velocity data, a
from-scratch collapsed Dirichlet-process sampler (Neal's Algorithm 3) integrates out
each Gaussian component analytically and lets a Chinese-restaurant allocation open and close
components — trans-dimensional moves with no Jacobian, the same recipe as Part 5. The
posterior on K peaks at 4 but spreads over 3–7: the honest answer is a
distribution over the number of components, not a single count. The collapsed sampler, a PyMC
stick-breaking DP, and R's mclust (EM + BIC) all agree on four — three philosophies,
one answer — completing the unknown-dimension pair (change-points and mixture
components) that caps the arc.
Free-Knot Splines — an Unknown Number of Knots
Part 7 of the variable-selection arc: nonparametric curve fitting, the third
unknown-dimension problem. A cubic spline bends at its knots — but how many, and
where? DiMatteo, Genovese & Kass (2001) let the number of knots itself
be inferred, the spline analogue of the change-point and
mixture-component samplers. A from-scratch free-knot RJMCMC
integrates out the coefficients and moves knots by birth / death / relocation, then
averages the curve over all sampled configurations. On a curve with a sharp peak the
posterior favours five knots (inferred, not assumed) and migrates them to the curvature
— capturing the peak (RMSE 0.014) where a PyMC/mgcv penalized spline, spreading one
smoothness everywhere, oversmooths it (0.067). Selecting where the flexibility goes beats
applying it uniformly — completing the unknown-dimension trilogy (change-points,
components, knots). Cross-checked against PyMC and R.
The Horseshoe at Scale — Selection when p ≫ n
Part 8 of the variable-selection arc: genomics-scale shrinkage. Real selection is
often thousands of predictors, tens of samples — where
enumeration (2⁴⁰⁸⁸ models) and discrete
spike-and-slab (4,088 indicators) both break down. The
horseshoe (Carvalho–Polson–Scott) instead shrinks: a global scale pulls
everything toward zero while heavy-tailed local scales let a few coefficients escape — continuous,
essentially parameter-free selection. On the riboflavin benchmark (n = 71,
p = 4,088 genes) a from-scratch Gibbs sampler — Makalic–Schmidt augmentation plus the
Bhattacharya O(n²p) coefficient draw — drives the median shrinkage weight to 0.998
and cleanly surfaces the driver genes YOAB and YXLE in ~9 seconds
(vs ~6 minutes for general NUTS). The lesson: when p ≫ n the algorithm, not just the prior,
is the enabling idea — and the reliable output is the shrinkage ranking, not individual near-zero
estimates. Cross-checked against PyMC and the glmnet lasso.
Logistic Variable Selection via Pólya–Gamma SSVS
Part 9 of the variable-selection arc: selection for a binary outcome. Most applied
problems have a yes/no response, which breaks the Gaussian conjugacy the earlier notebooks relied on.
Pólya–Gamma augmentation (Polson–Scott–Windle) fixes it — one latent weight per
observation makes the logistic likelihood conditionally Gaussian, so the exact
same spike-and-slab sweep carries over. On the
South African Heart Disease data both the SSVS search and an exact enumeration of all
2⁹ logistic models keep the same five risk factors — age, family history, tobacco,
type-A, LDL — dropping the body-size measures whose association is only marginal (selection
sees through the confounding). And selection is free: the five-variable model matches
the full model's predictions (r = 0.985) and generalizes slightly better out of sample. Bridges
Part 1 and Part 2 to binary
data. Cross-checked against PyMC and R's BAS.
Survival Variable Selection with Censored Data
Part 10 of the variable-selection arc: time-to-event outcomes, and the complication
no earlier notebook faced — censoring. When a study ends most patients are still
alive, so their survival time is known only to exceed their follow-up. A from-scratch
Weibull proportional-hazards spike-and-slab handles it through the
cumulative-hazard term — subtracted for every subject, so a censored patient
correctly contributes "survived at least this long" — while the
Part-1 inclusion step is untouched. On the Mayo
PBC data it recovers the hand-built natural-history model (bilirubin, age, edema,
albumin, prothrombin — plus stage and copper) and correctly declines the ineffective
treatment (inclusion ≈ 0.10). A parametric Weibull SSVS, PyMC, and R's Cox-based
bic.surv agree — the selection is a property of the proportional-hazards structure, not
the baseline form. This completes the arc's outcome types — continuous, counts,
binary, and now time-to-event. Cross-checked against PyMC and R.