Logistic Variable Selection via Pólya–Gamma SSVS

Python · R · Part 9 of the variable-selection arc — binary outcomes

Overview

Part 9 of the variable-selection arc: selection for a binary outcome. Most applied selection problems have a binary response — disease or not, default or not — which breaks the clean conjugacy the earlier notebooks relied on: the logistic likelihood is not Gaussian, so the spike-and-slab Gibbs of Part 1 does not apply directly. The elegant fix is Pólya–Gamma data augmentation (Polson, Scott & Windle 2013): introduce a latent weight $\omega_i$ per observation and the logistic likelihood becomes conditionally Gaussian in the coefficients — after which the exact same conjugate SSVS sweep works. The showcase is the South African Heart Disease data (the logistic example in Elements of Statistical Learning): 462 men, coronary heart disease against nine risk factors — which of the nine belong in the model? Because $p=9$ is small, we run the SSVS search and enumerate all $2^9=512$ logistic models exactly — bridging Part 1 and Part 2 for a binary outcome — cross-checked against PyMC and R's BAS.

\omega_i\mid\beta \sim \text{PG}(1,\ x_i'\beta) \;\Longrightarrow\; \beta\mid\omega,\gamma \sim \mathcal N \ \text{(weighted LS)}, \quad z_i = \tfrac{y_i-1/2}{\omega_i}

The Pólya–Gamma trick — a non-conjugate likelihood made Gaussian

Logistic regression with a spike-and-slab prior on the standardized predictors: $\Pr(y_i=1)=1/(1+e^{-x_i'\beta})$ with the usual narrow-spike / wide-slab prior and $\gamma_j\sim\text{Bernoulli}(w)$ . The likelihood has no conjugate prior, so $\beta\mid\text{data}$ is not a tidy normal and the Part-1 Gibbs stalls. Polson–Scott–Windle show that with a latent Pólya–Gamma variable $\omega_i\sim\text{PG}(1, x_i'\beta)$ per observation, the likelihood contribution becomes Gaussian in $\beta$ with pseudo-response $(y_i-\tfrac12)/\omega_i$ and weight $\omega_i$ . So the Gibbs sweep is: (1) draw $\omega_i\mid\beta\sim\text{PG}(1, x_i'\beta)$ with the Devroye/PSW exact sampler; (2) draw $\beta\mid\omega,\gamma$ as a weighted-least-squares normal with the spike-and-slab prior variances, exactly as in Part 1; (3) draw the indicators $\gamma_j$ from the spike/slab density ratio. And exactly: with only nine predictors we also enumerate all $2^9$ models, each logistic model's marginal likelihood by a Laplace approximation with the binomial Fisher information $X'WX$ — giving exact posterior model and inclusion probabilities. Engine: a self-contained logitsel.py.

Results

Five of nine risk factors survive. Both the Pólya–Gamma search and the exact enumeration keep age, family history, tobacco, type-A behaviour and LDL cholesterol, and drop systolic BP, adiposity, obesity and alcohol — exactly the model Elements of Statistical Learning reaches by stepwise selection, recovered here as a posterior over models with uncertainty attached. age is certain (inclusion ≈ 1), family history nearly so, the others firmly in but graded. A famous non-finding: selection is about conditional, not marginal, relevance — and here the two genuinely part ways. Adiposity, on its own about as correlated with heart disease as the retained LDL ( $r\approx0.25$ ), is dropped because its signal is already carried by the factors it travels with (age, obesity, LDL); type-A behaviour, only weakly correlated marginally ( $r\approx0.10$ ), is kept for the independent signal it adds — the logistic posterior sees through the confounding. The selected model is also interpretable: coronary risk climbs ~5% in odds per year of age, and each retained factor (family history most) shifts the whole age–risk curve up additively on the log-odds scale, tracking the observed proportions. Selection is free here: the five-variable model's predicted probabilities track the full nine-variable model's almost exactly ( $r=0.985$ ) and generalize slightly better out of sample (CV AUC 0.782 vs 0.770) — the dropped factors were contributing variance, not signal. Three ways agree: the from-scratch Gibbs, PyMC's native logistic SSVS, and exact enumeration all select the same five — the augmentation is what lets the whole Gaussian selection machinery of Part 1 carry over to binary data (and, unchanged, to counts and multinomial beyond).

Notebooks

Python Logistic SSVS from scratch — Pólya–Gamma augmentation turning the logistic likelihood conditionally Gaussian, run as both a spike-and-slab search and exact enumeration of all 2⁹ models, with a PyMC cross-check, on the South African heart data (Python) Open → R R cross-check — BAS Bayesian logistic model averaging (R) Open →

Downloads

logitsel.py The from-scratch logistic-selection engine — Pólya–Gamma SSVS Gibbs sampler and exact Laplace enumeration of all 2ᵖ logistic models (new self-contained module) saheart_data.csv South African Heart Disease — 462 men, 9 risk factors, binary CHD (bestglm::SAheart; Rousseauw et al. 1983)

References

Polson, N. G., Scott, J. G. & Windle, J. (2013). Bayesian inference for logistic models using Pólya–Gamma latent variables. Journal of the American Statistical Association 108(504), 1339–1349. — the Pólya–Gamma augmentation reproduced here
George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. Journal of the American Statistical Association 88(423), 881–889. — the spike-and-slab SSVS sweep carried over to the binary case
Hastie, T., Tibshirani, R. & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.), §4.4.2. Springer. — the South African heart-disease logistic example and stepwise selection
Rousseauw, J. et al. (1983). Coronary risk factor screening in three rural communities. South African Medical Journal 64, 430–436. — the source of the heart-disease data
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 cross-check
Barbieri, M. M. & Berger, J. O. (2004). Optimal predictive model selection. Annals of Statistics 32(3), 870–897. — the median-probability model and inclusion probabilities