Variable Selection

Definition

Variable selection is the problem of identifying which subset of $p$ candidate regressors $X_1,\ldots,X_p$ should enter a linear model $Y = X_\gamma\beta_\gamma + \varepsilon$ , trading off model fit against complexity across the $2^p$ possible subsets. Bayesian approaches assign a prior over subsets and use posterior inference to identify high-probability models.

Key Ideas

The model space has $2^p$ elements; for $p=60$ this exceeds $10^{18}$ , making direct enumeration impossible. Markov Chain Monte Carlo (MCMC) navigates the space by visiting high-probability regions.
A binary indicator vector $\gamma \in \{0,1\}^p$ encodes which variables are included ( $\gamma_i=1$ ) or excluded ( $\gamma_i=0$ ).
The prior on $\gamma$ encodes the researcher's beliefs about model size; the simplest choice $P(\gamma_i=1)=p_i=0.5$ gives $f(\gamma)=2^{-p}$ (indifference), implying equal prior probability to all models.
High-frequency $\gamma$ values in MCMC output identify high-posterior-probability subsets without requiring explicit enumeration.

How It Works

Stochastic Search Variable Selection (George-McCulloch 1993)

The SSVS formulation links each $\beta_i$ to a latent inclusion indicator $\gamma_i$ via a mixture prior:

$\beta_i|\gamma_i \sim (1-\gamma_i)\,\mathcal{N}(0,\,\tau_i^2) + \gamma_i\,\mathcal{N}(0,\,c_i^2\tau_i^2) \tag{eq. 2}$

with $\tau_i$ small (exclusion: spike) and $c_i \gg 1$ (inclusion: slab). This does not hard-zero $\beta_i$ when $\gamma_i=0$ ; instead it shrinks it toward zero by a factor of $c_i$ . The joint prior is:

$\beta|\gamma \sim \mathcal{N}_p(0,\, D_\gamma R D_\gamma) \tag{eq. 4}$

where $D_\gamma = \mathrm{diag}[a_1\tau_1,\ldots,a_p\tau_p]$ with $a_i = 1$ if $\gamma_i=0$ and $a_i=c_i$ if $\gamma_i=1$ (eq. 5). The error variance prior is $\sigma^2|\gamma \sim \mathrm{IG}(\nu_\gamma/2,\, \nu_\gamma\lambda_\gamma/2)$ (inverse-gamma (IG), eq. 6).

Three-block Gibbs sampler (given current draws $(\beta^{j-1},\sigma^{j-1},\gamma^{j-1})$ ):

Block 1 — $\beta^j$ : Conjugate Normal update (eq. 13): $\beta^j \sim \mathcal{N}_p\!\left(A_\gamma(\sigma^{j-1})^{-2}X'X\hat\beta_{LS},\; A_\gamma\right)$ $A_\gamma = \left((\sigma^{j-1})^{-2}X'X + D_\gamma^{-1}R^{-1}D_\gamma^{-1}\right)^{-1}$

Block 2 — $\sigma^j$ : Conjugate Inverse-Gamma update (eq. 14): $\sigma^j \sim \mathrm{IG}\!\left(\frac{n+\nu_\gamma}{2},\; \frac{\|Y-X\beta^j\|^2 + \nu_\gamma\lambda_\gamma}{2}\right)$

Block 3 — $\gamma_i^j$ : Independent Bernoulli draws componentwise (eqs. 15–16). For each $i$ : $P(\gamma_i=1|\beta^j, \gamma_{(i)}) = \frac{a}{a+b}$ where (eqs. 16c–d): $a = f(\beta^j|\gamma_{(i)},\, \gamma_i=1)\cdot p_i, \qquad b = f(\beta^j|\gamma_{(i)},\, \gamma_i=0)\cdot(1-p_i)$ Each $f(\cdot)$ is a univariate Normal density evaluation at $\beta_i^j$ under the spike or slab component — the $\gamma_i$ conditional does not depend on $Y$ directly, only through the already-drawn $\beta^j$ .

Tuning Parameters

The method requires choices for $(\tau_i, c_i, p_i)$ :

$\tau_i$ should be small enough that $\beta_i \in \mathcal{N}(0,\tau_i^2)$ is effectively zero in the regression.
$c_i>1$ controls the slab width relative to the spike; larger $c_i$ sharpens separation.
Calibration: $\hat\sigma_{\hat\beta_i}/\tau_i$ (signal-to-noise of each coefficient relative to the spike). Suggested choices: $(\hat\sigma_{\hat\beta_i}/\tau_i, c_i) \in \{(1,5),(1,10),(10,100),(10,500)\}$ .
Signal ratio: $r_i = \sqrt{(\hat\sigma^2_{\hat\beta_i}/\tau_i^2 + c_i^2)/(\hat\sigma^2_{\hat\beta_i}/\tau_i^2 + 1)}$ (eq. 9) — measures the Bayes factor separation between inclusion and exclusion for coefficient $i$ .
Correlation: $R=I$ (independence) or $R \propto (X'X)^{-1}$ (Zellner g-prior structure).

Unified Hierarchical Framework (George-McCulloch 1997)

George-McCulloch (1997) embeds SSVS within a unified hierarchical mixture model:

$f(Y|\beta,\sigma) = \mathcal{N}_n(X\beta,\sigma^2 I), \quad \pi(\beta,\sigma,\gamma) = \pi(\beta|\sigma,\gamma)\pi(\sigma|\gamma)\pi(\gamma)$

$\pi(\gamma) = \prod_i w_i^{\gamma_i}(1-w_i)^{1-\gamma_i} \tag{eq. 4}$

Two prior formulations emerge from this framework:

Nonconjugate prior (= SSVS): $\pi(\beta|\gamma) = \mathcal{N}_p(0,D_\gamma R_\gamma D_\gamma)$ (eq. 5); $\beta$ and $\sigma$ are independent given $\gamma$ . Point-mass limit ( $v_0 \to 0$ ): $\pi(\beta_i|\gamma) = (1-\gamma_i)I_0 + \gamma_i\mathcal{N}(0,v_{1_{\gamma(i)}})$ (eq. 9) — makes Gibbs reducible. Practical $\delta$ -threshold: $\log(v_1/v_0)/(v_0^{-1}-v_1^{-1}) = \delta^2$ (eq. 8).

Conjugate prior (new in 1997): $\pi(\beta|\sigma,\gamma) = \mathcal{N}_p(0,\sigma^2 D^*_\gamma R_\gamma D^*_\gamma)$ (eq. 15); $\beta$ and $\sigma$ are not independent given $\gamma$ . Marginal of $\beta_i$ is $t$ -distributed (eq. 19). Translation from nonconjugate scale: $v^*_0 = v_0/\hat\sigma^2$ , $v^*_1 = v_1/\hat\sigma^2$ (eq. 18). Point-mass limit: $\pi(\beta_i|\sigma,\gamma_{(i)}) = (1-\gamma_i)I_0 + \gamma_i T(\nu,0,\lambda_\gamma v^*_{1_{\gamma(i)}})$ (eq. 21).

Analytical Marginal Posterior Under Conjugate Prior (George-McCulloch 1997)

With the conjugate prior, $\beta$ and $\sigma^2$ can be integrated out analytically, leaving a posterior over the discrete $\gamma$ space alone. When $v^*_0 > 0$ :

$\pi(\gamma|Y) \propto g(\gamma) \equiv |\tilde X'\tilde X|^{-1/2}\,|D^*_\gamma R_\gamma D^*_\gamma|^{-1/2}\,(\nu\lambda+S^2_\gamma)^{-(n+\nu)/2}\,\pi(\gamma) \tag{eq. 24}$

$S^2_\gamma = Y'Y - Y'X\!\left(X'X+(D^*_\gamma R_\gamma D^*_\gamma)^{-1}\right)^{-1}X'Y \tag{eq. 25}$

The g-prior case ( $v^*_0 \equiv 0$ ) yields a simplified formula (eqs. 27–28). This reduces the three-block Gibbs to a one-block sampler over the $2^p$ discrete $\gamma$ space — a major computational simplification.

Fast Updating and Gray Code Exhaustive Evaluation (George-McCulloch 1997)

Chambers (1971) condition (29): When $v^*_0>0$ and $R_\gamma=I$ , flipping one $\gamma_i$ updates $g(\gamma)$ in $O(p^2)$ (Chambers rank-1 update):

$g(\gamma) = \left(\prod_i T^2_{ii}[(1-\gamma_i)v^*_0+\gamma_i v^*_1]\right)^{-1/2}(\nu\lambda+Y'Y-W'W)^{-(n+\nu)/2}\pi(\gamma) \tag{eq. 31}$

Dongarra et al. (1979) condition (30): When $D^*_{1\gamma}R_{1\gamma}D^*_{1\gamma}=c(X'_\gamma X_\gamma)^{-1}$ (g-prior structure), update is $O(q^2_\gamma)$ :

$g(\gamma) = (1+c)^{-q_\gamma/2}\!\left(\nu\lambda+Y'Y-(1+1/c)^{-1}W'W\right)^{-(n+\nu)/2}\pi(\gamma) \tag{eq. 32}$

Gray Code enumeration: Order all $2^p$ $\gamma$ vectors so consecutive ones differ in exactly one bit (Reingold-Nievergelt-Deo 1977; Diaconis-Holmes 1994). Apply fast $O(p^2)$ updating at each step to enumerate $g(\gamma)$ for all models in $O(2^p \cdot p^2)$ total operations (vs. brute force $O(2^p \cdot p^3)$ ). Feasible through $p \leq 25$ ( $p=15$ : 7 sec; $p=20$ : 5 min on a 1996 SPARCstation 10).

Three MCMC Algorithms and Normalization Constant Estimation (George-McCulloch 1997)

NG — Nonconjugate Gibbs (= SSVS): Full three-block Gibbs over $(\beta,\sigma^2,\gamma)$ . Fastest per iteration at large $p$ ( $\approx 30$ sec/1000 iters at $p=100$ ) but vulnerable to stuck- $\gamma$ failure.

CG — Conjugate Gibbs: Samples only $\gamma$ using the Bernoulli inclusion ratio:

$\frac{\pi(\gamma_i=1,\gamma_{(i)}|Y)}{\pi(\gamma_i=0,\gamma_{(i)}|Y)} = \frac{g(\gamma_i=1,\gamma_{(i)})}{g(\gamma_i=0,\gamma_{(i)})} \tag{eq. 40}$

Never samples $\beta$ or $\sigma$ directly. Slower per iteration at large $p$ ( $\approx 5$ min/1000 iters at $p=100$ ) but mixes correctly even under weak information.

CM — Conjugate Metropolis (= MC³ of Madigan-York 1995): Proposes flipping or swapping one $\gamma_i$ ; accepts via Metropolis ratio using $g(\gamma)$ (eqs. 43–44). Similar speed to CG; explores model space differently.

Normalization constant estimation: Since $\pi(\gamma|Y) = C \cdot g(\gamma)$ (eq. 34), estimate $\hat C$ from MCMC output via indicator function on a reference set $A$ :

$\hat C = \frac{1}{g(A) \cdot K}\sum_k I_A(\gamma^{(k)}) \tag{eq. 35}$

Improved individual posterior estimates $\hat\pi(\gamma|Y) = \hat C \cdot g(\gamma)$ (eq. 36); total visited probability $\hat\pi(B|Y) = \hat C \cdot g(B)$ (eq. 37) characterizes model space coverage.

Stuck- $\gamma$ Failure Mode of Nonconjugate Gibbs

Under weak information (e.g., $\sigma=200$ , $v_1/v_0=2500$ ), NG malfunctions: $\beta$ is drawn from $\mathcal{N}(0,v_0)$ (the spike) and is unlikely to resemble a draw from $\mathcal{N}(0,v_1)$ (the slab). Consequently $\gamma_i$ stays at 0 for long runs, occasionally jumps to 1 and stays there. The $\gamma$ chain is nearly non-ergodic. Monte Carlo SE for marginal inclusion probabilities: NG ≈ 0.10 vs. CG/CM ≈ 0.02 in this setting. The conjugate methods, which never sample $\beta$ , are immune.

Financial Index Tracking Application (George-McCulloch 1997)

$p=200$ S&P 500 constituent stocks, $n=362$ weekly returns. Two-phase strategy:

Phase 1 (NG screening): 10,000 draws, $w_i=0.05$ , $\delta=0.008$ , $v_1/v_0=625$ . Identifies top 50 stocks by marginal inclusion frequency.

Phase 2 (CG estimation): 50 stocks, 450,000 draws, $w_i=0.01$ . Identifies 789 high-probability models within $e^{-5}$ of the best model. Best model: 21 variables, $R^2=95.6\%$ , matches classical stepwise regression. The stepwise model's posterior probability is only $0.00176\times$ the maximum-probability model — demonstrating that a single "best" model severely understates posterior uncertainty.

Why It Matters

Provides a principled Bayesian approach to subset selection that avoids multiple comparisons issues and inherently accounts for model uncertainty.
The "all post-burn-in" strategy — using every Gibbs draw, not just the last — gives a direct read on the marginal posterior inclusion probability $P(\gamma_i=1|Y)$ for each variable.
Scales to $p=60$ (SSVS) and $p=200$ (two-phase CG strategy) by exploiting the J-shaped frequency distribution: most models are never visited because the chain concentrates on high-probability subsets.
Gray Code exhaustive evaluation delivers exact posterior probabilities for all $2^p$ models at $p \leq 25$ , making Bayesian model averaging exact in moderate dimensions.
The stuck- $\gamma$ failure mode illustrates a general principle: a discrete sampling block (Block 3) that depends on a continuous draw (Block 1 for $\beta$ ) can have its mixing controlled entirely by the scaling of the continuous prior — independent of the data.
Foundational precursor to modern spike-and-slab, horseshoe, and Bayesian LASSO priors for high-dimensional regression; conjugate formulation prefigures the integrated-likelihood approaches used in BAS and BMA software.

Open Questions

Tuning $(v_0, v_1, w_i)$ remains problem-specific; data-dependent calibration can compromise the prior's external validity.
The $O(p^2)$ fast update cost is manageable through $p \sim 100$ but becomes prohibitive for very large $p$ (e.g. $p \sim 10{,}000$ ).
CG/CM are slower per iteration than NG at large $p$ ; the optimal algorithm choice depends on the signal-to-noise ratio of the data, which is unknown in advance.
Modern continuous shrinkage priors (horseshoe, Dirichlet-Laplace) avoid the $\gamma$ block entirely, reducing tuning burden at some cost in interpretability and exact model probabilities.
Normalization constant estimation (§4.5) requires that the MCMC run visits a reference set frequently enough to estimate $\hat C$ reliably — an implicit convergence requirement that is hard to verify.

Bayesian Stochastic Search for VAR Restrictions (George-Sun-Ni 2008)

George, Sun, and Ni extend SSVS to VAR models, simultaneously searching restrictions on the $p \times (h+Lp)$ coefficient matrix $\Phi$ and the off-diagonal elements of the upper-triangular factor $\Psi$ in the precision decomposition $\Sigma^{-1} = \Psi\Psi'$ (eq. 9). Upper-triangular $\Psi$ is required for global identification (unique $\Sigma \to \Psi$ mapping); arbitrary $\Psi$ satisfies only local rank conditions (Bekker-Pollock 1986; Amisano-Giannini 1997).

Prior structure: Five-component mixture prior.

Always-included $\phi_{non}$ : $\mathcal{N}_{n-m}(\tilde\phi_{non}, M_{non})$ (eq. 11).
Restricted $\phi_{s_i}$ : $(1-\gamma_i)\mathcal{N}(0,\tau_{0i}^2)+\gamma_i\mathcal{N}(0,\tau_{1i}^2)$ ; $\gamma_i \sim \text{Bernoulli}(p_i)$ (eqs. 13–14).
Off-diagonal $\psi_{ij}$ : $(1-\omega_{ij})\mathcal{N}(0,\kappa_{0ij}^2)+\omega_{ij}\mathcal{N}(0,\kappa_{1ij}^2)$ ; $\omega_{ij} \sim \text{Bernoulli}(q_{ij})$ (eqs. 16–17).
Diagonal $\psi_{ii}^2 \sim \text{Gamma}(a_i,b_i)$ — always included, guaranteeing $\Sigma^{-1} = \Psi\Psi'$ is positive definite (eq. 18).
Semiautomatic tuning: $\tau_{ki} = c_k\hat\sigma_{\phi_{s_i}}$ , $c_0=0.1$ , $c_1=10$ (analogous to George-McCulloch (1997)).

Five-step Gibbs sampler (Appendix A.3): Every conditional is standard — Gamma for $\psi_{ii}^2$ , Normal for $\eta_j$ , Bernoulli for $\omega_{ij}$ and $\gamma_i$ , Normal for $\phi$ — with no Metropolis-Hastings steps. Key property inherited from George-McCulloch (1993): $\gamma_i$ does not depend on $Y$ given $\phi$ (eq. 37).

Key finding: Restricting $\Phi$ aids $\Psi$ selection and vice versa — the two searches are interdependent. Rao-Blackwell model-averaged forecasts improve over maximum likelihood estimator (MLE) by 20–67% in simulation.

Poisson Log-Linear Model Selection (Albert 1996)

Albert (1996) extends the SSVS idea to contingency tables via a Poisson log-linear likelihood. The design matrix $X_0$ encodes which interaction terms are included in model $M_k$ . The block-diagonal precision prior is:

$\Sigma^{-1} = \text{diag}(0I_{p_0},\, P_1I_{p_1},\,\ldots,\, P_sI_{p_s})$

where $\beta_0$ (always-included terms) receives an improper flat prior ( $0$ precision) and each group $j$ of testable parameters receives precision $P_j$ . Rather than a fixed spike-and-slab, a Gamma hyperprior $P_j \sim \text{Gamma}(\nu_j/2, b_j^2\nu_j/2)$ integrates out to a $t_{\nu_j}$ marginal for $\beta_j$ . Setting $\nu_j=1$ gives Cauchy marginals — heavier-tailed and more robust to extreme cell counts than the Normal ( $\nu_j\to\infty$ ).

Posterior model probabilities are approximated by a Laplace expansion around the posterior mode found by IWLS. The Gibbs sampler cycles over:

$\beta_j|\text{rest}$ : conjugate Normal given current $P_j$
$Z|\text{rest}$ : auxiliary IWLS working response
$P_j^*|\text{rest}$ : Gamma( $\nu_j/2+p_j/2$ ,; $b_j^2\nu_j/2+\sum_t\beta_{jt}^2/2$ )

Bayesian model averaging (BMA) with this sampler inflates posterior standard errors nearly 2× relative to the maximum-probability single model, illustrating the cost of ignoring model uncertainty.