Binary logistic selection — which risk factors matter?

Bayesian variable selection, Part 9 — Polya-Gamma SSVS

Most applied selection problems have a binary outcome — disease or not, default or not. That 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 Polya-Gamma data augmentation (Polson, Scott & Windle 2013): introduce a latent weight $\omega_i$ for each 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 (Rousseauw et al. 1983; the logistic-regression example in Hastie-Tibshirani-Friedman's Elements of Statistical Learning): 462 men, coronary heart disease (yes/no) against nine risk factors. Which of the nine belong in the model? Because $p=9$ is small, we can run the SSVS search and enumerate all $2^9=512$ logistic models exactly — bridging Parts 1 and 2 for a binary outcome — and cross-check both against PyMC and R's BAS.

The predictors: sbp (systolic BP), tobacco (cumulative), ldl (cholesterol), adiposity, famhist (family history), typea (type-A behaviour), obesity, alcohol, age.

The model and the Polya-Gamma trick

Logistic regression with a spike-and-slab prior on the standardized predictors' coefficients: $$ \Pr(y_i=1) = \frac{1}{1+e^{-x_i'\beta}},\qquad \beta_j\mid\gamma_j=0\sim\mathcal N(0,\tau^2),\quad \beta_j\mid\gamma_j=1\sim\mathcal N(0,c^2\tau^2),\quad \gamma_j\sim\text{Bernoulli}(w). $$

Why it is hard, and the fix. The logistic 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 Polya-Gamma variable $\omega_i\sim\text{PG}(1,\,x_i'\beta)$ per observation, the likelihood contribution becomes $$ \propto \exp\!\Big(-\tfrac{\omega_i}{2}\big(x_i'\beta - \tfrac{y_i-1/2}{\omega_i}\big)^2\Big), $$ i.e. Gaussian in $\beta$ with "response" $(y_i-\tfrac12)/\omega_i$ and weight $\omega_i$. So the Gibbs sweep is:

1. $\omega_i\mid\beta \sim \text{PG}(1, x_i'\beta)$ — drawn with the Devroye/PSW exact sampler;
2. $\beta\mid\omega,\gamma \sim \mathcal N$ — a weighted-least-squares draw with the spike-and-slab prior variances, exactly as in Part 1;
3. $\gamma_j\mid\beta_j \sim \text{Bernoulli}$ — the spike/slab density ratio.

The intercept is always included. Hyperparameters: spike $\tau=0.1$, slab-to-spike ratio $c=20$ (slab sd $=2$ on the log-odds scale), prior inclusion $w=0.5$.

And exactly. With only nine predictors we also enumerate all $2^9$ models: each logistic model's marginal likelihood by a Laplace approximation (Poisson-style, but with the binomial Fisher information $X'\!WX$, $W=\hat\mu(1-\hat\mu)$), giving exact posterior model and inclusion probabilities. Engine: logitsel.py.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
import logitsel as L

d = pd.read_csv("saheart_data.csv"); names = [c for c in d.columns if c != "chd"]
X = d[names].values; y = d["chd"].values                           # 462 men, 9 risk factors, binary CHD

ss = L.pg_ssvs(y, X, tau=0.1, c=20.0, w=0.5, ndraw=4000, burn=2000, seed=1)   # Polya-Gamma SSVS
en = L.enumerate_logit(y, X, names)                                # exact: all 2^9 = 512 models

print("PG-SSVS (Polya-Gamma Gibbs) inclusion probabilities:")
for nm, pi in sorted(zip(names, ss["incl"]), key=lambda t: -t[1]): print("  %-10s %.3f" % (nm, pi))
print("\nexact enumeration (512 logistic models):")
print("  " + "   ".join("%s %.2f" % (nm, pi) for nm, pi in sorted(zip(names, en["incl"]), key=lambda t: -t[1])[:5]) + "    | rest < 0.13")
print("top model:  %s   (posterior probability %.3f)" % (" + ".join(en["top"][0][0]), en["top"][0][1]))

mcorr = np.array([np.corrcoef(X[:, j], y)[0, 1] for j in range(len(names))])   # marginal point-biserial corr with CHD
print("\nmarginal corr with CHD vs. inclusion (selection is conditional, not marginal):")
for nm, mc, pi in sorted(zip(names, mcorr, en["incl"]), key=lambda t: -t[1]):
    print("  %-10s marginal r = %+.2f    inclusion = %.2f" % (nm, mc, pi))

# --- Fig: the ranked logistic model space ---
top = en["top"][:6]                                                # top 6 of 512 models, by posterior probability
labels = [" + ".join(nms) for nms, w in top]; probs = [w for nms, w in top]
ypos = np.arange(len(top))[::-1]
fig, ax = plt.subplots(figsize=(8.4, 4.4))
ax.barh(ypos, probs, color="firebrick", alpha=.85)
for yp, pr in zip(ypos, probs): ax.text(pr + 0.004, yp, "%.3f" % pr, va="center", fontsize=8)
ax.set_yticks(ypos); ax.set_yticklabels(labels, fontsize=9)
ax.set_xlim(0, max(probs) * 1.15); ax.set_xlabel("posterior model probability")
ax.set_title("Ranked logistic model space  (exact enumeration, top 6 of 512)")
ax.grid(axis="x", alpha=.25)
plt.tight_layout(); plt.show()

PG-SSVS (Polya-Gamma Gibbs) inclusion probabilities:
  age        0.998
  famhist    0.884
  typea      0.637
  tobacco    0.582
  ldl        0.553
  obesity    0.123
  adiposity  0.101
  sbp        0.092
  alcohol    0.069

exact enumeration (512 logistic models):
  age 1.00   famhist 1.00   tobacco 0.92   typea 0.90   ldl 0.86    | rest < 0.13
top model:  tobacco + ldl + famhist + typea + age   (posterior probability 0.485)

marginal corr with CHD vs. inclusion (selection is conditional, not marginal):
  age        marginal r = +0.37    inclusion = 1.00
  tobacco    marginal r = +0.30    inclusion = 0.92
  famhist    marginal r = +0.27    inclusion = 1.00
  ldl        marginal r = +0.26    inclusion = 0.86
  adiposity  marginal r = +0.25    inclusion = 0.08
  sbp        marginal r = +0.19    inclusion = 0.09
  typea      marginal r = +0.10    inclusion = 0.90
  obesity    marginal r = +0.10    inclusion = 0.12
  alcohol    marginal r = +0.06    inclusion = 0.05

Reading the selection

• Five of nine risk factors survive. Both the Polya-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. This is exactly the model Elements of Statistical Learning arrives at by stepwise selection — recovered here as a posterior over models, with uncertainty attached.
• Age is certain, the others graded. age enters every good model (inclusion $\approx1$); famhist nearly so; tobacco, typea, ldl are firmly in but with more posterior doubt. The four dropped factors sit near or below the neutral 0.5 line.
• A famous non-finding. Selection is about conditional, not marginal, relevance — and the two genuinely come apart here (see the marginal-vs-inclusion table above). adiposity, on its own about as correlated with CHD as the retained ldl ($r\approx0.25$), is dropped because its signal is already carried by the factors it travels with (age, obesity, ldl); meanwhile typea, only weakly correlated marginally ($r\approx0.10$), is kept for the independent signal it adds. The logistic posterior sees through the confounding.
• The model space is diffuse (figure). Even the top model, tobacco+ldl+famhist+typea+age, holds under half the posterior mass — as in Part 2, "which model" is genuinely uncertain, which is the case for reporting inclusion probabilities rather than a single winner.

import pymc as pm                                                    # cross-check: logistic SSVS in a PPL (no Polya-Gamma needed)
Xz = (X - X.mean(0)) / X.std(0); p = len(names)
with pm.Model() as mod:
    tau, c = 0.1, 20.0
    g  = pm.Bernoulli("g", 0.5, shape=p); b = pm.Normal("b", 0, 1, shape=p)
    beta = b * pm.math.switch(g, c*tau, tau); b0 = pm.Normal("b0", 0, 10)
    pm.Bernoulli("y", logit_p=b0 + pm.math.dot(Xz, beta), observed=y)
    idata = pm.sample(2000, tune=1500, chains=4, step=[pm.BinaryGibbsMetropolis([g])],
                      target_accept=0.9, random_seed=1, progressbar=False)
pm_incl = idata.posterior["g"].mean(("chain","draw")).values
print("PyMC logistic SSVS inclusion: ", "  ".join("%s %.2f" % (n, v)
      for n, v in sorted(zip(names, pm_incl), key=lambda t: -t[1])[:5]), " | rest < 0.12")
print("-> all three engines select the SAME five risk factors and drop sbp, adiposity, obesity, alcohol.")

# --- Fig: inclusion probabilities across all three engines ---
order = np.argsort(-en["incl"])                                    # rank predictors by exact inclusion
labs = [names[j] for j in order]; x = np.arange(len(names)); wbar = 0.27
fig, ax = plt.subplots(figsize=(8.6, 5.0))
ax.bar(x - wbar, ss["incl"][order], wbar, color="firebrick", label="PG-SSVS (Polya-Gamma)")
ax.bar(x,        en["incl"][order], wbar, color="navy",      label="exact enumeration (512)")
ax.bar(x + wbar, pm_incl[order],    wbar, color="gray",      label="PyMC SSVS")
ax.axhline(0.5, color="k", ls="--", lw=1, alpha=.6)               # neutral inclusion reference
ax.set_xticks(x); ax.set_xticklabels(labs, rotation=45, ha="right")
ax.set_ylim(0, 1.05); ax.set_ylabel("inclusion probability")
ax.set_title("Inclusion probabilities across three engines")
ax.legend(); ax.grid(axis="y", alpha=.25)
plt.tight_layout(); plt.show()

Multiprocess sampling (4 chains in 4 jobs)

CompoundStep

>BinaryGibbsMetropolis: [g]

>NUTS: [b, b0]

Sampling 4 chains for 1_500 tune and 2_000 draw iterations (6_000 + 8_000 draws total) took 11 seconds.

C:\Users\user\anaconda3\envs\pymc-env\Lib\site-packages\arviz_stats\base\diagnostics.py:90: RuntimeWarning: invalid value encountered in scalar divide
  (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

PyMC logistic SSVS inclusion:  age 1.00  famhist 0.94  tobacco 0.61  typea 0.58  ldl 0.57  | rest < 0.12
-> all three engines select the SAME five risk factors and drop sbp, adiposity, obesity, alcohol.

Does dropping four factors cost anything? The predictions say no.

A selection result is only worth having if the smaller model is as good. So compare the predicted CHD probabilities of the best 5-variable model against the full 9-variable model, person by person, and check how each generalizes out of sample (10-fold cross-validation).

import statsmodels.api as sm
from sklearn.model_selection import StratifiedKFold
from sklearn.metrics import roc_auc_score, log_loss
best = ["tobacco","ldl","famhist","typea","age"]; bi = [names.index(b) for b in best]
def fit_pred(cols, Xtr, ytr, Xte):                                  # fit a logistic model, return predicted probs
    m = sm.GLM(ytr, sm.add_constant(Xtr[:, cols]), family=sm.families.Binomial()).fit()
    return m.predict(sm.add_constant(Xte[:, cols], has_constant="add"))

pf = fit_pred(list(range(9)), Xz, y, Xz); pb = fit_pred(bi, Xz, y, Xz)    # in-sample predicted P(CHD)
print("predicted-probability correlation, best vs full: %.3f" % np.corrcoef(pf, pb)[0, 1])
skf = StratifiedKFold(10, shuffle=True, random_state=1); af, ab, lf, lb = [], [], [], []
for tr, te in skf.split(Xz, y):                                    # out-of-sample: 10-fold CV
    af.append(roc_auc_score(y[te], fit_pred(list(range(9)), Xz[tr], y[tr], Xz[te])))
    ab.append(roc_auc_score(y[te], fit_pred(bi, Xz[tr], y[tr], Xz[te])))
print("10-fold CV AUC:  full %.3f   best %.3f   (parsimony is free -- indeed slightly better)"
      % (np.mean(af), np.mean(ab)))

# --- Fig: best (5-var) vs full (9-var) predictions, in-sample and cross-validated ---
pfv = np.asarray(pf); pbv = np.asarray(pb)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12.4, 4.7))
ax1.scatter(pfv, pbv, s=16, color="firebrick", alpha=.5, edgecolor="none")     # person-by-person P(CHD)
lim = [0, max(pfv.max(), pbv.max()) * 1.02]
ax1.plot(lim, lim, "k--", lw=1, alpha=.6, label="y = x")
ax1.set_xlim(lim); ax1.set_ylim(lim)
ax1.set_xlabel("full 9-variable  P(CHD)"); ax1.set_ylabel("best 5-variable  P(CHD)")
ax1.set_title("In-sample predictions agree  (r = %.3f)" % np.corrcoef(pfv, pbv)[0, 1])
ax1.legend(); ax1.grid(alpha=.25)
bp = ax2.boxplot([af, ab], widths=.55, patch_artist=True)                       # 10-fold CV AUC spread
for patch, col in zip(bp["boxes"], ["navy", "firebrick"]): patch.set_facecolor(col); patch.set_alpha(.35)
for md in bp["medians"]: md.set_color("k")
ax2.scatter(np.full(len(af), 1), af, color="navy", zorder=3, s=20)
ax2.scatter(np.full(len(ab), 2), ab, color="firebrick", zorder=3, s=20)
ax2.axhline(np.mean(af), color="navy", ls=":", lw=1, alpha=.7)
ax2.axhline(np.mean(ab), color="firebrick", ls=":", lw=1, alpha=.7)
ax2.set_xticks([1, 2]); ax2.set_xticklabels(["full (9)", "best (5)"])
ax2.set_ylabel("out-of-sample AUC")
ax2.set_title("10-fold CV AUC  (full %.3f  vs  best %.3f)" % (np.mean(af), np.mean(ab)))
ax2.grid(axis="y", alpha=.25)
plt.tight_layout(); plt.show()

predicted-probability correlation, best vs full: 0.985
10-fold CV AUC:  full 0.770   best 0.782   (parsimony is free -- indeed slightly better)

Reading the selected model: how risk rises with age

Inclusion probabilities say which factors matter; the fitted model says how. Since age enters every good model, a natural summary is the predicted CHD probability as a function of age, drawn separately for low / median / high values of each of the other four selected factors (the remaining factors held at their median). Fitting on the original (un-standardized) scale keeps the axes interpretable.

m = sm.GLM(y, sm.add_constant(d[best]), family=sm.families.Binomial()).fit()   # best model, ORIGINAL scale
b = m.params; med = {c: d[c].median() for c in best}; ages = np.linspace(d["age"].min(), d["age"].max(), 100)
def prob(age, **over):                                             # P(CHD) at a given age, other factors fixed
    v = {**med, **over, "age": age}
    return 1.0 / (1.0 + np.exp(-(b["const"] + sum(b[k]*v[k] for k in best))))
print("age effect: odds x%.2f per year;  P(CHD) rises %.2f -> %.2f from age 30 to 60 (median profile, no family history)"
      % (np.exp(b["age"]), prob(30, famhist=0), prob(60, famhist=0)))

# figure: P(CHD) vs age per other factor (low/median/high), with age-binned OBSERVED proportions
bins = np.linspace(d["age"].min(), d["age"].max(), 9); mids = 0.5 * (bins[:-1] + bins[1:])
def binned(mask):                                                  # observed CHD proportion + binomial SE per age bin
    a = d["age"].values[mask]; yy = d["chd"].values[mask]
    which = np.clip(np.digitize(a, bins) - 1, 0, len(mids) - 1)
    prop = np.full(len(mids), np.nan); se = np.full(len(mids), np.nan)
    for k in range(len(mids)):
        s = which == k
        if s.sum() > 0: pk = yy[s].mean(); prop[k] = pk; se[k] = np.sqrt(pk * (1 - pk) / s.sum())
    return prop, se

fig, axes = plt.subplots(2, 2, figsize=(11.5, 8.0), sharex=True, sharey=True); axes = axes.ravel()
# panel 1: family history (curves absent/present; observed points SPLIT by family history)
for val, col, lab in [(0, "navy", "no family history"), (1, "firebrick", "family history")]:
    axes[0].plot(ages, [prob(a, famhist=val) for a in ages], color=col, lw=2, label=lab)
    pr, se = binned(d["famhist"].values == val)
    axes[0].errorbar(mids, pr, yerr=se, fmt="o", mfc="none", mec=col, ecolor=col, ms=6, capsize=2, lw=1)
axes[0].set_title("varying famhist")
# panels 2-4: continuous factors (curves low/median/high; observed points POOLED)
pr, se = binned(np.ones(len(d), bool))
for ax, f in zip(axes[1:], ["tobacco", "ldl", "typea"]):
    for lab, val, col in [("low (Q1)", d[f].quantile(.25), "navy"),
                          ("median",   d[f].median(),      "gray"),
                          ("high (Q3)", d[f].quantile(.75), "firebrick")]:
        ax.plot(ages, [prob(a, **{f: val}) for a in ages], color=col, lw=2, label="%s = %.1f" % (f, val))
    ax.errorbar(mids, pr, yerr=se, fmt="o", mfc="none", mec="k", ecolor="k", ms=6, capsize=2, lw=1, label="observed")
    ax.set_title("varying %s" % f)
for ax in axes: ax.grid(alpha=.25); ax.legend(fontsize=8, loc="upper left")
axes[0].set_ylabel("P(CHD)"); axes[2].set_ylabel("P(CHD)")
axes[2].set_xlabel("age"); axes[3].set_xlabel("age")
plt.tight_layout(); plt.show()

age effect: odds x1.05 per year;  P(CHD) rises 0.11 -> 0.36 from age 30 to 60 (median profile, no family history)

Reading the curves

• The model tracks the data. The open circles are the observed CHD proportion in each age bin (with binomial standard-error bars) — split by family history in the first panel, pooled in the others. They fall along the fitted curves, so the smooth logistic is a fair summary of the empirical age-risk relationship, not an artefact of the functional form.
• Age is the backbone of risk. In every panel the probability climbs steadily with age — the odds rise about 5% per year, so a median-profile man's CHD probability roughly triples between 30 and 60.
• Each selected factor lifts the whole curve. A family history shifts the age curve up the most — at any age it roughly doubles the odds, and the split observed points confirm it (the family-history bins sit clearly above the no-history bins). Heavy tobacco, high LDL and high type-A each raise risk at every age too, but by less. Because these are additive on the log-odds scale, a risk factor shifts the age curve vertically rather than changing its shape.
• This is what a selected model buys you: interpretability. Five factors, each with a clear direction and magnitude, tell a coherent clinical story — age × heredity × smoking × cholesterol × behaviour — that the nine-factor model would only muddy with four near-zero, statistically indistinguishable coefficients.

Results — logistic selection, three ways

• The Polya-Gamma trick delivers. A non-conjugate logistic likelihood is turned into a sequence of Gaussian SSVS draws by one latent variable per observation — and the from-scratch Gibbs, PyMC's native logistic SSVS, and the exact enumeration all select the same five risk factors. The augmentation is what lets the entire Gaussian selection machinery of Part 1 carry over to binary data (and, unchanged, to counts, multinomial and beyond).
• Search $\approx$ enumeration, again. The two SSVS samplers give softer inclusion probabilities than the exact enumeration (a prior effect, as in Parts 1 vs 2), but the ranking and the selected set are identical. Where enumeration is possible it is the gold standard; where it is not ($2^p$ too large), the Polya-Gamma Gibbs scales to it.
• Selection is free here. The five-variable model's predicted probabilities track the full nine-variable model's almost exactly ($r=0.985$), and it generalizes slightly better out of sample (CV AUC 0.782 vs 0.770, lower log-loss) — the four dropped factors were contributing variance, not signal. That is the practical payoff of selection: a simpler, more interpretable model at no predictive cost — indeed a small gain.
• A clinically meaningful answer. Coronary risk here is driven by age, heredity, smoking, behaviour type and cholesterol — not by the body-size measures (adiposity, obesity), whose apparent link to CHD is explained away by the factors they travel with. That is exactly the distinction — conditional over marginal relevance — that Bayesian variable selection is built to draw.

References

• Polson, N. G., Scott, J. G. & Windle, J. (2013). Bayesian inference for logistic models using Polya-Gamma latent variables. JASA 108, 1339–1349.
• George, E. I. & McCulloch, R. E. (1993). Variable selection via Gibbs sampling. JASA 88, 881–889.
• Hastie, T., Tibshirani, R. & Friedman, J. (2009). The Elements of Statistical Learning (2nd ed.), Section 4.4. Springer.
• Rousseauw, J. et al. (1983). Coronary risk factor screening in three rural communities. South African Medical Journal 64, 430–436.

Data: saheart_data.csv — South African Heart Disease: 462 men, 9 risk factors, binary CHD (bestglm::SAheart). Next: logit_R.ipynb — the R cross-check (BAS Bayesian logistic model averaging).