Gamerman (1997) Sampling from the Posterior Distribution in Generalized Linear Mixed Models

Summary

Gamerman (1997) proposes a unified Markov chain Monte Carlo (MCMC) scheme for Bayesian generalized linear mixed models (GLMMs) by repurposing the iteratively weighted least squares (IWLS) linearization from classical GLMM fitting as a Metropolis proposal. At each iteration, the current linear predictor $\eta^{(t)}$ is used to compute pseudo-responses and weights, yielding a working Normal posterior $N(m^{(t)}, C^{(t)})$ that serves as the Metropolis-Hastings (MH) proposal for both fixed effects $\beta$ and subject-level random effects $\gamma_i$; variance components $\Sigma$ are drawn directly from a conjugate inverse-Wishart full conditional. The resulting Metropolis-within-Gibbs cycle achieves acceptance rates of 98% (binomial seeds example) and above 85% (Poisson epilepsy example) with no tuning, and extends readily to nested random effects, non-normal priors via scale mixtures, and unknown link functions.

Key Claims

• Weighted least squares (WLS) proposal is self-calibrating: Proposal mean $m^{(t)} = (R^{-1} + X'WX)^{-1}(R^{-1}a + X'W\tilde{y})$ and covariance $C^{(t)} = (R^{-1} + X'WX)^{-1}$ update automatically as $\eta^{(t)}$ evolves, adapting to the posterior curvature without manual tuning.
• Blocking structure: Three block types per iteration — (1) $\beta$ via Metropolis with WLS proposal; (2) each $\gamma_i$ separately via Metropolis with subject-level WLS proposal treating $X\beta$ as offset; (3) $\Sigma$ via conjugate inverse-Wishart direct draw. The multivariate blocking avoids slow mixing from high intra-block correlations.
• Seeds example (binomial): Crowder (1978) data, $n=21$ groups, logistic link; acceptance rate 98%; posterior mean $\hat\beta_0 \approx -0.55$, $\hat\beta_1 \approx 0.08$ (bacteria effect), $\hat\sigma^2 \approx 0.10$.
• Epilepsy example (Poisson, nested random effects): Thall-Vail (1990) data, 59 patients × 4 visits, log-linear link, two-level random effects: subject $\gamma_i \sim N(0,\sigma_1^2)$ and visit $\delta_{ij} \sim N(0,\sigma_2^2)$; estimated $\hat\sigma_1 \approx 0.50$, $\hat\sigma_2 \approx 0.36$ under Gaussian; acceptance above 85%.
• $t_4$ extension: Modeling $\delta_{ij}$ as a scale mixture of normals (a Student-$t$ prior with 4 degrees of freedom, $t_4$) reduces $\hat\sigma_2$ from 0.36 to 0.26, indicating the Gaussian model over-attributes outlier visits to variance rather than heavy tails.
• Non-normal priors via scale mixture: Student-t prior $\beta|\lambda \sim N(a, \lambda R)$ with $\lambda \sim \text{IG}(\nu/2, \nu/2)$ (inverse-gamma) adds one scalar Metropolis step for $\lambda$ with negligible overhead; unknown scale parameters handled by the same mechanism.
• Extensions: Unknown link functions, survival/frailty models (conjugate double generalized linear model, GLM), and hierarchical dynamic GLMs (DGLMs; Gamerman-Migon 1993) all fit within the same IWLS-proposal framework.

Concepts Introduced or Extended

• Generalized Linear Mixed Model — introduced; WLS proposal algorithm for full Bayesian inference
• Mixed-Effects Regression — extended to non-Gaussian exponential-family responses
• Metropolis-Hastings Algorithm — extended with WLS-derived working-model proposal
• Gibbs Sampler — extended; Metropolis-within-Gibbs blocking for GLMMs
• Bayesian Hierarchical Model — extended; nested random effects with non-conjugate likelihoods

Entities Mentioned

• Dani Gamerman

Quotes

"The main idea of the proposed sampling scheme is to use the normal approximation obtained from a first-order Taylor series expansion of the log-likelihood about the current values of the parameters as the candidate generating distribution."

"The proposal distribution for $\beta$ is a multivariate normal distribution centred at the WLS estimate, with variance proportional to the inverse of the information matrix of the working model."

My Take

The paper's core insight is elegant: the IWLS linearization used for decades in frequentist GLMM fitting already produces a locally accurate Normal approximation to the posterior — so why not use it as a Metropolis proposal? The resulting algorithm is almost parameter-free, adapts each iteration to the current state, and works across all exponential-family likelihoods without modification. This is a practical and theoretically sound shortcut that would have been practically unavailable without the Metropolis-Hastings formalism developed by Chib-Greenberg (1995).

Limitations worth noting: (1) The WLS proposal updates $\beta$ and $\gamma_i$ in separate blocks, ignoring the posterior correlation between them — for highly correlated fixed and random effects this may still lead to slow mixing. (2) The method assumes the working-model Normal approximation is adequate; for very sparse binary data (near-zero or near-one probabilities) the linearization can be poor. (3) The paper predates adaptive MCMC and Hamiltonian Monte Carlo / no-U-turn sampler (HMC/NUTS), which now dominate GLMM Bayesian computation in practice (Stan, integrated nested Laplace approximation, INLA). Nonetheless, the WLS-proposal concept directly influenced subsequent MCMC literature for non-conjugate hierarchical models, including Chib-Carlin (1999) for hierarchical longitudinal models.