Pólya-Gamma Sampling for Logistic Random-Effects Models

Why Pólya-Gamma Augmentation Helps

In a logistic random-effects model, the log likelihood is not quadratic in the coefficients. A Gaussian prior therefore does not produce a Gaussian posterior conditional on the observed outcomes alone. Pólya-Gamma augmentation introduces latent positive variables that make the likelihood quadratic in the linear predictors. Conditional on those variables, the fixed effects and random effects have a multivariate Gaussian distribution and can be updated jointly.

This is a conditional simplification, not a claim that the marginal posterior is Gaussian. The latent variables must be sampled, and variance and correlation parameters generally still require their own updates. The construction applies to the usual logistic link and is particularly useful when Gaussian priors are specified for regression coefficients and random effects.

Model and the Pólya-Gamma Identity

Suppose there are \(m\) conditionally independent binomial observations. For observation \(i\), \(y_i\) of \(n_i\) trials are successes, and the success probability is \(p_i\). Conditional on the model parameters,

Here \(o_i\) is a known offset, \(\mathbf c_i^\mathsf{T}\) is row \(i\) of a design matrix \(\mathbf C\), and \(\boldsymbol\theta\) collects the fixed-effect coefficients and random effects. The binomial observations are independent conditional on these quantities. For a GLMM, the random effects have a Gaussian prior conditional on variance and correlation parameters. In the derivations below, those parameters are held fixed, and the resulting prior for \(\boldsymbol\theta\) is written \(N(\mathbf m,\mathbf Q^{-1})\), with positive-definite precision matrix \(\mathbf Q\).

Set \(\kappa_i=y_i-n_i/2\). The binomial likelihood for observation \(i\) is

The Pólya-Gamma distribution \(\operatorname{PG}(b,c)\) is a distribution on positive values. The identity needed here is its Laplace transform:

Since \(1+\exp(\psi)=2\exp(\psi/2)\cosh(\psi/2)\), the logistic likelihood can be rewritten using this transform. The full binomial likelihood, not just its kernel, satisfies

The conditional distribution follows because the joint density of \(\omega_i\) and the data, as a function of \(\omega_i\) given \(\eta_i\), is proportional to \(\exp(-\omega_i\eta_i^2/2)p(\omega_i\mid n_i,0)\), the exponential tilting that defines \(\operatorname{PG}(n_i,\eta_i)\). Thus the augmentation has two alternating steps: draw each \(\omega_i\) given the current \(\eta_i\), then draw the coefficient vector given the \(\omega_i\)'s.

Deriving the Gaussian Conditional Update

Write \(\boldsymbol\kappa=(\kappa_1,\ldots,\kappa_m)^\mathsf{T}\), \(\boldsymbol\omega=(\omega_1,\ldots,\omega_m)^\mathsf{T}\), and \(\boldsymbol\Omega=\operatorname{diag}(\omega_1,\ldots,\omega_m)\). Conditional on the latent variables, the part of the joint likelihood that depends on \(\boldsymbol\theta\) is

Multiplying by the Gaussian prior density adds the exponent \(-(\boldsymbol\theta-\mathbf m)^\mathsf{T}\mathbf Q(\boldsymbol\theta-\mathbf m)/2\). Collecting the terms that involve \(\boldsymbol\theta\) gives a quadratic form with precision and precision-weighted mean

Completing the square therefore gives the conditional distribution

This derivation requires a Gaussian prior with a positive-definite precision matrix conditional on the variance and correlation parameters, and the logistic binomial likelihood specified above. Because every \(\omega_i>0\), the matrix \(\mathbf C^\mathsf{T}\boldsymbol\Omega\mathbf C\) is positive semidefinite. Adding it to positive-definite \(\mathbf Q\) leaves \(\mathbf Q_\ast\) positive definite, so the Gaussian conditional is well-defined even if the design matrix alone is rank deficient.

Including Random Effects

The Gaussian update applies to fixed and random effects together, provided their joint prior is Gaussian conditional on the variance and correlation parameters. For example, let \(\boldsymbol\theta=(\boldsymbol\beta^\mathsf{T},\mathbf b^\mathsf{T})^\mathsf{T}\), where \(\mathbf b\) contains group-level effects. The design matrix \(\mathbf C\) then contains both fixed-effect and random-effect columns. Its cross-products automatically include the dependence induced by observations sharing a group. A random effect is not independently updated observation by observation; its posterior information accumulates across every observation to which it contributes.

What the Method Does—and Does Not—Solve

Pólya-Gamma sampling is an exact augmentation: integrating out each \(\omega_i\) recovers the original logistic likelihood. Conditional Gaussianity makes the coefficient update convenient, but it does not eliminate all computational work. The \(\omega_i\)'s must be generated from their tilted distributions, and random-effect variance and correlation parameters do not generally become conjugate merely because the logistic likelihood has been augmented. Their updates must still respect their specified priors and the model's hierarchical structure.

Compared with the latent-variable augmentation for binary GLMMs, which commonly represents binary outcomes through thresholded latent variables and a probit link, Pólya-Gamma augmentation preserves the logistic link directly. It also differs from a random-walk Metropolis update for coefficients: here the coefficient block has an exact Gaussian conditional, rather than a proposal accepted or rejected using a likelihood ratio. Neither comparison means that the augmentation is universally faster; performance depends on model dimension, data, and the quality of the latent-variable sampler.

Several practical checks help identify errors. Verify that \(\kappa_i=y_i-n_i/2\), that the Pólya-Gamma shape is \(n_i\) for a binomial observation, and that the offset term enters as \(-\mathbf C^\mathsf{T}\boldsymbol\Omega\mathbf o\) in \(\mathbf h_\ast\). Check that \(\mathbf Q_\ast\) is symmetric positive definite and that simulated \(\omega_i\)'s are positive. Comparing a small model's likelihood-based posterior calculations with the augmented calculation can expose missing constants or incorrect signs.

A common mistake is to treat the identity as an exact representation of the likelihood kernel without its multiplier, or to say that the likelihood is exactly \(\exp(\kappa_i\eta_i-\omega_i\eta_i^2/2)\) mixed over \(\operatorname{PG}(n_i,0)\). The correct mixture for the kernel includes \(2^{-n_i}\); the full likelihood additionally includes \(\binom{n_i}{y_i}\). Another mistake is to call the marginal posterior Gaussian. Only the conditional distribution of the coefficient vector, given the latent variables and variance parameters, is Gaussian. The marginal distribution averages over those latent variables and is generally not Gaussian.