Bayesian Logistic Regression

Definition

Bayesian logistic regression models a binary response $y_i \in \{0,1\}$ via $P(y_i=1) = (1+e^{-x_i\beta})^{-1}$ , placing a prior $\pi(\beta)$ over coefficients. Because the logistic likelihood is not conjugate to any standard prior for $\beta$ , 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: $\varepsilon_i \sim \text{Logistic}(0,1)$ is represented as $\varepsilon_i \mid \lambda_i \sim N(0, \lambda_i)$ with $\lambda_i = (2\psi_i)^2$ , $\psi_i \overset{\text{iid}}{\sim} \text{KS}$ (Kolmogorov-Smirnov distribution). This is exact — not the $t_8$ approximation obtained by $\lambda_i \sim \text{Gamma}(4,4)$ (Albert-Chib suggestion), which departs from the logistic in the tails.
Conditional conjugacy: Given latent $z$ and mixing variances $\lambda$ , the normal prior $\pi(\beta) = N(b, v)$ yields $\beta \mid z, \lambda \sim N(B, V)$ with $V = (v^{-1} + x'Wx)^{-1}$ , $B = V(v^{-1}b + x'Wz)$ , $W = \text{diag}(\lambda_i^{-1})$ . The logistic model is obtained from $V$ by updating $W$ at each iteration.
Two joint-update schemes:
- Scheme A (recommended): Update $\{z, \lambda\}$ jointly given $\beta$ ; $z_i \mid \beta, y_i \sim \text{Logistic}(x_i\beta, 1)$ truncated (sampled by inversion); then $\lambda_i \mid z_i, \beta$ by rejection with a generalized inverse Gaussian GIG $(0.5, 1, r_i^2)$ envelope where $r_i = z_i - x_i\beta$ .
- Scheme B: Update $\{z, \beta\}$ jointly given $\lambda$ (as in the probit joint update), then update $\lambda \mid z, \beta$ . Yields larger step distances but more autocorrelated samples.
Variable selection extension: Auxiliary variables allow joint updates of $(\gamma, \beta)$ by conditioning on $\{z, \lambda\}$ ; the MH acceptance reduces to a Bayes-factor form with $\beta$ 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 $\beta_{-j}$ — the conditional likelihood for $\beta_j$ becomes binary logistic with offset $C_{ij} = \log\sum_{k\neq j}\exp(x_i\beta_k)$ . Scheme A generalises by Gibbs loop over $j=1,\ldots,Q-1$ .
Rao-Blackwellization: The marginal posterior of any subset $\beta^*$ is estimated as the Gaussian mixture $\hat{p}(\beta^*\mid y) = M^{-1}\sum_{j=1}^M N(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

$y_i = \mathbf{1}[z_i > 0], \quad z_i = x_i\beta + \varepsilon_i, \quad \varepsilon_i \mid \lambda_i \sim N(0, \lambda_i), \quad \lambda_i = (2\psi_i)^2, \quad \psi_i \sim \text{KS}$

so $\varepsilon_i$ marginally follows $\text{Logistic}(0,1)$ .

Scheme A Algorithm (one iteration)

Draw $\beta \mid z, \lambda \sim N(B, V)$ with $V = (v^{-1}I + x'Wx)^{-1}$ , $B = Vx'Wz$ .
For each $i$ : draw $z_i \mid \beta, y_i$ from truncated logistic (by inversion).
For each $i$ : set $r_i = z_i - x_i\beta$ ; draw $\lambda_i \mid r_i$ using rejection sampler with GIG $(0.5, 1, r_i^2)$ proposal (A4 in Holmes-Held appendix). Acceptance $\alpha(\lambda) = \exp(0.5\lambda)\pi_{\text{KS}}(\lambda)$ using alternating-series squeezing for the KS density.

Rejection Sampler for $\lambda_i$ (A4)

Proposal $\lambda \sim \text{GIG}(0.5, 1, r^2)$ drawn via $r/\text{IG}(1,r)$ inversion; accept with $\alpha(\lambda) = \exp(0.5\lambda)\pi_{\text{KS}}(\lambda)$ where the KS density is evaluated through monotone alternating series with breakpoint $4/3$ . Acceptance rate $\approx 0.25$ for $r^2 \to 0$ , rising to $\approx 1$ for $r^2 > 10$ .

Overdispersion: Beta-Binomial Hierarchical Logit (Kahn-Raftery 1996)

Grouped binary data (successes out of $n_i$ trials per unit) often show extra-binomial variation — most visibly an excess of units with a zero rate — that plain logistic regression cannot accommodate. Kahn-Raftery (1996) generalize the logit to a beta-binomial hierarchical model: $y_i\sim\text{Bin}(n_i,p_i)$ with $p_i\sim\text{Beta}$ whose parameters depend on covariates, so each unit has its own latent rate around a covariate-driven mean. The posterior mean of $p_i$ is approximately a shrinkage estimator — a weighted average of the logistic-regression ensemble estimate and the unit's observed rate, weighting the ensemble more when $n_i$ is small (little own-information) — i.e. empirical-Bayes borrowing strength (see James-Stein / hierarchical shrinkage). Fully Bayesian inference over the hyperparameters is obtained by Laplace approximation with locally uniform priors (Raftery's approximate-Bayes-factor machinery for GLMs), which also yields a natural test for overdispersion. This is the pre-MCMC route to the same overdispersed-binary modeling that count-data regression handles for counts.

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 $t_8$ 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 $\lambda$ rejection sampling; $V$ recomputed at each $\lambda$ 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\times 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 $\lambda_i$ degrades as $r_i^2 \to 0$ ; Polson-Scott-Windle (2013) Pólya-Gamma scale mixture later provided an alternative with $\text{PG}(1,0)$ auxiliary variables and a simpler exact sampler.