Bayesian Random-Effects Panel Logit

Python · R  ·  Download Gibbs sampler

Model

A mixed (multilevel) logistic regression for binary outcomes measured repeatedly per subject. Each subject's repeated measures are correlated, so every subject gets a random intercept bjb_j shared across all of their visits — which is what models the within-subject correlation (distinct from a per-row random effect). It is the binary, longitudinal twin of the Epil Poisson GLMM: same subject-random-intercept structure, Bernoulli–logit likelihood. Unlike the Seeds binomial GLMM (random effect per plate), here the effect is shared across each subject's panel of observations. It is the binary-longitudinal member of a recurring random-effects GLMM family that runs across the catalog — Gaussian (Rats), Poisson (Epil), binomial (Seeds), and binary-panel (here) — one random-intercept idea, different likelihoods.

yit∼Bernoulli(pit),logit pit=xit′β+bs(i),bj∼N(0,σ2)y_{it} \sim \text{Bernoulli}(p_{it}), \qquad \text{logit}\,p_{it} = x_{it}'\beta + b_{s(i)}, \qquad b_j \sim N(0, \sigma^2)

Sampler — Pólya–Gamma data augmentation

This dataset is famous for a large, hard-to-estimate σ≈4\sigma \approx 4 (Lesaffre & Spiessens 2001) that defeats naive samplers — a from-scratch RW-Metropolis mixes too slowly and its estimates depend on the step size. The fix is Pólya–Gamma data augmentation (Polson, Scott & Windle 2013): augmenting latent ωit∼PG(1,ηit)\omega_{it} \sim \mathrm{PG}(1, \eta_{it}) makes the logit likelihood conditionally Gaussian, giving exact, tuning-free Gibbs steps for β\beta and the bjb_j (conjugate Normals) and an Inverse-Gamma for σ2\sigma^2 — the same trick as the binary logit and hierarchical logit examples. The PG draw here adds an analytic tail-mean correction so few series terms suffice even for near-separated subjects (large ∣η∣|\eta|).

BlockDrawFull conditional
(0)ωit∣⋅\omega_{it} \mid \cdotPólya–Gamma augmentation, PG(1,ηit)\mathrm{PG}(1, \eta_{it}) (tail-mean corrected)
(1)β∣⋅\beta \mid \cdotNormal conjugate, precision X′ΩX+B0−1X'\Omega X + B_0^{-1}
(2)bj∣⋅b_j \mid \cdotPer-subject Normal, precision ∑i∈jωi+1/σ2\sum_{i\in j}\omega_i + 1/\sigma^2; re-centred into β0\beta_0
(3)σ2∣⋅\sigma^2 \mid \cdotInverse-Gamma conjugate from the random effects

Notebooks

Applied to the toenail / onychomycosis trial (De Backer 1998): 294 patients randomised to two antifungals (itraconazole vs terbinafine), each evaluated at up to 7 visits over ~12 months (1908 observations), with a binary onycholysis outcome. Onycholysis falls steeply over the year in both arms (37% → 8% severe): the month coefficient is ≈−0.40\approx -0.40 per month (odds down ~33%/month), and terbinafine clears it faster — the treatment×time interaction is ≈−0.14\approx -0.14 (steeper decline), with the main treatment effect ≈0\approx 0 (arms equal at baseline by randomisation). The headline statistical feature is σ≈4\sigma \approx 4 — very large between-patient heterogeneity. Because σ\sigma is so large, the subject-specific (conditional) curve is much steeper than the population-averaged (marginal) one — the classic GLMM-vs-marginal attenuation, so these β\beta's are conditional within-patient effects, several times larger than a marginal/GEE model reports. Cross-checked three ways — from-scratch PG-Gibbs, PyMC (non-centered NUTS), and R lme4::glmer, which exposes the Lesaffre–Spiessens caveat that σ\sigma is numerically unstable at low quadrature — swinging from 4.56 to 3.69 as the Gauss–Hermite node count changes before settling at ≈4.0\approx 4.0 — because a large random-effect variance is genuinely hard to integrate with few nodes. The two exact MCMC samplers, which never approximate that integral, both land at σ≈4.1\sigma \approx 4.1.

Python Mixed logistic from scratch — Pólya–Gamma Gibbs, subject random intercept, onychomycosis trial (Python) Open → Python PyMC cross-check — non-centered Bernoulli GLMM, NUTS (Python) Open → R R cross-check — lme4::glmer with the quadrature-point (nAGQ) caveat (R) Open →

Downloads

longlogit_gibbs.py From-scratch mixed logistic — Pólya–Gamma Gibbs with tail-corrected PG draw (NumPy / scipy) toenail.csv Onychomycosis trial — 294 patients, up to 7 visits each (1908 obs)

Gibbs Sampler — Source Code

"""
Longitudinal / mixed logistic regression -- binary GLMM with a SUBJECT random intercept shared across each
subject's repeated measurements -- from-scratch Metropolis-within-Gibbs. The toenail / onychomycosis trial
(De Backer 1998; Lesaffre & Spiessens 2001): a binary outcome measured repeatedly per patient over time.

This is the BINARY analogue of the Poisson GLMM hpois_gibbs.py (Epil): same group-random-intercept
structure (one b per subject, applied to all of that subject's rows), Bernoulli-logistic likelihood
instead of Poisson. It is NOT seeds_gibbs (whose random effect is per-row) -- here the effect is shared
across the subject's visits, which is what models the within-patient correlation of repeated measures.

Model
-----
  y_it ~ Bernoulli(p_it),  logit p_it = x_it' beta + b_{s(i)},  b_j ~ N(0, sigma^2)
  x_it = [1, treat, time, treat*time]  (treat constant within patient; time varies across visits)

Sampler: POLYA-GAMMA data augmentation (Polson, Scott & Windle 2013) -> EXACT Gibbs, no tuning, mixes
well even when sigma is large (the toenail sigma ~4 is notoriously hard for RW/quadrature).
  augment  omega_it ~ PG(1, eta_it),  eta_it = x_it'beta + b_{s(i)};   kappa_it = y_it - 1/2
  1. beta   | omega, b   ~ Normal (precision X'Omega X + prior_prec)
  2. b_j    | omega, beta ~ Normal, per subject (precision sum_{i in j} omega_i + 1/sigma^2)  + centre -> beta_0
  3. sigma^2 | b          ~ Inverse-Gamma
(A from-scratch RW-Metropolis version mixes too slowly here -- the estimates depend on the step size.)
"""
import numpy as np
import pandas as pd
from scipy.special import expit


def _rpg1(c, rng, n_terms=100):
    """Sample PG(1, c) (vectorised) via the sum-of-Gammas representation with an analytic TAIL-MEAN
    correction: omega = (1/2pi^2) [ full_mean + sum_{k<=K} (g_k - 1)/denom_k ], where full_mean =
    (pi/2a) tanh(pi a), a = |c|/2pi. Replacing the truncated tail by its exact expectation removes the
    truncation bias that otherwise inflates the variance estimate for near-separated units (large |c|),
    so a small K suffices. g_k ~ Exp(1)."""
    c = np.abs(np.asarray(c, float)); a = c / (2 * np.pi)
    k = np.arange(1, n_terms + 1)
    denom = (k - 0.5) ** 2 + a[:, None] ** 2                     # (n, K)
    g = rng.standard_gamma(1.0, size=(len(c), n_terms))
    full = np.where(a > 1e-8, (np.pi / (2 * np.where(a > 1e-8, a, 1.0))) * np.tanh(np.pi * a), np.pi ** 2 / 2)
    fluct = np.sum((g - 1.0) / denom, axis=1)                    # mean-zero fluctuation of the leading terms
    return (1.0 / (2 * np.pi ** 2)) * (full + fluct)


def toenail_data(csv='toenail.csv'):
    d = pd.read_csv(csv)
    uniq = np.unique(d['patient']); idx = {p: i for i, p in enumerate(uniq)}
    group = d['patient'].map(idx).to_numpy()
    X = np.column_stack([np.ones(len(d)), d['treat'], d['month'], d['treat'] * d['month']])
    return d['y'].to_numpy(float), X, group


def simulate_longlogit(J=300, ni=7, beta=(-1.5, -0.2, -0.4, -0.15), sigma=2.5, seed=0):
    rng = np.random.default_rng(seed); beta = np.asarray(beta, float)
    treat = np.repeat(rng.integers(0, 2, J), ni).astype(float)
    time = np.tile(np.arange(ni), J).astype(float)
    group = np.repeat(np.arange(J), ni)
    X = np.column_stack([np.ones(J*ni), treat, time, treat*time])
    b = rng.normal(0, sigma, J)
    y = (rng.random(J*ni) < expit(X @ beta + b[group])).astype(float)
    return y, X, group, b


def _logit_mle(X, y, maxit=100):
    beta = np.zeros(X.shape[1])
    for _ in range(maxit):
        p = expit(X @ beta); w = p * (1 - p)
        step = np.linalg.solve((X * w[:, None]).T @ X + 1e-8*np.eye(X.shape[1]), X.T @ (y - p)); beta += step
        if np.max(np.abs(step)) < 1e-10:
            break
    return beta


def longlogit_gibbs(y, X, group, R=8000, burn=2000, prior_sd=10.0, a0=0.001, b0=0.001,
                    seed=0, n_terms=200):
    rng = np.random.default_rng(seed)
    y = np.asarray(y, float); n, k = X.shape; group = np.asarray(group, int); J = group.max() + 1
    kappa = y - 0.5; Xt = X.T; P0 = np.eye(k) / prior_sd ** 2
    beta = _logit_mle(X, y); b = np.zeros(J); sigma = 1.0
    keep = R - burn
    B = np.zeros((keep, k)); SIG = np.zeros(keep); Bm = np.zeros(J); BRE = np.zeros((keep, J))
    for it in range(R):
        eta = X @ beta + b[group]
        om = _rpg1(eta, rng, n_terms)                            # Polya-Gamma augmentation
        # 1. beta | omega, b  (Gaussian)
        prec = (X * om[:, None]).T @ X + P0
        cov = np.linalg.inv(prec); mb = cov @ (Xt @ (kappa - om * b[group]))
        beta = mb + np.linalg.cholesky(cov) @ rng.standard_normal(k)
        eta0 = X @ beta
        # 2. b_j | omega, beta  (per-subject Gaussian, vectorised)
        prec_b = np.bincount(group, weights=om, minlength=J) + 1.0 / sigma ** 2
        rhs_b = np.bincount(group, weights=(kappa - om * eta0), minlength=J)
        b = rhs_b / prec_b + rng.standard_normal(J) / np.sqrt(prec_b)
        m = b.mean(); b -= m; beta[0] += m                        # centre subject effects -> intercept
        # 3. sigma^2 | b  (Inverse-Gamma)
        sigma = np.sqrt(1.0 / rng.gamma(a0 + J / 2.0, 1.0 / (b0 + 0.5 * np.sum(b ** 2))))
        if it >= burn:
            kk = it - burn; B[kk] = beta; SIG[kk] = sigma; Bm += b; BRE[kk] = b
    return dict(beta=B, sigma=SIG, b=Bm / keep, bdraws=BRE)


def summary(draws, names=None):
    q = draws.shape[1]; names = names or [f'b{j}' for j in range(q)]
    return [(nm, draws[:, j].mean(), draws[:, j].std(), *np.percentile(draws[:, j], [2.5, 97.5]))
            for j, nm in enumerate(names)]

References