Simulation-Based Residuals for GLMMs

Why Compare Observed Outcomes With Simulations?

A response residual compares an observed outcome with a fitted value, but in a generalized linear mixed model (GLMM) the response can be discrete, its variability depends on its mean, and observations within a cluster are related. A residual measured in raw response units may therefore be hard to compare across observations. Simulation-based residuals address this by asking where each observed outcome falls among outcomes the fitted model could have generated.

The central object is not a single simulated data set. It is a collection of replicate outcomes generated under the fitted model, using the same design and preserving the model’s dependence structure. The observed value is then compared with that simulated reference. For discrete outcomes, many replicates can equal the observed value, so the residual construction must handle ties explicitly.

Generating Replicates From a Fitted GLMM

A GLMM specifies a distribution for the responses conditional on latent random effects and a distribution for those effects. For a cluster \(j\), write \(Y_{ij}\) for the response of unit \(i\), \(b_j\) for the cluster’s random effect, and \(\theta\) for the model parameters. A full-model replicate is generated according to the fitted model: generate the random effects under their fitted distribution and then generate the responses conditional on those effects. The same simulated \(b_j\) is shared by observations in cluster \(j\), so the replicate retains the model’s within-cluster dependence.

The design must also be held fixed: use the same covariates, cluster membership, offsets, and observation pattern as in the data. If an outcome is a binomial count with a known number of trials, for example, retain that number of trials in each replicate. A collection of independent draws for each row that ignores shared random effects is not a valid replicate of a GLMM’s joint response structure.

This tutorial uses full-model predictive replicates to illustrate the construction. A diagnostic must always state its simulation target: a different target can answer a different question. The rank formula itself does not decide whether the random effects were generated or held fixed; the simulation scheme does.

Randomized Ranks for Discrete Outcomes

Let \(L_i\) be the number of simulated outcomes strictly less than \(y_i\), and let \(E_i\) be the number equal to \(y_i\). Introduce \(V_i\sim\operatorname{Uniform}(0,1)\), independently for the randomization, and define

The \(E_i+1\) counts the observed outcome as well as all simulated ties. Thus \(U_i\) is randomized across the whole rank interval occupied by the observed value and its ties. This prevents a discrete outcome with many repeated values from producing only a few possible residual values. The corresponding normal-score residual is

where \(\Phi\) is the standard normal cumulative distribution function. Values of \(U_i\) near \(0.5\) give normal scores near zero; small and large \(U_i\) give negative and positive scores, respectively. The sign describes whether the observation falls low or high relative to its simulated reference, not whether a covariate has a positive or negative effect.

The rank has finite resolution: it depends on \(B\), and randomization avoids exact endpoints. With small \(B\), the residuals are coarse and their extreme values are poorly resolved. Increasing \(B\) improves Monte Carlo resolution but does not repair a misspecified model or remove uncertainty caused by estimating its parameters.

Why the Rank Is Calibrated Under an Exact Model

Suppose the model and its parameters are known, and the observed response and the \(B\) simulated responses are exchangeable under the predictive distribution used for the comparison. Exchangeable means that permuting their labels does not change their joint distribution. If all values are distinct, each of the \(B+1\) possible ranks for the observed value is then equally likely. With ties, random ordering within the tied group gives the same result: the observed value’s randomized rank is equally likely to occupy any of the \(B+1\) rank positions.

Adding an independent uniform jitter within the rank’s interval makes \(U_i\) uniform on \((0,1)\) under this exact exchangeability condition. The tie formula above is equivalent to randomizing continuously across the tied rank interval. Applying \(\Phi^{-1}\) then gives a standard normal marginal distribution for \(R_i\). This is a calibration result under the stated assumptions; it does not imply that residuals from a fitted GLMM are independent of each other.

In practice, the model parameters are usually estimated from the same data being diagnosed. Conditional on the fitted estimates, the observed outcomes are not generally exchangeable with new outcomes generated as if those estimates were fixed and known. Consequently, plug-in simulation-based residuals are useful diagnostics, but their exact uniform or normal calibration is not guaranteed. Estimation, dependence, and Monte Carlo error should all be considered when interpreting their distribution.

Using Residuals and Checking Model Predictions

Simulation-based residuals are useful because the reference distribution can reflect a GLMM’s response distribution and dependence, rather than imposing a constant-variance normal approximation. They extend the idea behind randomized quantile residuals: locate an observation relative to its fitted distribution, with randomization for discrete outcomes. The practical difference is that the reference distribution can be constructed from simulated data, which is especially useful when the GLMM’s predictive distribution is not available in a convenient closed form.

Individual residuals are only one diagnostic view. A useful complementary check compares a data summary with the same summary calculated for every complete replicate data set. The summary should target a specific concern, such as the total number of zeros or the maximum cluster count. Because the replicate data sets retain the fitted clustering, their summary distribution reflects the dependence implied by the model.

Assumptions, Diagnostics, and Limitations

A residual’s apparent calibration depends on how well the simulations reproduce the intended predictive distribution. Before interpreting a pattern, check that the response family, link, random-effects structure, covariates, offsets, and observation design used to simulate are the ones specified by the fitted model. Verify that shared random effects are shared in each replicate and that all observations are generated together when the diagnostic targets the joint response. If replicates are generated independently by row despite cluster dependence in the fitted GLMM, the reference distribution is wrong.

Also check Monte Carlo adequacy. Repeat the simulation with more replicates when ranks or tail proportions are unstable, and recognize that an estimate based on few simulations has substantial Monte Carlo variation. A visibly uneven residual distribution can signal a predictive mismatch, but it can also reflect finite simulation error, parameter estimation, or dependence among residuals. Likewise, a seemingly well-behaved set of ranks does not prove that every part of a GLMM is correct: a diagnostic summary can be insensitive to some misspecifications.

Compared with response, Pearson, and deviance residuals, simulation-based ranks provide a common reference scale without relying on a normal approximation to the response. Compared with analytic randomized quantile residuals, they use an empirical reference generated from the fitted model. Their cost is computational: simulation is required, rank resolution is limited by \(B\), and uncertainty from model fitting is not automatically captured by fixed-parameter simulation. These residuals complement, rather than replace, checks of model structure and targeted predictive summaries.