George-McCulloch (1997) Approaches For Bayesian Variable Selection

Summary

A comprehensive treatment of Bayesian variable selection in linear regression that extends George-McCulloch (1993) by (1) embedding stochastic search variable selection (SSVS) within a unified hierarchical mixture model, (2) introducing a conjugate prior formulation that enables closed-form analytical integration of $\beta$ and $\sigma$, (3) deriving Gray Code exhaustive evaluation for small $p$, and (4) comparing three Markov chain Monte Carlo (MCMC) algorithms (nonconjugate Gibbs NG, conjugate Gibbs CG, conjugate Metropolis CM) on simulated problems and a $p=200$ financial index-tracking application. The conjugate prior links $\beta$ and $\sigma^2$ so their marginal may be computed analytically; this eliminates two Gibbs blocks and prevents the stuck-$\gamma$ pathology that afflicts NG under weak information. Published in Statistica Sinica 7 (1997): 339–373; working paper (revised November 1996) by George at UT Austin (Management Science and Information Systems) and McCulloch at U Chicago (Graduate School of Business).

Key Claims

• Unified hierarchical model (§2): $f(Y|\beta,\sigma) = \mathcal{N}_n(X\beta, \sigma^2 I)$; prior $\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}$ (eq. 4). This unifies all George-McCulloch (1993) and (1997) variants under one framework.
• Nonconjugate prior = SSVS (§3): $\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 significance threshold $\delta$ satisfies $\log(v_1/v_0)/(v_0^{-1}-v_1^{-1}) = \delta^2$ (eq. 8).
• Conjugate prior (§4): $\pi(\beta|\sigma,\gamma) = \mathcal{N}_p(0, \sigma^2 D^*_\gamma R_\gamma D^*_\gamma)$ (eq. 15); marginal of $\beta_i$ is $t$ (eq. 19). Translation from nonconjugate: $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 integration (§4.2): With conjugate prior, $\beta$ and $\sigma$ may be integrated out. 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)$ (eq. 24), where $S^2_\gamma = Y'Y - Y'X(X'X+(D^*_\gamma R_\gamma D^*_\gamma)^{-1})^{-1}X'Y$ (eq. 25). G-prior case ($v^*_0 \equiv 0$) gives eqs. 27–28.
• Fast updating (§4.3): Under condition (29) ($v^*_0>0$, $R_\gamma=I$): Chambers (1971) algorithm updates $g(\gamma)$ in $O(p^2)$ when one $\gamma_i$ flips; formula $g(\gamma) = (\prod_i T^2_{ii}[(1-\gamma_i)v^*_0+\gamma_i v^*_1])^{-1/2}(\nu\lambda+Y'Y-W'W)^{-(n+\nu)/2}\pi(\gamma)$ (eq. 31). Under condition (30) (g-prior, $D^*_{1\gamma}R_{1\gamma}D^*_{1\gamma}=c(X'_\gamma X_\gamma)^{-1}$): Dongarra et al. (1979) update in $O(q^2_\gamma)$; $g(\gamma) = (1+c)^{-q_\gamma/2}(\nu\lambda+Y'Y-(1+1/c)^{-1}W'W)^{-(n+\nu)/2}\pi(\gamma)$ (eq. 32).
• Gray Code exhaustive evaluation (§4.4): Order all $2^p$ $\gamma$ values so consecutive ones differ in exactly one bit; apply fast updating at each step. Total $O(2^p \cdot p^2)$ operations vs. brute-force $O(2^p \cdot p^3)$. Feasible for $p \leq 25$ (p=15: 7 sec; p=20: 5 min on 1996 SPARCstation 10; Diaconis-Holmes 1994; Reingold-Nievergelt-Deo 1977).
• Three MCMC algorithms (§5): NG (nonconjugate Gibbs = SSVS), CG (conjugate Gibbs, samples $\gamma_i$ via Bernoulli ratio $g(\gamma_i=1,\gamma_{(i)})/g(\gamma_i=0,\gamma_{(i)})$ eq. 40, never samples $\beta$ or $\sigma$), CM (conjugate Metropolis = Markov chain Monte Carlo model composition, MC³, of Madigan-York 1995, proposes flipping or swapping one $\gamma_i$, eq. 43–44). Speed at $p=100$: NG ≈ 30 sec/1000 iter; CG/CM ≈ 5 min/1000 iter. Normalization constant estimation improves individual probabilities after a MCMC run (eqs. 34–38).
• Stuck-γ pathology of NG (§5.3, simulation 3): Under weak information ($\sigma=200$, $v_1/v_0=2500$), NG gets trapped: $\beta$ draws from $\mathcal{N}(0,v_0)$ are unlikely to look like $\mathcal{N}(0,v_1)$, so $\gamma_i$ stays at 0 for long runs then jumps to 1 and stays. Monte Carlo standard error (SE) for NG ≈ 0.10 vs. CG/CM ≈ 0.02 in this setting.
• Financial application (§6): $p=200$ stocks, $n=362$ weekly S&P 500 returns. Two-phase: Phase 1 (NG, 10K draws, $w_i=0.05$, $\delta=0.008$, $v_1/v_0=625$) identifies top 50 stocks; Phase 2 (CG, 50 stocks, 450K draws, $w_i=0.01$) identifies 789 high-probability models within $e^{-5}$ of best. Best model: 21 variables, $R^2=95.6\%$, matches stepwise regression. Stepwise model posterior probability = $0.00176 \times$ max-probability model.

Concepts Introduced or Extended

• Spike-and-Slab Prior — conjugate formulation, analytical marginalisation of β and σ², stuck-γ pathology, three-algorithm comparison
• Variable Selection — conjugate prior, analytical integration, fast updating, Gray Code, three-algorithm comparison, stuck-γ failure mode, normalization constant estimation

Entities Mentioned

• Edward I. George
• Robert E. McCulloch

Quotes

"We find that the most successful and general-purpose strategy often involves a two-phase approach: a fast initial screening phase using the nonconjugate method to identify a manageable set of candidate predictors, followed by a more accurate conjugate phase for detailed model comparison." (§6, paraphrased)

"When information is weak, the nonconjugate Gibbs sampler can get trapped in a mode: $\beta$ is drawn from the spike and is unlikely to resemble a draw from the slab, so $\gamma_i$ stays at 0 for long stretches." (§5.3, paraphrased)

My Take

The paper's central contribution is showing that the conjugate–nonconjugate distinction matters computationally and inferentially, not just aesthetically. The analytical integration result is elegant: marginalizing $\beta$ and $\sigma$ collapses a three-block Gibbs into a one-block sampler over the discrete $\gamma$ space, and the Gray Code enumeration for $p \leq 25$ delivers exact posterior probabilities for all $2^p$ models. The stuck-$\gamma$ failure mode is underappreciated: it shows that a seemingly reasonable prior can make a Gibbs chain nearly non-ergodic, with the chain's behavior masking the problem. The two-phase S&P 500 application (p=200→50→21 variables) is a convincing demonstration of scalability. One limitation: the conjugate approach requires $O(p^2)$ or $O(q^2_\gamma)$ updates which are slow for very large $p$; the nonconjugate approach trades mixing quality for speed, making the choice of method problem-dependent.