Tutorials › Generalized Linear Mixed Models (GLMM) › Pearson Dispersion Statistics: Interpretation and Limits

Diagnostics and Predictive Checks for Sparse and Overdispersed Outcomes · Tutorial 624 of 698

Pearson Dispersion Statistics: Interpretation and Limits

Learn how to calculate Pearson-based dispersion summaries in GLMMs and interpret them cautiously when random effects are estimated or outcomes are sparse.

Intermediate 10 min read

What You'll Learn

  • Calculate Pearson residuals and the raw Pearson statistic for conditional GLMM fits.
  • Interpret a Pearson statistic divided by a nominal residual degrees of freedom as a rough summary, not a universal test.
  • Distinguish conditional variance from marginal variance when random effects are present.
  • Explain how estimating random effects can make conditional residuals appear smaller.
  • Recognize why sparse outcomes and finite samples weaken reference-value interpretations.
  • Compare Pearson summaries with residual plots, deviance summaries, and simulation-based checks.

Why a Pearson Dispersion Summary Needs Context

A Pearson dispersion summary compresses many residuals into one number. It is useful for noticing whether residual variation looks large or small relative to a fitted model, but in a generalized linear mixed model (GLMM) the number does not have a universal reference distribution. The model may include estimated random effects, the outcomes may be sparse, and the observations may be independent only after conditioning on those random effects. Each feature affects what the summary says.

The key distinction is between calculating a Pearson statistic and calibrating it. The calculation is usually direct once the variance used in the residual is specified. Calibration—deciding whether a value is unusually large or small—requires care about the model, estimation, sample size, and the chosen degrees of freedom.

Definition and Model Conditions

Let \(Y_{ij}\) denote an outcome for observation \(i\) in cluster \(j\), and let \(b_j\) denote that cluster's random effect or random-effects vector. A common GLMM specification assumes that, conditional on the random effects and covariates, the outcomes are independent and follow a specified response distribution. Write the conditional mean as \(\mu_{ij}(b_j)=E(Y_{ij}\mid b_j)\), and the conditional variance as \(\operatorname{Var}(Y_{ij}\mid b_j)=\phi V(\mu_{ij}(b_j))\), where \(V\) is the variance function and \(\phi\) is a dispersion parameter when the family allows one. For standard conditional Poisson and binomial models, \(\phi\) is conventionally fixed at 1.

Definition: Given fitted conditional means \(\hat{\mu}_{ij}\), a Pearson residual is \(r_{ij}^{P}=(y_{ij}-\hat{\mu}_{ij})/\sqrt{V(\hat{\mu}_{ij})}\) when the fitted conditional dispersion is 1. If a specified dispersion \(\phi\) is included in the residual scale, the denominator is \(\sqrt{\phi V(\hat{\mu}_{ij})}\). The raw Pearson statistic is \(Q_P=\sum_{j}\sum_i (r_{ij}^{P})^2\). A commonly reported summary is \(Q_P/df_{\mathrm{resid}}\), with the reported residual degrees of freedom stated explicitly.

For a grouped binomial observation with \(m_{ij}\) trials and fitted success probability \(\hat p_{ij}\), the conditional mean is \(m_{ij}\hat p_{ij}\), and the conditional variance is \(m_{ij}\hat p_{ij}(1-\hat p_{ij})\). Its Pearson residual is therefore \((y_{ij}-m_{ij}\hat p_{ij})/\sqrt{m_{ij}\hat p_{ij}(1-\hat p_{ij})}\). For an individual Bernoulli outcome, set \(m_{ij}=1\).

These definitions require fitted means and variances that correspond to the same target. In particular, a conditional Pearson residual uses a conditional mean and a conditional variance; a marginal residual uses a marginal mean and marginal variance. Substituting a marginal mean into a conditional variance formula, or the reverse, does not generally produce a correctly standardized residual.

Why Unit Dispersion Is Only a Starting Point

If the true conditional mean and variance are known and the conditional distribution is correctly specified, then each squared Pearson residual has conditional expectation 1. For one observation,

$$ E\left[\left.\frac{(Y_{ij}-\mu_{ij}(b_j))^2}{V(\mu_{ij}(b_j))}\right|b_j\right] = \frac{\operatorname{Var}(Y_{ij}\mid b_j)}{V(\mu_{ij}(b_j))} =1 $$

when \(\phi=1\). Thus the expected sum is the number of observations if the parameters and random effects are known. In an ordinary GLM, parameter estimation uses information from the same observations and typically reduces the residual variation available for the sum. Under suitable large-sample conditions, the Pearson statistic is often compared approximately with a chi-square distribution having residual degrees of freedom, commonly \(n-p\) for \(n\) independent observations and \(p\) fitted fixed-effect parameters.

That familiar approximation is not an automatic GLMM result. Conditional residuals use estimated random effects, which are typically shrunk toward the population mean, and the fixed effects and random-effects covariance parameters are estimated too. The resulting residuals are not independent, known-parameter residuals. There is no single universally correct rule for subtracting degrees of freedom from \(Q_P\) across all GLMMs. Software may report a nominal value such as the number of observations minus the number of fixed-effect coefficients; treat that convention as a reference summary, not as an exact finite-sample calibration.

Worked Calculations

Worked Example: A Pearson Statistic for Grouped Binomial Data

Suppose a hypothetical conditional binomial GLMM gives fitted probabilities \(0.20, 0.30, 0.50,\) and \(0.70\) for four clusters, each with \(m=10\) trials. The observed success counts are \(2,4,6,\) and \(8\). Treat these as fitted conditional probabilities for the calculation.

For the first cluster, the fitted count is \(10(0.20)=2\), the variance is \(10(0.20)(0.80)=1.6\), and the residual is \((2-2)/\sqrt{1.6}=0\). Its contribution to \(Q_P\) is 0.

For the second cluster, the fitted count is 3, the variance is \(10(0.30)(0.70)=2.1\), and the squared residual is \((4-3)^2/2.1=0.4762\). For the third, the fitted count is 5, the variance is \(10(0.50)(0.50)=2.5\), and the contribution is \((6-5)^2/2.5=0.4000\). For the fourth, the fitted count is 7, the variance is \(10(0.70)(0.30)=2.1\), and the contribution is \((8-7)^2/2.1=0.4762\).

Therefore \(Q_P=0+0.4762+0.4000+0.4762=1.3524\). This is the raw conditional Pearson statistic for these four grouped observations. A ratio cannot be interpreted without specifying how residual degrees of freedom were defined, and even a reported ratio near 1 would not by itself establish that the GLMM is adequate.

Worked Example: Sparse Bernoulli Outcomes Can Dominate the Sum

Consider a hypothetical set of four binary outcomes with fitted conditional probabilities all equal to \(0.10\), and observed outcomes \(0,0,1,0\). For a Bernoulli observation, the squared Pearson residual is \((y-\hat p)^2/\{\hat p(1-\hat p)\}\).

For each of the three zeros, the contribution is \((0-0.10)^2/\{0.10(0.90)\}=0.01/0.09=0.1111\). For the one, the contribution is \((1-0.10)^2/\{0.10(0.90)\}=0.81/0.09=9\). Thus \(Q_P=3(0.1111)+9=9.333\), allowing for rounding.

One event accounts for most of the total. The calculation is correct, but a large total in such a small sparse sample is not enough to diagnose the cause: it may reflect an unusual event, a poorly calibrated fitted probability, or a model feature not captured by this summary. A simple large-sample chi-square comparison is particularly fragile here because the data provide few events with which to assess the outcome distribution.

Worked Example: Conditional and Marginal Variances Differ

Suppose a hypothetical Poisson random-intercept model has \(Y\mid b\sim\operatorname{Poisson}(2e^b)\), where \(b\sim N(0,0.25)\). Conditional on \(b\), the variance equals the mean. Marginally, however, the outcome variance includes variation in the conditional mean across clusters.

For a lognormal random effect, \(E(e^b)=e^{0.25/2}=e^{0.125}\). The marginal mean is therefore \(\mu_M=2e^{0.125}\approx2.2663\). By the law of total variance, \(\operatorname{Var}(Y)=E\{\operatorname{Var}(Y\mid b)\}+\operatorname{Var}\{E(Y\mid b)\}\). Here this becomes \(\mu_M+\mu_M^2(e^{0.25}-1)\).

$$ \operatorname{Var}(Y) \approx 2.266297+(2.266297)^2(e^{0.25}-1) \approx 2.266297+5.1361(0.284025) \approx 3.7251 $$

The conditional variance for a cluster with \(b=0\) is 2, while the marginal variance is about 3.725. For an observed count \(y=4\), the conditional Pearson residual at \(b=0\) is \((4-2)/\sqrt{2}=1.414\), whereas standardizing the marginal difference from the marginal mean gives \((4-2.2663)/\sqrt{3.7251}\approx0.898\). These calculations answer different questions; one variance cannot be substituted for the other without changing the diagnostic target.

Estimated Random Effects, Degrees of Freedom, and Diagnostics

A fitted conditional mean often plugs an estimated random effect into the linear predictor. Because the estimate is informed by the cluster's own outcomes and is usually shrunk toward zero, the fitted mean can move toward those outcomes. This can make conditional residuals smaller than residuals evaluated at the unknown true random effects. The resulting \(Q_P\) may therefore look reassuring even when the random-effects distribution or the overall model structure is inadequate.

A useful derivation of the limitation is to compare the known-parameter identity with the fitted calculation. The identity uses \(Y_{ij}-\mu_{ij}(b_j)\), where \(b_j\) and the model parameters are fixed at their true values. A fitted residual instead uses \(Y_{ij}-\hat\mu_{ij}(\hat b_j)\). Since \(\hat b_j\) depends on the data, the numerator is no longer a centered outcome conditional on a fixed mean. Consequently, the expectation-one calculation above does not establish that fitted squared residuals sum to a chi-square variable with a particular degrees of freedom.

When reviewing a Pearson summary, inspect the calculation and its surrounding evidence rather than applying a mechanical cutoff. Check that the model family and variance function are plausible, that fitted means and denominators are on the same conditional or marginal scale, and that convergence and boundary estimates do not undermine the fit. Then examine residual patterns by fitted value, covariate, and cluster, and check whether a small number of sparse outcomes dominates the sum. A global statistic cannot reveal where a mismatch occurs.

The Pearson statistic is one diagnostic among several. Pearson residual plots retain information about direction and patterns that a sum of squares discards. Deviance summaries use a different residual construction, while simulation-based residual checks can assess a statistic against replicated data generated under a specified fitted model. None removes the need to define the prediction target and account for estimation; the methods answer related but nonidentical questions.

Common Mistakes and Limits

  • Calling \(Q_P/df_{\mathrm{resid}}\) a definitive overdispersion test. It is a rough scale summary unless a valid reference distribution and degrees of freedom are justified for the fitted model.
  • Assuming the GLM rule \(n-p\) automatically applies. A GLMM estimates random-effects parameters and uses estimated cluster effects; nominal software degrees of freedom do not generally yield an exact chi-square reference.
  • Mixing conditional and marginal quantities. A conditional mean must be paired with its conditional variance; marginalization over random effects changes both the mean and variance.
  • Reading a large value as proof of one particular cause. The statistic flags a discrepancy in aggregate but cannot distinguish an incorrect mean structure, variance structure, random-effects specification, or a few influential observations.
  • Reading a value near 1 as proof of fit. Sparse data, few clusters, estimated random effects, or offsetting residual contributions can conceal important problems.

The principal limitation is therefore interpretive, not arithmetic. Pearson-based summaries are easy to calculate and useful for organizing a diagnostic review, but their reference values are not universal for fitted GLMMs. In sparse or finite samples, describe the statistic, the fitted scale, and the degrees-of-freedom convention, then support any conclusion with more detailed diagnostics.

Check Your Understanding

Answer each question, then compare your work with the solutions.

  1. A Bernoulli observation has \(y=0\) and fitted probability \(\hat p=0.20\). Calculate its squared Pearson residual.
  2. Why does the identity \(E\{(Y-\mu)^2/V(\mu)\mid b\}=1\) not prove that a fitted GLMM Pearson statistic has an exact chi-square distribution?
  3. In the random-intercept example, why is the marginal variance larger than the conditional variance at \(b=0\)?
  4. A fitted model has a large Pearson statistic, most of which comes from one rare event. State one appropriate conclusion and one conclusion that is not justified.

Solutions. (1) \((0-0.20)^2/\{0.20(0.80)\}=0.04/0.16=0.25\). (2) The identity assumes fixed true means and parameters conditional on the random effects. Fitted means use estimated parameters and data-dependent estimated random effects, so the residuals need not have that reference distribution. (3) Marginal variance adds the variance of the conditional mean across random effects to the expected conditional variance; clusters differ in their conditional rates. (4) It is appropriate to report that the aggregate statistic is strongly affected by a rare event and warrants investigation. It is not justified to declare a particular model component wrong, or to claim a definitive test result, from that statistic alone.

Clinical Biostats