Bayesian inference for logistic models using Polya-Gamma latent variables

We propose a new data-augmentation strategy for fully Bayesian inference in models with binomial likelihoods. The approach appeals to a new class of Polya-Gamma distributions, which are constructed in detail. A variety of examples are presented to show the versatility of the method, including logistic regression, negative binomial regression, nonlinear mixed-effects models, and spatial models for count data. In each case, our data-augmentation strategy leads to simple, effective methods for posterior inference that: (1) circumvent the need for analytic approximations, numerical integration, or Metropolis-Hastings; and (2) outperform other known data-augmentation strategies, both in ease of use and in computational efficiency. All methods, including an efficient sampler for the Polya-Gamma distribution, are implemented in the R package BayesLogit. In the technical supplement appended to the end of the paper, we provide further details regarding the generation of Polya-Gamma random variables; the empirical benchmarks reported in the main manuscript; and the extension of the basic data-augmentation framework to contingency tables and multinomial outcomes.

💡 Research Summary

The paper introduces a novel data‑augmentation scheme for fully Bayesian inference in models with binomial (and multinomial) likelihoods, based on a newly defined class of Polya‑Gamma (PG) distributions. The key insight is that the logistic likelihood can be expressed as a mixture of Gaussian kernels with respect to a PG latent variable, thereby converting the non‑conjugate logistic link into a conditionally conjugate Gaussian form. This transformation enables a simple Gibbs sampler that alternates between drawing the regression coefficients (or other linear predictors) from a multivariate normal distribution and sampling the PG latent variables from their exact PG(b, c) distribution.

The authors first construct the PG distribution, derive its moment‑generating function, and provide an efficient exact sampler using a series representation combined with an accept‑reject step. This sampler is implemented in the R package BayesLogit, ensuring practical usability.

Four model families are examined: (1) standard logistic regression, where each observation contributes a PG(1, ψi) term; (2) negative‑binomial regression for over‑dispersed count data, achieved by augmenting the NB likelihood with PG variables; (3) nonlinear mixed‑effects models, where random effects are incorporated into the linear predictor and still retain conditional Gaussianity; and (4) spatial count models that embed a Gaussian random field in the predictor, allowing spatial correlation to be handled without additional approximations. In each case, the PG augmentation eliminates the need for Metropolis‑Hastings steps, Laplace approximations, or variational methods.

Extensive empirical benchmarks compare the PG‑augmented Gibbs sampler against Hamiltonian Monte Carlo, Metropolis‑Adjusted Langevin Algorithms, and earlier data‑augmentation approaches such as Albert‑Chib and scale‑mixture schemes. Performance metrics include effective sample size per second, convergence diagnostics, and total runtime. The PG method consistently achieves higher effective sample sizes and faster convergence, especially in high‑dimensional settings or when predictors are highly collinear.

The technical supplement extends the framework to multinomial outcomes and contingency tables. By assigning independent PG variables to each category (or cell), the same conjugate Gaussian updating step applies, preserving the simplicity of the Gibbs loop for multi‑category problems.

Overall, the paper provides a rigorous theoretical foundation for the Polya‑Gamma distribution, a practical algorithm for its exact simulation, and demonstrates that the resulting data‑augmentation strategy yields a universally applicable, computationally efficient, and easy‑to‑implement solution for Bayesian inference in a broad class of logistic‑type models.