Bayesian Hierarchical Binary Logit

Python · R  ·  Download sampler module

Model

The hierarchical (multilevel) binary logit lets each group jj — here a U.S. state — have its own intercept αj\alpha_j and slope βj\beta_j, which are themselves regressed on a group-level predictor zj=(1,uj)z_j = (1, u_j) (state mean income). This decomposes a single income variable into an individual effect and a contextual effect — the structure behind the rich-state / poor-state voting paradox.

Pr⁡(yi=1)=logit−1 ⁣(αj[i]+βj[i] xi)\Pr(y_i = 1) = \text{logit}^{-1}\!\big(\alpha_{j[i]} + \beta_{j[i]}\, x_i\big)
(αjβj)∼N ⁣(Δ′zj, Vb)\binom{\alpha_j}{\beta_j} \sim N\!\big(\Delta' z_j,\ V_b\big)

The paradox. The within-state slope βj\beta_j is positive — richer individuals vote Republican more — while the between-state coefficient Δu→α\Delta_{u\to\alpha} is negative — richer states vote Republican less (Connecticut vs. Mississippi). A pooled single-level logit carries no state-level term, so it can report only one income coefficient: on these data it returns +0.17+0.17, the individual effect attenuated by the opposing contextual one. It does not flip sign — but the contextual effect is not merely mis-estimated, it is inexpressible. The multilevel model separates the individual income effect from the contextual one — the textbook resolution of Simpson's paradox (Gelman, Shor, Bafumi & Park 2007).

Sampler — hierarchical Pólya–Gamma Gibbs

Pólya–Gamma augmentation (Polson, Scott & Windle 2013) makes every level Gaussian, so the whole sampler is conjugate Gibbs — the multilevel extension of the flat-logit PG-Gibbs in the Bayesian Binary Logit example. Drawing latent ωi∼PG(1, ηi)\omega_i \sim \mathrm{PG}(1,\, \eta_i) renders the state effects, the group-level regression Δ\Delta, and the random-effect covariance VbV_b each a standard conjugate update.

BlockDrawFull conditional
(1) ωi∣⋅\omega_i \mid \cdotPG(1, αj+βjxi)\mathrm{PG}(1,\, \alpha_j + \beta_j x_i) — Pólya–Gamma augmentation, one per observation
(2) (αj,βj)∣⋅(\alpha_j,\beta_j) \mid \cdot N ⁣(Pj−1(Xj′κj+Vb−1Δ′zj), Pj−1)N\!\big(P_j^{-1}(X_j'\kappa_j + V_b^{-1}\Delta'z_j),\ P_j^{-1}\big), Pj=Xj′ΩjXj+Vb−1P_j = X_j'\Omega_j X_j + V_b^{-1}
(3) Δ∣⋅\Delta \mid \cdotMatrix-normal — group-level regression of the random effects on zjz_j
(4) Vb∣⋅V_b \mid \cdotInverse-Wishart on the random-effect residuals E=B−ZΔE = B - Z\Delta

with κ=y−12\kappa = y - \tfrac12 and Ωj=diag(ω)\Omega_j = \mathrm{diag}(\omega).

Notebooks

Python notebook. Section 1 validates the from-scratch sampler on synthetic 2-level data (50 states, n=5000n=5000): it recovers the known Δ\Delta and VbV_b and reproduces the opposite-sign within/between pattern (Δ01=+0.42\Delta_{01}=+0.42, Δ10=−0.86\Delta_{10}=-0.86). A sweep in the same section isolates what actually identifies a hierarchy: the group-level parameters Δ\Delta and VbV_b are learned from the JJ states, not the 5,000 individuals, so their precision scales with J\sqrt{J} — piling more respondents into each state leaves VbV_b stuck near 0.12, while adding states drives it to the true 0.10. Section 2 turns to the Gelman–Shor–Bafumi–Park 2004 exit poll (Harvard Dataverse; n=67,381n=67{,}381 across 50 states, Bush vs. Kerry): a positive within-state income slope (Δ01=+0.21\Delta_{01}=+0.21, all 50 state slopes positive) coexists with a negative contextual effect (Δ10=−0.31\Delta_{10}=-0.31, richer states lean Democratic) — the canonical Gelman picture, now with credible intervals. Section 3 refits the model in PyMC (LKJ covariance prior, non-centered NUTS) and lands within 0.002 of the hand-written Gibbs.

R notebook. bayesm::rhierBinLogit fits the same multilevel logit but samples the logit likelihood by random-walk Metropolis instead of Pólya–Gamma augmentation, inside the same conjugate hierarchy. Its posterior means match the Python fit to 0.001. Its credible intervals, however, come out 1.2–1.8× wider — and the notebook traces that to the prior, not to any disagreement: the two fits put different inverse-Wishart priors on VβV_\beta (from-scratch IW(ν=4, S=I)\mathrm{IW}(\nu=4,\ S=I), prior mean II; bayesm's default IW(ν=5, V=5I)\mathrm{IW}(\nu=5,\ V=5I), prior mean 2.5 I2.5\,I), and with only J=50J=50 states that prior still has real pull. Together: three samplers, two languages, one paradox.

Python Hierarchical Pólya–Gamma Gibbs, synthetic + Gelman 2004 data, PyMC cross-check (Python) Open → R bayesm::rhierBinLogit — RW-Metropolis hierarchical logit, R vs Python (R) Open →

Downloads

hier_logit_gibbs.py From-scratch hierarchical Pólya–Gamma Gibbs sampler (NumPy only) gelman_2004.csv 2004 exit poll — n=67,381 respondents (Bush/Kerry, income bracket, state) gelman_2004_states.csv State-level predictor — survey-derived mean income per state

Sampler Module — Source Code

"""
Bayesian hierarchical (multilevel) logistic regression — from-scratch Gibbs (Pólya–Gamma).

Model (varying intercept + varying slope, with a group-level predictor):
    y_i ~ Bernoulli(sigmoid(alpha_{j[i]} + beta_{j[i]} * x_i)),    i in group j[i]
    b_j = (alpha_j, beta_j)' ,   b_j ~ N(Delta' z_j, V_b)
where x_i is the individual covariate (e.g. income), z_j the group-level covariates
(e.g. [1, state_mean_income]), Delta (q x p) the group-level regression, V_b (p x p)
the random-effect covariance.

Pólya–Gamma augmentation (Polson–Scott–Windle 2013) makes every block conjugate:
  1. omega_i      ~ PG(1, eta_i)                              (data augmentation)
  2. b_j | .      ~ N( P_j^{-1}(X_j'kappa_j + Vb^{-1} Delta'z_j), P_j^{-1} ),
                    P_j = X_j' Omega_j X_j + Vb^{-1},  kappa = y - 1/2
  3. Delta | .    ~ MatrixNormal( (Z'Z+A0)^{-1} Z'B, (Z'Z+A0)^{-1}, V_b )
  4. V_b | .      ~ InverseWishart( nu0 + J, S0 + E'E ),  E = B - Z Delta

Reuses the PG sampler from binary_logit_mcmc. NumPy only. Companion to
binary_logit_mcmc.py (flat logit) and mts_gibbs.py.
"""
import numpy as np
from binary_logit_mcmc import rpg, sigmoid


def simulate_hier_logit(J=50, n_per=40, Delta=None, Vb=None, seed=0):
    """Simulate a 2-level varying-intercept/slope logit that exhibits the
    'rich state / poor state' pattern: positive within-group slope, negative
    group-level intercept–predictor relationship."""
    rng = np.random.default_rng(seed)
    u = rng.normal(size=J)                          # group-level predictor (e.g. state mean income)
    Z = np.column_stack([np.ones(J), u])            # (J, q=2)
    if Delta is None:
        # rows = group covariates (intercept, u); cols = REs (alpha, beta)
        Delta = np.array([[-0.10, 0.40],            # baseline alpha,   baseline beta (+ income slope)
                          [-0.70, 0.00]])           # effect of u on alpha (NEGATIVE), on beta (~0)
    if Vb is None:
        Vb = np.array([[0.30, 0.05], [0.05, 0.10]])
    B = Z @ Delta + rng.normal(size=(J, 2)) @ np.linalg.cholesky(Vb).T   # (J, 2) true REs
    state  = np.repeat(np.arange(J), n_per)
    income = rng.normal(size=J * n_per)             # individual predictor (centered)
    eta = B[state, 0] + B[state, 1] * income
    y = (rng.random(len(eta)) < sigmoid(eta)).astype(float)
    return dict(y=y, income=income, state=state, u=u, Z=Z,
                B_true=B, Delta_true=Delta, Vb_true=Vb)


def _riwishart(nu, S, rng):
    """Draw V ~ Inverse-Wishart(nu, scale=S) via a Bartlett-decomposition Wishart."""
    p = S.shape[0]
    L = np.linalg.cholesky(np.linalg.inv(S))        # chol of S^{-1}
    A = np.zeros((p, p))
    for i in range(p):
        A[i, i] = np.sqrt(rng.chisquare(nu - i))
        for k in range(i):
            A[i, k] = rng.standard_normal()
    W = L @ A @ A.T @ L.T                            # Wishart(nu, S^{-1})
    return np.linalg.inv(W)                           # Inverse-Wishart(nu, S)


def hier_logit_gibbs(y, income, state, Z, R=4000, burn=1000, seed=0,
                     A0=0.01, nu0=None, S0=None, n_terms=128):
    """Hierarchical PG-Gibbs. Returns posterior draws of B (J×2), Delta (q×2), V_b (2×2)."""
    rng = np.random.default_rng(seed)
    y = np.asarray(y, float); income = np.asarray(income, float); state = np.asarray(state, int)
    J, q = Z.shape; p = 2
    X = np.column_stack([np.ones(len(y)), income])   # individual design (1, income)
    kappa = y - 0.5
    idx = [np.where(state == j)[0] for j in range(J)]
    A0 = A0 * np.eye(q)
    nu0 = p + 2 if nu0 is None else nu0
    S0 = np.eye(p) if S0 is None else S0
    ZtZ = Z.T @ Z

    B = np.zeros((J, p)); Delta = np.zeros((q, p)); Vb = np.eye(p)
    keep = R - burn
    Bd = np.zeros((keep, J, p)); Dd = np.zeros((keep, q, p)); Vd = np.zeros((keep, p, p))

    for g in range(R):
        # 1. Polya-Gamma augmentation
        eta = (X * B[state]).sum(1)
        omega = rpg(1.0, eta, rng, n_terms)
        # 2. random effects b_j (Gaussian full conditional, per group)
        Vbinv = np.linalg.inv(Vb)
        for j in range(J):
            ij = idx[j]; Xj = X[ij]; oj = omega[ij]
            Pj = (Xj * oj[:, None]).T @ Xj + Vbinv
            rhs = Xj.T @ kappa[ij] + Vbinv @ (Delta.T @ Z[j])
            cov = np.linalg.inv(Pj)
            B[j] = cov @ rhs + np.linalg.cholesky(cov) @ rng.standard_normal(p)
        # 3. group-level regression Delta (matrix-normal full conditional)
        U = np.linalg.inv(ZtZ + A0)
        M = U @ (Z.T @ B)
        Delta = M + np.linalg.cholesky(U) @ rng.standard_normal((q, p)) @ np.linalg.cholesky(Vb).T
        # 4. RE covariance V_b (inverse-Wishart full conditional)
        E = B - Z @ Delta
        Vb = _riwishart(nu0 + J, S0 + E.T @ E, rng)
        if g >= burn:
            k = g - burn; Bd[k] = B; Dd[k] = Delta; Vd[k] = Vb
    return dict(B=Bd, Delta=Dd, Vb=Vd)


def summary(out):
    return dict(Delta=out['Delta'].mean(0), Vb=out['Vb'].mean(0), B=out['B'].mean(0))

References