Cluster-Aggregated Residual Checks

Why Aggregate Residuals by Cluster?

A GLMM can fit the overall outcome distribution reasonably well while systematically underpredicting some clusters and overpredicting others. A cluster-aggregated residual check asks whether the observed outcomes in a particular cluster, taken together, are unusually high or low relative to what the fitted model predicts. This shifts attention from the individual residual to the cluster as a unit, while retaining a clear comparison with the model’s predictive variability.

A raw sum alone is not enough. A cluster with many observations will generally have more opportunity for a large sum than a small cluster, and binary or count outcomes have discrete, often asymmetric distributions. The key step is therefore to compare each observed cluster discrepancy with discrepancies generated under the fitted model, using simulations that reproduce its assumptions and the intended prediction target.

Definition and Model Setup

Index clusters by $j=1,\ldots,J$, and observations within cluster $j$ by $i=1,\ldots,n_j$. Let $Y_{ij}$ be the outcome, $b_j$ the cluster random effect, and $x_{ij}$ and $z_{ij}$ the fixed- and random-effect design vectors. A common GLMM specification is

g(\mu_{ij}(b_j))=x_{ij}^{\mathsf T}\beta+z_{ij}^{\mathsf T}b_j, \qquad Y_{ij}\mid b_j \sim \mathcal{F}(\mu_{ij}(b_j),\phi).

Here $\mu_{ij}(b_j)=E(Y_{ij}\mid b_j)$, $g$ is the link function, and $\mathcal{F}$ denotes the chosen conditional outcome distribution. The random effects have a specified distribution, often $b_j\sim N(0,D)$. Conditional on the random effects and model parameters, outcomes are usually assumed independent within a cluster. If the model specifies dependence across clusters, that joint dependence must be retained in the simulations.

Definition: Given fitted means $\hat\mu_{ij}$, the cluster-aggregated observed-minus-fitted discrepancy for cluster $j$ is $$ D_j=\sum_{i=1}^{n_j}(Y_{ij}-\hat\mu_{ij}). $$ A positive $D_j$ means the cluster’s observed total exceeds its fitted total; a negative $D_j$ means it falls below. The units are outcome units, such as events or successes.

The fitted mean must match the diagnostic question. A conditional fitted mean uses an estimated cluster effect, such as $\hat b_j$, and describes the outcome conditional on that cluster’s fitted effect. A marginal fitted mean averages over the random-effect distribution and describes the population-average prediction. These are not generally equal in nonlinear GLMMs. The earlier material on conditional and marginal fitted values distinguishes the targets; whichever is chosen, the simulation must use the same target and reproduce the corresponding fitting or prediction procedure.

For a conditional check, one may compare observed outcomes with means conditional on estimated effects, but simply simulating new outcomes at fixed $\hat b_j$ treats those estimates as known. For a new-cluster or population-level check, simulate random effects from the fitted distribution and then simulate outcomes. For an existing-cluster check, a principled predictive simulation can draw effects from their fitted conditional distribution given the cluster data. A parametric bootstrap may also refit the model in each replicate. The method of simulation changes the reference distribution and must be reported.

Constructing the Reference Distribution

For each simulated replicate $r=1,\ldots,B$, generate data under the fitted model, calculate the fitted means using the selected procedure, and compute the cluster discrepancy $D_j^{(r)}$ by the same formula as for the observed data. Comparing $D_j$ with the simulated values accounts for cluster size, discreteness, and model-implied variability more directly than comparing all raw sums to a common cutoff.

If the fitted means and parameters were known, and observations were conditionally independent given $b_j$, then the conditional expectation of the discrepancy would be zero and its conditional variance would be the sum of the outcome variances. For example, for conditionally Bernoulli outcomes that variance is $\sum_i \mu_{ij}(b_j)\{1-\mu_{ij}(b_j)\}$. In practice, fitted parameters and random effects are estimated, and the marginal outcomes may be dependent. Simulation is useful precisely because a simple variance formula may not capture those features.

Key principle: The observed discrepancy and every simulated discrepancy must be computed in the same way, under a simulation scheme that represents the model and prediction target being checked. Compare $D_j$ with its own simulated reference distribution, not with an arbitrary universal threshold.

Worked Example: Aggregating Poisson Count Residuals

Suppose a hypothetical fitted Poisson GLMM predicts event counts for one clinic over four periods as $1.2,\ 0.8,\ 1.5,\ 0.5$. The observed counts are $3,\ 0,\ 2,\ 1$. Calculate the cluster discrepancy.

The observed total is $3+0+2+1=6$, and the fitted total is $1.2+0.8+1.5+0.5=4.0$. Therefore

D_j=\sum_i(Y_{ij}-\hat\mu_{ij}) =(3-1.2)+(0-0.8)+(2-1.5)+(1-0.5)=2.0.

The clinic has two more events than its fitted total. This is not yet evidence of lack of fit: the total is random, and its variability depends on the fitted means and the random-effect treatment. Suppose a simulation of 2,000 replicates under the stated fitted-model procedure gives a central 95% interval from $-4$ to $5$ events for this clinic’s discrepancy. The observed value $2.0$ lies inside that interval, so this check does not flag the clinic as unusually high or low at that reference level.

The comparison is specific to this clinic’s covariates, exposure, cluster size, and simulation scheme. It does not establish that the model is correct; it says only that this particular aggregate discrepancy is not unusual under the chosen simulations.

Scaling and Simulation-Based Interpretation

A raw discrepancy is easy to interpret in outcome units but may be hard to compare across clusters. One option is a simulation-standardized discrepancy:

Z_j^{\mathrm{sim}}= \frac{D_j-\bar D_j^{\mathrm{sim}}}{s_j^{\mathrm{sim}}}, \qquad \bar D_j^{\mathrm{sim}}=\frac{1}{B}\sum_{r=1}^{B}D_j^{(r)}, \qquad s_j^{\mathrm{sim}}= \sqrt{\frac{1}{B-1}\sum_{r=1}^{B}(D_j^{(r)}-\bar D_j^{\mathrm{sim}})^2}.

This scale measures how far the observed discrepancy is from the simulated mean in simulated standard-deviation units. It is a descriptive standardization, not necessarily a normally distributed statistic. When the simulated distribution is skewed or discrete, its quantiles are usually more informative than a rule such as “flag values beyond 2.”

Worked Example: Unequal Cluster Sizes in a Binary GLMM

Consider two hypothetical clinics with fitted success probabilities of $0.20$ for each person. Clinic A has $n_A=5$ people and observes 3 successes; clinic B has $n_B=20$ and observes 7 successes. Compute their observed-minus-fitted discrepancies.

For A, the fitted total is $5(0.20)=1.0$, so $D_A=3-1.0=2.0$ successes. For B, the fitted total is $20(0.20)=4.0$, so $D_B=7-4.0=3.0$ successes. The raw discrepancies show that B has the larger excess in success counts, but they do not show which is more unusual relative to its own expected variation.

Under an illustrative independent Bernoulli reference with probability $0.20$, the variance of a sum of $n$ outcomes is $n(0.20)(0.80)$. Thus the standard deviations are $\sqrt{5(0.20)(0.80)}=\sqrt{0.8}=0.894$ for A and $\sqrt{20(0.20)(0.80)}=\sqrt{3.2}=1.789$ for B. The variance-standardized discrepancies are approximately $2.0/0.894=2.24$ for A and $3.0/1.789=1.68$ for B.

This calculation is only an illustration under independent, known-probability assumptions. A GLMM simulation should instead reproduce the fitted random-effect structure and the actual fitting procedure. The example shows why a larger raw sum need not be the more extreme cluster check.

Reading Simulated Comparisons

For a prespecified upper-tail concern, an empirical simulation tail proportion can be calculated as

\hat q_j= \frac{1+\sum_{r=1}^{B} I(D_j^{(r)}\ge D_j)} {B+1}.

The added one prevents a reported value of zero when no simulated discrepancy exceeds the observed one. This proportion estimates how often the simulated discrepancy is at least as large as the observed discrepancy under the simulation scheme. It is not the probability that the model is true, and it is not automatically a calibrated frequentist p-value, especially when parameters were estimated from the same observations.

Worked Example: Interpreting an Upper-Tail Simulation Check

For a hypothetical clinic, the observed aggregate discrepancy is $D_j=8.6$ excess events. In $B=1{,}000$ simulations, only 3 replicated discrepancies are at least $8.6$. If an unusually high total was the prespecified concern, calculate the empirical upper-tail proportion.

\hat q_j=\frac{1+3}{1{,}000+1} =\frac{4}{1{,}001}\approx 0.00400.

Only about 0.4% of the simulated discrepancies are at least as large as the observed one, using the plus-one estimate. This is a reason to investigate possible high-side lack of fit for this clinic under the specified model and simulation procedure. It does not prove that the model is false or identify the cause. If the concern was two-sided, the assessment should instead compare the observation with both tails or with a suitable central simulation interval; a one-sided calculation cannot answer a two-sided question.

Workflow, Diagnostics, and Limitations

State the target.
Decide whether the check is conditional on existing cluster effects, marginal over random effects, or predictive for new clusters. Identify the outcome scale and the meaning of a high or low cluster total.

Calculate observed discrepancies.
For each cluster, sum outcomes and subtract the corresponding sum of fitted means. Use consistent observations, offsets, and exposure definitions.

Generate model-based replicates.
Simulate random effects and outcomes as required by the target. If the intended procedure re-estimates parameters or effects, repeat that operation for each replicate.

Compare and investigate.
Use cluster-specific simulated intervals, quantiles, or tail proportions. Examine flagged clusters for data problems or systematic features, then use other diagnostics to investigate possible explanations.

A plot of observed $D_j$ against cluster-specific simulated intervals can reveal clusters whose aggregate outcome is unusual under the model. A directional pattern across many clusters may be more informative than a single flag. However, examining many clusters increases the chance of seeing at least one extreme value even when the model is adequate. The individual cluster comparisons are generally dependent when clusters share fitted parameters or correlated random effects. Do not interpret a collection of nominal tail proportions as independent tests or claim a family-wide error rate without an appropriate multiplicity procedure.

Several failures can make the comparison misleading. Simulating outcomes at estimated random effects without reflecting their estimation uncertainty can give an overly narrow reference distribution. Simulating new independent random effects when the model specifies spatial or temporal correlation fails to preserve its dependence structure. Comparing raw sums across clusters without cluster-specific references ignores differences in size and fitted risk. Finally, a discrepancy can be driven by a data error, an omitted predictor, a misspecified random-effect distribution, or another feature; the aggregate check alone does not distinguish among them.

Cluster aggregation complements rather than replaces individual residual displays. Individual residuals can show opposing within-cluster departures that cancel in a sum; an aggregate check can expose a coherent cluster-level shift that is less apparent when inspecting observations separately. Conversely, an unusual sum does not reveal which observation or pattern produced it. The main alternative is therefore not a single competing statistic but a set of diagnostics aimed at different scales. Cluster-aggregated checks are descriptive model checks, not a substitute for assessing the model’s assumptions or validating its predictions.

Common Mistakes

Calling a large sum automatically unusual. A discrepancy of 5 may be ordinary in a large cluster and surprising in a small one. Compare it with the cluster’s own simulated distribution.
Mixing conditional and marginal predictions. A sum based on conditional fitted means should not be judged against simulations for a different prediction target.
Treating a tail proportion as the probability the model is true. It is a model-based simulation comparison under specified assumptions, not a posterior probability that the null model is correct.
Using only a one-sided check after seeing the data. Decide whether high, low, or either direction matters; if either direction is relevant, use a two-sided comparison.
Interpreting one flagged cluster as a model-wide diagnosis. A flag identifies a discrepancy worth investigating, not its cause, and multiple cluster checks create multiplicity.

Key takeaway: Sum observed-minus-fitted outcomes within each cluster, then judge each sum against a simulated reference distribution that matches the GLMM, the fitting procedure, and the prediction target. The result identifies unusual aggregate discrepancies, not definitive proof of a particular model failure.

Check Your Understanding

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

A cluster has observed total 12 and fitted total 9.5. What is its aggregate discrepancy, and what does its sign mean?
Why can raw aggregate discrepancies be difficult to compare between clusters of different sizes?
In 499 simulations, 9 replicated discrepancies are at least as large as the observed discrepancy. Calculate the plus-one upper-tail proportion.
Explain why simulating new random effects may be inappropriate for a check conditional on particular existing clusters.

Solutions. (1) $D_j=12-9.5=2.5$ outcomes. It is positive, so the observed total exceeds the fitted total by 2.5. (2) Larger clusters generally have more outcome variation and more opportunity for a large sum; the simulated reference distribution accounts for cluster-specific size and fitted risk. (3) $\hat q_j=(1+9)/(499+1)=10/500=0.0200$, or 2.00%. This is a simulation tail proportion under the specified procedure, not automatically a calibrated frequentist p-value. (4) Drawing new effects targets variation for new clusters or a population-level prediction. For an existing-cluster conditional check, the simulation should reflect the information about that cluster’s effect, for example through an appropriate conditional predictive distribution or a procedure that repeats the relevant estimation steps.