Common Misinterpretations of GLMM Residuals

What a GLMM residual can tell you

A residual is a discrepancy between an observed outcome and a model-based expectation. In a generalized linear mixed model (GLMM), that discrepancy depends on which expectation is used. A fitted value conditional on estimated random effects answers a different question from a prediction averaged over the random-effects distribution. Neither kind of residual is automatically a normally distributed error.

This matters especially for binary and sparse-count outcomes. Such observations have a restricted, discrete set of possible values, while a normal distribution is continuous. A residual plot can still reveal patterns that deserve investigation, but “the residuals are not normal” is not, by itself, evidence that a GLMM is wrong. The relevant question is whether the observed discrepancies are compatible with the model and with the diagnostic target.

For a Bernoulli observation with fitted probability \(\hat p_{ij}\), the conditional variance is \(\hat p_{ij}(1-\hat p_{ij})\). For a Poisson observation with fitted mean \(\hat\mu_{ij}\), it is \(\hat\mu_{ij}\). These variances define the scale of Pearson residuals under the assumed conditional response model. They do not turn a discrete outcome into a continuous normal variable.

These interpretations assume that the response family and link are appropriate, the mean structure is adequately specified, and the fitted values correspond to the conditioning level being discussed. They also depend on how uncertainty from estimated parameters and random effects is handled. Ordinary residual plots typically treat fitted quantities as if known, so they are descriptive diagnostics rather than exact tests with universally calibrated reference distributions.

Why normal-looking residuals are not the target

The distribution of a residual follows from the distribution of the outcome and the fitted mean. For a Bernoulli outcome, \(Y_{ij}\) can only be 0 or 1. With a fixed fitted probability \(p\), its standardized Pearson residual can therefore take only two possible values: $$ r^P= \begin{cases} -\sqrt{\dfrac{p}{1-p}}, & Y=0,\\[6pt] \sqrt{\dfrac{1-p}{p}}, & Y=1. \end{cases} $$ There is no reason for these two values to form a bell-shaped distribution. For a sparse Poisson outcome, residuals are also tied to integer counts and may show a concentration near zero with a long positive tail.

The standardized Bernoulli residual has mean zero and variance one when \(p\) is the true probability and is treated as known. To see this, the residual is \((Y-p)/\sqrt{p(1-p)}\), so its expectation is zero because \(E(Y)=p\). Its variance is \(\operatorname{Var}(Y)/[p(1-p)]=1\). But mean zero and variance one do not imply normality: the residual still has just two possible values. This distinction is one reason that a normal quantile plot can look poor even when the response distribution is correctly specified.

Conditional residuals are not population-level evidence

In a GLMM, a conditional fitted value uses the estimated random effect for the observation’s cluster. A population-average fitted mean instead averages the response mean over the random-effects distribution. For a nonlinear link, these values generally differ: averaging \(g^{-1}(x^\mathsf{T}\beta+z^\mathsf{T}b)\) over \(b\) is generally not the same as evaluating \(g^{-1}(x^\mathsf{T}\beta)\). Consequently, a residual based on conditional fitted values assesses the model’s fit to observations given estimated cluster effects; it does not directly show how well the model predicts an average or a new cluster in the population.

Estimated cluster effects are informed by the observations in those clusters and are typically shrunk toward the overall distribution. That data use makes conditional residuals useful for examining within-cluster fit, but it also means they are not independent, untouched evidence about the population distribution of cluster effects. A plot of conditional residuals should not be used as a substitute for a population-level predictive check or a direct assessment of the random-effects distribution.

Dispersion ratios are summaries, not automatic verdicts

For count or grouped-binomial GLMMs, a common descriptive summary is the Pearson statistic divided by a residual degrees-of-freedom value: $$ \hat\phi=\frac{X^2_P}{df}, \qquad X^2_P=\sum_{ij}(r^P_{ij})^2. $$ A ratio near one is compatible with the assumed variance scale under some settings. A value appreciably above one can flag extra variation relative to that scale; a value below one can flag less variation. But GLMM residual degrees of freedom are not always uniquely defined by simply subtracting the number of fixed-effect coefficients from the number of observations. Random effects, estimation, grouping, and model structure complicate that calculation. The ratio is often best treated as a rough diagnostic summary, not a universal test statistic with a guaranteed reference distribution.

Even a clearly elevated ratio does not identify its cause. Possible explanations include a misspecified mean structure, unmodeled dependence, omitted heterogeneity, or an unsuitable response distribution. It should be considered alongside residual patterns, the study design, and other relevant checks. Conversely, a ratio near one does not demonstrate that every part of the model is correct: positive and negative discrepancies can be summarized into a single number, concealing structure.

Using residuals without overreading them

A useful residual review begins by stating the target: conditional fit for observed clusters, or prediction averaged over random effects. Then use summaries suited to the response scale and inspect whether discrepancies show structure rather than demanding normality. For example, residuals plotted against fitted values or predictors can reveal systematic trends; grouping them by cluster can reveal persistent within-cluster patterns. Such plots are clues, not proofs, and sparse data can make apparent patterns unstable.

For a Poisson observation with conditional fitted mean \(\hat\mu=0.40\), the Pearson residual for a zero count is $$ r^P=\frac{0-0.40}{\sqrt{0.40}}=-0.6325. $$ The deviance residual for a Poisson zero is \(-\sqrt{2\hat\mu}\), giving \(-\sqrt{0.80}=-0.8944\). The different values reflect different residual definitions, not conflicting calculations. A residual plot’s appearance can depend on the chosen definition, so conclusions should not hinge on a single plot or threshold.

The main alternatives answer related but distinct questions. Response residuals retain the outcome’s original units. Pearson residuals scale discrepancies by the model-based conditional standard deviation. Deviance residuals use likelihood contributions and can highlight observations poorly represented by the fitted response distribution. Predictive checks compare replicated or predicted outcomes with observed features and can target a broader distributional or population question. None of these tools makes residual normality a general GLMM requirement; selection depends on the diagnostic question.

A practical failure mode is to treat a single outlying residual as a diagnosis. A large residual may reflect a rare but model-compatible outcome, parameter estimation, or an influential observation; it may also point to a genuine mismatch. Investigate the observation and its context, but avoid removing it solely because it is inconvenient. Another failure mode is to use conditional residuals to claim that new clusters will be predicted well. That claim needs a diagnostic target that represents new-cluster prediction, not merely fit to clusters whose outcomes contributed to their estimated random effects.