Pearson Residuals and Their Scaling

Why Scale a Response Residual?

A response residual is the observed outcome minus its fitted mean. But a difference of one has different significance in different parts of a GLMM: for a rare binary outcome it may be large, while for a count with a large expected value it may be small. Pearson residuals address this by expressing the difference in units of the model-implied standard deviation.

In a GLMM, that scale depends on the chosen conditioning. A conditional fitted mean uses a specified value of the random effects; a marginal fitted mean averages over their distribution. The corresponding variances are not generally the same. This tutorial focuses on how to derive and interpret Pearson residuals on either clearly stated scale, and on why discrete residuals need not resemble continuous normal data.

For a conventional GLMM, the response distribution is specified conditionally on random effects. Write \(b\) for the random effects and \(\mu_{ij}(b)=E(Y_{ij}\mid b)\) for the conditional mean. If the conditional variance is \(\operatorname{Var}(Y_{ij}\mid b)=\phi V(\mu_{ij}(b))\), the conditional Pearson residual evaluated at fitted parameters and fitted random effects is

Here \(\hat b\) denotes the fitted random effects, \(V(\cdot)\) is the family’s variance function, and \(\phi\) is a dispersion factor when the model includes one. For Bernoulli, binomial, and Poisson models, the standard variance specifications usually have \(\phi=1\). A different scale is appropriate if the model explicitly allows a different conditional variance.

Deriving the Common Discrete-Outcome Forms

The derivation is direct: specify the conditional mean and variance for the outcome, substitute them in the Pearson-residual definition, and simplify. These formulas apply when the stated response distribution and its variance are appropriate conditional on the random effects. In a standard GLMM, observations are assumed conditionally independent given those effects, unless the model specifies additional dependence.

For grouped binomial data, \(y_{ij}\) is the observed number of successes out of \(m_{ij}\) trials, not the observed proportion. If the data are instead recorded as a proportion \(y_{ij}/m_{ij}\), both the residual and its variance must be expressed on the proportion scale. Mixing a count numerator with a proportion-scale denominator gives an incorrect residual.

Raw, Variance-Scaled, and Leverage-Adjusted Residuals

The term “standardized residual” is used in more than one way. To avoid ambiguity, distinguish the following quantities.

Thus, a Pearson residual is already scaled by the fitted response variance. Some sources call it a standardized Pearson residual; others reserve “standardized” for a further adjustment. Always state which definition is being used.

In a simple GLM, a commonly used approximate leverage adjustment is \(r_{P,i}/\sqrt{1-h_{ii}}\), where \(h_{ii}\) is the \(i\)th diagonal element of an appropriate hat matrix. It attempts to account for the fact that fitted values are estimated from the same data, which can make ordinary fitted residuals less variable than residuals computed with known parameters. In a GLMM, the analogue depends on how the random effects and covariance parameters are estimated and on the fitting approximation. There is no single leverage adjustment that applies unchanged to every GLMM. Use such a correction only when its definition and approximation match the fitting method.

Why Discreteness Changes the Apparent Spread

If the parameters and random effects were known and the conditional model were correct, a Pearson residual would have conditional mean zero and variance one: $$ E\left[\frac{Y_{ij}-\mu_{ij}(b)} {\sqrt{\operatorname{Var}(Y_{ij}\mid b)}}\middle|b\right]=0, \qquad \operatorname{Var}\left[\frac{Y_{ij}-\mu_{ij}(b)} {\sqrt{\operatorname{Var}(Y_{ij}\mid b)}}\middle|b\right]=1. $$ This follows by subtracting the conditional mean and dividing by the conditional standard deviation. It does not imply that each residual is normally distributed, nor that the residuals computed after fitting are independent with exactly unit variance.

For a Bernoulli observation, only two residual values are possible at a given fitted probability: $$ r_P(1)=\sqrt{\frac{1-p}{p}}, \qquad r_P(0)=-\sqrt{\frac{p}{1-p}}. $$ When \(p\) is small, the positive value can be large, while the much more common zero outcome has a small negative residual. The residual distribution is therefore discrete and asymmetric even though its theoretical mean is zero and variance is one. For sparse Poisson counts, possible residual values are also discrete; with fitted mean \(\mu\), the smallest possible value is \(-\sqrt{\mu}\), attained at zero, while positive values increase in steps of \(1/\sqrt{\mu}\).

These features matter in practice. A histogram of binary Pearson residuals cannot become a smooth bell-shaped distribution merely because the sample is large: many observations can share one of a small number of values, especially when fitted probabilities are similar. In sparse counts, a few uncommon large positive residuals can coexist with many small negative residuals. Visual departure from normality is not, on its own, proof that the fitted mean or variance model is wrong.

Choosing the Scale and Reading Diagnostics

A conditional Pearson residual uses the conditional mean and response variance at fitted random effects. It asks whether an observation is surprising relative to the fitted response distribution at that random-effect value. A marginal residual instead uses a mean averaged over random effects and a variance that includes both response variation and variation in the conditional mean. By the law of total variance, $$ \operatorname{Var}(Y_{ij}) = E_b\{\operatorname{Var}(Y_{ij}\mid b)\} + \operatorname{Var}_b\{E(Y_{ij}\mid b)\}. $$ Consequently, a marginal Pearson residual is not obtained by simply reusing the conditional denominator with a marginal mean. The mean and variance must both refer to the same target distribution. These scales answer different questions; neither is universally preferable.

Fitted conditional residuals also use estimated fixed effects, variance parameters, and random effects. The theoretical mean-zero, variance-one result assumes the relevant parameters and conditioning values are known. Estimation, shrinkage of fitted random effects, and dependence among observations within clusters can alter the residuals’ distribution. Therefore, do not treat a collection of fitted Pearson residuals as independent standard normal observations without a specific justification.

When using Pearson residuals diagnostically, check that each residual was calculated on the intended response scale, with the correct denominator and conditioning convention. Look for repeated large values, systematic signs, or patterns against fitted means or covariates, while remembering that discreteness itself creates visible bands and asymmetry. An isolated large residual is a prompt to check the observation, model specification, and data processing; it is not a stand-alone model-selection rule.

The raw response residual remains useful because it is in the outcome’s original units: events, successes, or counts. Pearson scaling makes residual magnitudes more comparable across observations with different fitted variances. The scaling does not remove discreteness, account automatically for leverage, or repair a misspecified conditional variance. These distinctions are especially important for sparse outcomes, where a visually dramatic Pearson residual may arise from a rare but model-possible observation.