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.

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.