Generalized Linear Mixed Model

Definition

A generalized linear mixed model (GLMM) extends the generalized linear model (GLM) to clustered or longitudinal data by incorporating subject-specific random effects. For subject $i$ with observations $y_{ij}$ , the conditional mean $\mu_{ij} = E[y_{ij}|\gamma_i]$ satisfies $g(\mu_{ij}) = x_{ij}'\beta + z_{ij}'\gamma_i$ , where $g(\cdot)$ is a link function, $\beta$ are fixed effects common to all subjects, and $\gamma_i \sim N(0,\Sigma)$ are subject-specific random deviations. The response $y_{ij}|\gamma_i$ follows an exponential-family distribution with mean $\mu_{ij}$ , accommodating binary, count, multinomial, and survival outcomes.

Key Ideas

Link function and exponential family: The mean-linear predictor relationship is mediated by $g(\mu_{ij})$ ; common choices are logit (binary), log (Poisson/count), and probit (binary). The variance function $b''(\theta_{ij})$ determines how spread scales with the mean.
Marginal intractability: Unlike the normal Mixed Regression Model (MRM), integrating out $\gamma_i$ does not yield a closed form for the marginal likelihood $p(y_i|\beta,\Sigma)$ — this forces numerical or sampling-based estimation.
Conditional independence: Given $\gamma_i$ , all observations within subject $i$ are independent; the correlation structure arises entirely through the shared random effect.
Overdispersion: Random effects absorb extra-binomial or extra-Poisson variation not explained by covariates; without them, standard GLMs underestimate standard errors for clustered data.
Nested random effects: Multiple levels of grouping (subjects within clusters, visits within subjects) add independent random effect terms at each level, each with its own variance component.

How It Works

Iterative Weighted Least Squares (IWLS) Linearization (Gamerman 1997)

At iteration $t$ , linearize around current linear predictor $\eta^{(t)} = X\beta^{(t)} + Z\gamma^{(t)}$ :

Pseudo-response: $\tilde{y}_{ij} = \eta_{ij}^{(t)} + (y_{ij} - \mu_{ij}^{(t)}) g'(\mu_{ij}^{(t)})$
Weights: $W_{ij}^{-1} = b''(\theta_{ij}^{(t)})\,[g'(\mu_{ij}^{(t)})]^2$

This converts the non-Gaussian likelihood into a locally equivalent weighted Normal model.

Weighted Least Squares (WLS) Proposal for β

With prior $\beta \sim N(a, R)$ , the working-model posterior (treating $\gamma$ fixed) is:

$m^{(t)} = (R^{-1} + X'WX)^{-1}(R^{-1}a + X'W\tilde{y})$ $C^{(t)} = (R^{-1} + X'WX)^{-1}$

Draw candidate $\beta^* \sim N(m^{(t)}, C^{(t)})$ and accept via the Metropolis ratio. The acceptance rate exceeds 95% in well-specified models because the proposal closely tracks the true posterior.

WLS Proposal for γᵢ

Treat $X\beta$ as an offset and apply the same linearization to obtain a subject-level working model. The prior $\gamma_i \sim N(0,\Sigma)$ plays the role of a Gaussian "regression" prior in the working model.

Variance Components Σ

Given all $\gamma_i$ , the full conditional for $\Sigma$ under conjugate Inverse-Wishart (IW) prior $\Sigma \sim \mathrm{IW}(\nu_0, S_0)$ is:

$\Sigma|\{\gamma_i\} \sim IW\!\left(\nu_0 + n,\ S_0 + \sum_{i=1}^n \gamma_i\gamma_i'\right)$

A direct Gibbs draw — no Metropolis-Hastings (MH) step needed.

Nested Random Effects

For two-level nesting (e.g., visit-level effects $\delta_{ij} \sim N(0,\Sigma_2)$ within subjects), add a third Metropolis block for $\delta_{ij}$ and a second IW Gibbs block for $\Sigma_2$ . Non-normal priors on $\delta_{ij}$ (e.g., Student-t via scale mixture) add one scalar auxiliary variable block.

General Design: Splines, Kriging, and GAMs (Zhao-Staudenmayer-Coull-Wand 2006)

Decompose the linear predictor as $\eta = X_R\beta_R + Z_R u_R + X_G\beta_G + Z_G u_G + Z_C u_C$ , separating (R) grouped/longitudinal, (G) general-design (splines, Generalized Additive Models (GAMs)), and (C) spatial (conditional autoregressive, CAR) components. The key unification: penalized spline coefficients $u_k \sim N(0, \sigma^2)$ are random effects with diagonal covariance — the smoothing parameter $\lambda$ equals $\sigma^2$ and is estimated from the data like any variance component. This brings GAMs, generalized additive mixed models, varying-coefficient models, and low-rank kriging all inside the GLMM framework. Full conditionals $[v_k|v_{-k}, G, y]$ are log-concave for canonical exponential families → adaptive rejection sampling (Gilks-Wild 1992) or slice sampling handles the non-conjugate blocks. WinBUGS implements this without user-written MCMC code. Hierarchical centering for grouped random effects: reparameterise $\gamma_i = \beta_R + u_i^R$ to eliminate the $\beta_R$ – $u_R$ posterior correlation that slows mixing. Variance component priors (Inverse-Gamma [IG], folded-Cauchy, Uniform) produce posteriors differing by $\le 5\%$ in this setting.

Quasilikelihood Approximation and the MQL/PQL Bias (Goldstein-Rasbash 1996)

The multilevel-modeling (IGLS/MLwiN) tradition estimates GLMMs by linearizing the link and iterating a generalized-least-squares fit of the resulting approximate Normal multilevel model (Goldstein 1986, 1991). How the random part is linearized defines two approximations: marginal quasi-likelihood (MQL) expands it about zero, while predictive/penalized quasi-likelihood (PQL) expands about the current estimated residuals $\hat u_j$ and retains second-order terms. Rodriguez-Goldman (1995) showed by simulation that first-order MQL is severely downward-biased — in both fixed effects and variance components — exactly when random-effect variances are large, clusters are small, and the response is binary. Goldstein-Rasbash (1996) show that moving to (second-order) PQL largely removes this bias while keeping the fast IGLS machinery. PQL is thus the frequentist-quasilikelihood counterpart of the Bayesian routes below; both are approximations, and the residual bias for very sparse binary data is what ultimately motivates adaptive quadrature, Laplace/INLA, and full MCMC.

Bayesian Inference via INLA (Fong-Rue-Wakefield 2010)

Pseudo-likelihood (the SAS/GLIMMIX lineage) (Wolfinger-O'Connell 1993)

The same linearize-and-iterate idea appears in the mixed-model-software tradition as pseudo-likelihood (PL). Expand the link about the current fit to form a modified dependent variable $y^*=\hat\eta+(y-\hat\mu)g'(\hat\mu)$ , which (approximately) obeys a Gaussian linear mixed model; then alternate solving the mixed-model equations for the fixed and random effects with re-estimating the variance components by REML on the current pseudo-data, rebuilding $y^*$ each pass. Because the inner engine is a general LMM, arbitrary structured covariances go on both the random effects (G-side) and the correlated errors (R-side), an overdispersion parameter is optional, and both subject-specific and population-averaged effects are recovered. This is essentially PQL by another name (and the method behind SAS PROC GLIMMIX); it shares PQL's weakness — bias for sparse binary/clustered data, dramatized by the salamander-mating benchmark — which is what motivates the quadrature/Laplace/MCMC routes.

Because a GLMM is a latent Gaussian model (Gaussian random effects on the linear predictor, non-Gaussian likelihood), its posterior can be computed by integrated nested Laplace approximations (INLA) rather than MCMC. INLA deterministically approximates the marginal posteriors of the hyperparameters and of each latent-field component by nested Laplace approximations, returning accurate posteriors in seconds-to-minutes with no chains and no convergence diagnostics. Fong-Rue-Wakefield (2010) rework the full set of Breslow-Clayton (1993) penalized-quasi-likelihood examples in this framework and stress prior specification on interpretable scales: put the prior on a random effect's marginal standard deviation, or, for a penalized-spline term, choose the smoothing-variance prior to match a target implied degrees of freedom — the reflexive $\mathrm{IG}(\epsilon,\epsilon)$ prior on a variance is not non-informative and can distort the fit or even give an improper posterior. The one regime where these deterministic approximations (like PQL itself) are least reliable is clustered binary data, where each cluster carries little information and checking against MCMC is advisable.

Why It Matters

Clustered non-Gaussian data: Binary outcomes from repeated measurements (clinical trials, surveys) and longitudinal count data (epilepsy seizures, health-service utilization) violate the normality assumption of the MRM but are naturally handled by GLMMs.
Correct inference for clustered designs: Ignoring within-subject correlation in a GLM deflates standard errors; the GLMM restores proper uncertainty quantification.
Subject-specific vs. population-average interpretation: GLMM coefficients describe changes for a specific individual (conditional on $\gamma_i$ ), which differs from the marginal (population-average) interpretation of Generalized Estimating Equations (GEE) when $g$ is nonlinear.
Full Bayesian uncertainty: The Gamerman (1997) MCMC approach provides full posteriors over $\beta$ , $\Sigma$ , and all $\gamma_i$ , including uncertainty about variance components that frequentist Restricted Maximum Likelihood (REML) treats as fixed once estimated.

Open Questions

Proposal quality at extreme responses: When binary observations are very sparse (near 0 or 1), the IWLS linearization poorly approximates the likelihood and acceptance rates can drop sharply.
High-dimensional random effects: For large $n$ (many subjects), storing and updating $\{\gamma_i\}$ individually may be slow; blocked updates or auxiliary-variable augmentation (Albert-Chib 1993b style) can accelerate sampling.
Marginal vs. conditional interpretation: Applied researchers must choose between GLMM (conditional) and GEE (marginal) based on the scientific question; the two are not interchangeable when $g$ is nonlinear (logit, log).
MCMC vs. INLA vs. HMC: Gamerman's MCMC approach has been largely superseded in practice by R-INLA (Integrated Nested Laplace Approximation (INLA), very fast) and Stan (Hamiltonian Monte Carlo (HMC)/No-U-Turn Sampler (NUTS)); the WLS-proposal idea remains instructive for hand-crafted samplers.
CAR spatial random effects: When the random effect covariance is specified via a conditional autoregressive (CAR) prior (see Conditional Autoregressive Model), the GLMM posterior may be improper; Sun-Tsutakawa-Speckman (1999) give necessary and sufficient rank conditions for propriety even when the CAR precision matrix is singular.