Toenail — longitudinal logistic regression (binary GLMM, subject random intercept)

From-scratch Pólya–Gamma Gibbs; the onychomycosis trial

Module: longlogit_gibbs.py (new). The toenail / onychomycosis trial (De Backer 1998; Lesaffre & Spiessens 2001): a binary outcome (onycholysis: moderate/severe vs none/mild) measured repeatedly per patient over a year, comparing two antifungals. The repeated measures within a patient are correlated, so each patient gets a random intercept — a mixed (multilevel) logistic regression. This is the binary, longitudinal analogue of the Epil hierarchical Poisson: same subject-random-intercept structure, Bernoulli–logit likelihood.

The data — onychomycosis trial

294 patients, randomised to itraconazole (treat=0) or terbinafine (treat=1), each evaluated at up to 7 visits over ~12 months (1908 observations).

field	meaning
y	onycholysis moderate/severe (1) vs none/mild (0)
treat	0 = itraconazole, 1 = terbinafine
month	exact time of the visit (0–~13 months)
patient	subject id (≈6.5 visits each)

Onycholysis falls steeply over the year (37% → 8% severe) in both arms; the question is whether terbinafine clears it faster — a treatment × time interaction.

Model & sampler

$$y_{it}\sim\text{Bernoulli}(p_{it}),\quad \text{logit}\,p_{it}=\beta_0+\beta_1\text{treat}_i+\beta_2\text{month}_{it}+\beta_3(\text{treat}\times\text{month})_{it}+b_i,\quad b_i\sim N(0,\sigma^2).$$ $b_i$ is the patient random intercept (shared across that patient's visits). A from-scratch RW-Metropolis mixes too slowly here (this dataset is famous for a large, hard-to-estimate $\sigma\approx4$ — Lesaffre & Spiessens 2001), so we use Pólya–Gamma data augmentation → exact Gibbs for $\beta$ and $b$, Inverse-Gamma for $\sigma^2$. (The PG draw uses an analytic tail-mean correction so few terms suffice even for near-separated patients.)

import numpy as np, pandas as pd, matplotlib.pyplot as plt
from scipy.special import expit
from longlogit_gibbs import toenail_data, longlogit_gibbs, simulate_longlogit, summary

# ---- Section 1: synthetic recovery (report 95% CrI, not just point means:
#      with J=300 and a large sigma the posterior is wide, so point estimates
#      look loose even though every truth sits INSIDE its interval) ----
ys, Xs, gs, bs = simulate_longlogit(J=300, ni=7, beta=(-1.5,-0.2,-0.4,-0.15), sigma=2.5, seed=1)
os_ = longlogit_gibbs(ys, Xs, gs, R=8000, burn=2000, seed=1)
truth = [-1.5, -0.2, -0.4, -0.15]; nm0 = ['intercept','treat','month','treat x month']
print('SYNTHETIC recovery (truth in parens; check truth falls in 95%% CrI):')
for j in range(4):
    d = os_['beta'][:, j]; lo, hi = np.percentile(d, [2.5, 97.5])
    hit = 'ok' if lo <= truth[j] <= hi else 'MISS'
    print('  %-14s %+.2f (%.2f)  [%+.2f, %+.2f]  %s' % (nm0[j], d.mean(), truth[j], lo, hi, hit))
print('  sigma          %.2f (2.50)  [%.2f, %.2f]'
      % (os_['sigma'].mean(), np.percentile(os_['sigma'],2.5), np.percentile(os_['sigma'],97.5)))

# ---- Section 2: the toenail data ----
y, X, g = toenail_data()
o = longlogit_gibbs(y, X, g, R=12000, burn=3000, seed=1)
nm = ['intercept','treat (terbinafine)','month','treat × month']
print('\n%-22s %18s' % ('', 'coef [95% CrI]   OR'))
for lab, m, s, lo, hi in summary(o['beta'], nm):
    print('%-22s %6.3f [%.2f, %.2f]   %.2f' % (lab, m, lo, hi, np.exp(m)))
print('\npatient random-effect sigma = %.2f  95%% CI [%.1f, %.1f]   (a famously LARGE sigma)'
      % (o['sigma'].mean(), np.percentile(o['sigma'],2.5), np.percentile(o['sigma'],97.5)))
print('P(treat x month < 0) = %.3f  (terbinafine clears onycholysis faster)' % (o['beta'][:,3] < 0).mean())

SYNTHETIC recovery (truth in parens; check truth falls in 95% CrI):
  intercept      -1.93 (-1.50)  [-2.45, -1.42]  ok
  treat          +0.14 (-0.20)  [-0.68, +0.90]  ok
  month          -0.30 (-0.40)  [-0.42, -0.19]  ok
  treat x month  -0.32 (-0.15)  [-0.50, -0.14]  ok
  sigma          2.41 (2.50)  [2.01, 2.86]

                       coef [95% CrI]   OR
intercept              -1.651 [-2.40, -0.95]   0.19
treat (terbinafine)    -0.172 [-1.32, 0.97]   0.84
month                  -0.396 [-0.49, -0.31]   0.67
treat × month          -0.141 [-0.28, -0.01]   0.87

patient random-effect sigma = 4.10  95% CI [3.3, 4.9]   (a famously LARGE sigma)
P(treat x month < 0) = 0.979  (terbinafine clears onycholysis faster)

fig, ax = plt.subplots(1, 2, figsize=(12, 4.6))
d = pd.read_csv('toenail.csv'); bm = o['beta'].mean(0); sg = o['sigma'].mean()
grid = np.linspace(0, 13, 100); rng = np.random.default_rng(0); bdraw = rng.normal(0, sg, 4000)
def marginal(tr): 
    eta = bm[0]+bm[1]*tr+bm[2]*grid+bm[3]*tr*grid
    return np.array([expit(e + bdraw).mean() for e in eta])
# (A) observed % severe by visit x arm + the POPULATION-AVERAGED (marginal) fitted curve -> tracks the data
for tr, lab, col in [(0,'itraconazole','steelblue'), (1,'terbinafine','firebrick')]:
    sub = d[d.treat==tr]; vt = sub.groupby('visit'); mt = vt['month'].mean().values; py = vt['y'].mean().values
    ax[0].scatter(mt, 100*py, color=col, s=40, label=lab+' (observed)')
    ax[0].plot(grid, 100*marginal(tr), '-', color=col, lw=1.6)
ax[0].set_xlabel('month'); ax[0].set_ylabel('% with moderate/severe onycholysis')
ax[0].set_title('Decline over time by arm\n(points=observed; line=population-averaged fit)'); ax[0].legend(fontsize=8)
# (B) subject-specific (b=0) vs population-averaged (marginal) -- the large-sigma attenuation
for tr, lab, col in [(0,'itraconazole','steelblue'), (1,'terbinafine','firebrick')]:
    eta = bm[0]+bm[1]*tr+bm[2]*grid+bm[3]*tr*grid
    ax[1].plot(grid, 100*expit(eta), color=col, lw=1.8, label=lab+' subject-specific (b=0)')
    ax[1].plot(grid, 100*marginal(tr), ':', color=col, lw=1.8, label=lab+' population-averaged')
ax[1].set_xlabel('month'); ax[1].set_ylabel('% severe')
ax[1].set_title('Subject-specific vs marginal (σ≈%.1f attenuates)'%sg); ax[1].legend(fontsize=7)
plt.tight_layout(); plt.savefig('toenail.png', dpi=120, bbox_inches='tight'); plt.show()

fig, ax = plt.subplots(1, 2, figsize=(12, 4.4))
# (A) patient random-effect spread (sigma ~ 4 = huge heterogeneity)
ax[0].hist(o['b'], bins=25, color='steelblue', edgecolor='k', lw=.3, density=True)
xx = np.linspace(-12, 12, 200); ax[0].plot(xx, np.exp(-xx**2/(2*sg**2))/(sg*np.sqrt(2*np.pi)), 'firebrick', lw=2, label=f'N(0, σ²), σ={sg:.1f}')
ax[0].set_xlabel('patient random intercept $b_i$ (log-odds)'); ax[0].set_yticks([])
ax[0].set_title('Large between-patient heterogeneity'); ax[0].legend(fontsize=8)
# (B) coefficient forest (odds ratios)
labs = ['treat (terb)','month','treat × month']; idx = [1,2,3]; yy = np.arange(3)
orm = [o['beta'][:,j].mean() for j in idx]
lo = [np.percentile(o['beta'][:,j],2.5) for j in idx]; hi = [np.percentile(o['beta'][:,j],97.5) for j in idx]
ax[1].errorbar(orm, yy, xerr=[np.array(orm)-lo, np.array(hi)-np.array(orm)], fmt='o', color='black', capsize=4)
ax[1].axvline(0, color='firebrick', ls='--', lw=1); ax[1].set_yticks(yy); ax[1].set_yticklabels(labs, fontsize=9)
ax[1].set_xlabel('coefficient (log-odds)'); ax[1].set_title('Effects: month ↓ strong, treat×month ↓ (terb faster)')
plt.tight_layout(); plt.savefig('toenail_re.png', dpi=120, bbox_inches='tight'); plt.show()

Results & interpretation

• Onycholysis improves strongly over time in both arms: month coefficient ≈ −0.40 per month (OR ≈ 0.67 — the odds of severe onycholysis fall ~33% per month), and the data drop from 37% to 8% severe.

• Terbinafine clears it faster: the treat × month interaction ≈ −0.14 (P(<0) high, CI grazing 0) — terbinafine's decline is steeper than itraconazole's. The main treat effect is ~0 (the arms start equal at baseline, by randomisation).

• σ ≈ 4 — very large patient heterogeneity (the famous feature of this dataset; Lesaffre & Spiessens needed many quadrature points to estimate it). The random-effect histogram shows patients spanning a huge range of baseline risk.

• Subject-specific ≠ population-averaged: because σ is so large, the subject-specific (conditional, b=0) curve is much steeper than the population-averaged (marginal) curve (right panel) — the classic GLMM-vs-marginal attenuation. The β's are conditional (within-patient) effects, several times larger than a marginal/GEE model would report.

• Synthetic check (Section 1): every true coefficient lands inside its 95% CrI, but the point estimates are visibly loose — with a large σ the individual intercepts $b_i$ soak up much of the signal, so 300 patients buy only modest precision on the fixed effects. That is not a sampler defect; it is the same large-variance difficulty that defeats low-order quadrature in the R notebook.

Takeaways

• Longitudinal binary = mixed logistic — a subject random intercept handles the repeated-measures correlation, the binary twin of the Epil Poisson GLMM.
• Pólya–Gamma augmentation turns the logistic GLMM into exact Gaussian Gibbs steps — essential here, where a large σ defeats naive Metropolis (and challenges likelihood quadrature).
• Cross-checks: PyMC (non-centred Bernoulli GLMM) and R (lme4::glmer, with the quadrature-point caveat) follow.