Definition
Bayesian logistic regression models a binary response yi∈{0,1} via P(yi=1)=(1+e−xiβ)−1, placing a prior π(β) over coefficients. Because the logistic likelihood is not conjugate to any standard prior for β, posterior inference requires Markov Chain Monte Carlo (MCMC). The key algorithmic insight of Holmes-Held (2006) is that the logistic error distribution admits an exact scale mixture of normals representation via the Kolmogorov-Smirnov (KS) distribution, enabling a fully automatic Gibbs sampler with no Metropolis-Hastings (MH) accept/reject steps.
Key Ideas
- Scale mixture of normals for logistic errors: εi∼Logistic(0,1) is represented as εi∣λi∼N(0,λi) with λi=(2ψi)2, ψi∼iidKS (Kolmogorov-Smirnov distribution). This is exact — not the t8 approximation obtained by λi∼Gamma(4,4) (Albert-Chib suggestion), which departs from the logistic in the tails.
- Conditional conjugacy: Given latent z and mixing variances λ, the normal prior π(β)=N(b,v) yields β∣z,λ∼N(B,V) with V=(v−1+x′Wx)−1, B=V(v−1b+x′Wz), W=diag(λi−1). The logistic model is obtained from V by updating W at each iteration.
- Two joint-update schemes:
- Scheme A (recommended): Update {z,λ} jointly given β; zi∣β,yi∼Logistic(xiβ,1) truncated (sampled by inversion); then λi∣zi,β by rejection with a generalized inverse Gaussian GIG(0.5,1,ri2) envelope where ri=zi−xiβ.
- Scheme B: Update {z,β} jointly given λ (as in the probit joint update), then update λ∣z,β. Yields larger step distances but more autocorrelated samples.
- Variable selection extension: Auxiliary variables allow joint updates of (γ,β) by conditioning on {z,λ}; the MH acceptance reduces to a Bayes-factor form with β marginalized out (equation 12). Efficient dimension-jumping; Pima Indian example: glucose and body mass index (BMI) >99% posterior inclusion, blood pressure <1%.
- Polychotomous regression: For Q unordered classes, condition on β−j — the conditional likelihood for βj becomes binary logistic with offset Cij=log∑k=jexp(xiβk). Scheme A generalises by Gibbs loop over j=1,…,Q−1.
- Rao-Blackwellization: The marginal posterior of any subset β∗ is estimated as the Gaussian mixture p^(β∗∣y)=M−1∑j=1MN(B∗(j),V∗(j)) — the normalising constant is always known within the auxiliary variable framework, enabling exact simultaneous credible-region computation (Held 2004).
How It Works
State-Space Form
yi=1[zi>0],zi=xiβ+εi,εi∣λi∼N(0,λi),λi=(2ψi)2,ψi∼KS
so εi marginally follows Logistic(0,1).
Scheme A Algorithm (one iteration)
- Draw β∣z,λ∼N(B,V) with V=(v−1I+x′Wx)−1, B=Vx′Wz.
- For each i: draw zi∣β,yi from truncated logistic (by inversion).
- For each i: set ri=zi−xiβ; draw λi∣ri using rejection sampler with GIG(0.5,1,ri2) proposal (A4 in Holmes-Held appendix). Acceptance α(λ)=exp(0.5λ)πKS(λ) using alternating-series squeezing for the KS density.
Rejection Sampler for λi (A4)
Proposal λ∼GIG(0.5,1,r2) drawn via r/IG(1,r) inversion; accept with α(λ)=exp(0.5λ)πKS(λ) where the KS density is evaluated through monotone alternating series with breakpoint 4/3. Acceptance rate ≈0.25 for r2→0, rising to ≈1 for r2>10.
Why It Matters
- Provides a fully automatic, exact Gibbs sampler for Bayesian logistic regression — no accept/reject steps, no user-set tuning, embeddable in generic software.
- The KS scale mixture is exact; the t8 approximation (Albert-Chib suggestion) introduces quantifiable tail error.
- The variable-selection extension via auxiliary conditioning turns Reversible Jump MCMC (RJMCMC) for logistic regression into a clean Bayes-factor MH step, avoiding the usual instability of dimension-changing moves.
- Rao-Blackwellization yields analytic normalising constants, enabling simultaneous credible intervals without quadrature or numerical integration.
Open Questions
- Computational cost vs. probit: Each logistic iteration is slower (per-observation λ rejection sampling; V recomputed at each λ draw). Whether Scheme A or B is better depends on the language and hardware; the recommendation may not hold in compiled code.
- High-dimensional p: The p×p matrix V must be computed at each iteration; for large p or sparse x, alternative approaches (stochastic gradient Langevin, No-U-Turn Sampler [NUTS]) are competitive.
- GIG sampler efficiency: The alternating-series rejection sampler for λi degrades as ri2→0; Polson-Scott-Windle (2013) Pólya-Gamma scale mixture later provided an alternative with PG(1,0) auxiliary variables and a simpler exact sampler.
Related