Why Randomize a Discrete Residual?
For a continuous outcome with a correctly specified distribution, its cumulative probability evaluated at the observation is uniform on \((0,1)\). Transforming that probability through the inverse standard normal CDF gives a standard normal value. A discrete outcome behaves differently: its cumulative distribution function (CDF) jumps, so evaluating the CDF at an observation produces only a limited set of possible values. Ties and gaps are unavoidable, and a normal-quantile transformation of those values is not generally normally distributed.
Randomized quantile residuals address this discreteness by selecting a random value inside the probability interval associated with the observed outcome. The transformation is useful in a GLMM because it puts outcomes from discrete response families on a common, approximately normal scale. Its diagnostic value depends on using the fitted outcome distribution appropriate to the question, including a clear choice about how random effects enter that distribution.
Definition and Probability-Integral-Transform Justification
Let \(Y_i\) be a discrete response with fitted CDF \(\hat F_i(y)\), and let the observed value be \(y_i\). Write \(\hat F_i(y_i^-)=P_{\hat F_i}(Y_i<y_i)\), the fitted probability accumulated strictly below the observation. The probability assigned to the observed value is \(\hat F_i(y_i)-\hat F_i(y_i^-)\).
For a GLMM, the fitted CDF must be specified. With cluster \(j\), observation \(i\), random effect \(b_j\), and fitted parameters, a conditional choice is \(\hat F_{ij}^{\,c}(y\mid\hat b_j)\), the fitted response CDF given the estimated random effect. A marginal choice is the CDF after integrating over the random-effects distribution. These answer different calibration questions; the CDF used to construct the residual should be stated.
The reason randomization works can be seen first when the distribution is known. Suppose its support values are ordered \(y_1,y_2,\ldots\), with \(c_k=P(Y\le y_k)\) and \(c_0=0\). If \(Y=y_k\), draw \(U\) uniformly on \((c_{k-1},c_k)\). For any \(t\) within that interval,
The same argument applies in each probability interval, so \(U\) is uniform on \((0,1)\). Therefore \(\Phi^{-1}(U)\) has a standard normal distribution. For a GLMM, this exact result holds for the conditional responses when the conditional distribution and its parameters are correct and known, conditional independence holds, and the conditional CDF is used. In practice parameters and random effects are estimated, so residuals are diagnostic approximations, not an exact standard-normal sample.
Worked Calculations
Worked Example: A Bernoulli Outcome
Suppose a hypothetical GLMM gives a fitted conditional probability of success \(0.20\) for one observation, conditional on its covariates and estimated random effect. Thus \(Y\) is Bernoulli with fitted success probability \(0.20\). If the observed outcome is \(0\), the fitted probability below it is \(0\), and the fitted CDF at it is \(P(Y=0)=0.80\). Suppose the random draw for this residual is \(U=0.34\).
The residual is \(r=\Phi^{-1}(0.34)\approx-0.412\). This is a moderately negative residual: the randomized probability lies below \(0.50\), so it maps below zero on the normal scale.
If the outcome instead were \(1\), its probability interval would be \((0.80,1)\). For a draw \(U=0.91\), the residual would be \(r=\Phi^{-1}(0.91)\approx1.341\). These calculations show why a discrete observation does not determine one residual by itself: the random draw within its probability interval matters.
Worked Example: A Binomial Outcome
For a hypothetical observation with four trials, suppose the fitted conditional success probability is \(0.30\), and \(Y=2\) successes were observed. Assume the fitted binomial model is the intended conditional response distribution. The probabilities for \(Y=0,1,2\) are \(0.7^4=0.2401\), \(4(0.3)(0.7^3)=0.4116\), and \(\binom{4}{2}(0.3^2)(0.7^2)=0.2646\).
Consequently, the fitted probability below the observed value is \(\hat F(1)=0.2401+0.4116=0.6517\), and the CDF at the observation is \(\hat F(2)=0.6517+0.2646=0.9163\). If the randomized draw is \(U=0.80\), it lies within \((0.6517,0.9163)\). The residual is \(r=\Phi^{-1}(0.80)\approx0.842\).
This is a positive residual, but it is not a raw difference of \(2-4(0.30)\). Its scale comes from the fitted probability assigned to all possible binomial outcomes and the selected position within the observed outcome's probability interval.
Worked Example: A Sparse Poisson Outcome
Suppose a hypothetical count observation has a fitted conditional Poisson mean of \(1\), given its covariates and estimated random effect. For a Poisson variable with mean \(1\), \(P(Y=k)=e^{-1}/k!\). If the observed count is \(4\), then
A randomized draw \(U=0.989\) lies in that interval, giving \(r=\Phi^{-1}(0.989)\approx2.290\). The count is in the upper tail of this fitted distribution and produces a large positive residual for this particular draw.
Now consider an observed count of \(0\), whose interval is \((0,e^{-1})\), or approximately \((0,0.3679)\). A draw \(U=0.12\) gives \(r\approx-1.175\); a draw \(U=0.31\) gives \(r\approx-0.496\). The observed zero is unchanged, but the residual differs because each draw selects a different location within the interval. A single large residual is a signal to investigate, not proof that the model is wrong.
Choosing the GLMM Distribution and Reading the Residuals
In a GLMM, a response distribution can be considered conditional on its random effect or marginalized over the random-effects distribution. With a conditional residual, the relevant interval is built from \(F_{ij}(y\mid b_j)\) evaluated at the fitted random effect. Under a correctly specified conditional model with known parameters and random effects, and conditionally independent responses, these residuals are standard normal and conditionally independent. Estimated random effects and parameters weaken that exact reference: residuals can be dependent within a cluster, and their distribution can be affected by estimation.
Alternatively, define the marginal CDF by integrating the conditional CDF over the random-effects distribution. This gives an unconditional fitted probability for an individual response. Under a correctly specified marginal distribution, each randomized residual has a standard normal marginal distribution when the parameters are known. But shared random effects can still induce dependence between residuals from the same cluster. Marginal standard-normal calibration does not imply that the whole set of residuals is independent.
A useful first check is whether the residual distribution is broadly compatible with a standard normal reference: centered near zero, with a spread near one and no conspicuous excess of extreme values. These are approximate expectations, not exact sample requirements. For a finite fitted sample, random variation and parameter estimation can create visible differences. In addition, patterns in residuals associated with covariates, fitted probabilities, or clusters may reveal a problem that an overall histogram obscures.
State whether the residual uses the fitted response distribution conditional on estimated random effects or a distribution marginalized over random effects.
For each observation, compute the fitted probability strictly below its observed value and the fitted CDF at that value.
Draw one uniform value between the two bounds and apply the inverse standard normal CDF.
Look for broad departures from normality and structured patterns, while recognizing estimation effects, within-cluster dependence, and randomization variability.
Comparisons, Diagnostics, and Limitations
Response, Pearson, and deviance residuals summarize discrepancy on scales tied to the response or its fitted variance and deviance contribution. Randomized quantile residuals instead use the full fitted discrete CDF and map its probability scale to normal quantiles. This makes their scale comparable across discrete response families, but does not make them universally superior: probability calculations may be difficult for complicated response distributions, and the normal reference is only approximate after fitting.
Randomness is an important practical feature, not a nuisance to ignore. A fixed random seed makes a diagnostic reproducible, but one seed can give a misleading impression when many observations have wide probability intervals. Repeating the randomization can show whether a suspected pattern is stable. Randomization does not repair a misspecified mean, link, dispersion, random-effects distribution, or dependence structure; it only spreads discrete probability mass across its corresponding CDF interval.
Other warning signs include an observed value assigned zero fitted probability, which makes the interval empty and the residual undefined; numerical probabilities that fail to sum accurately; and extreme-tail probabilities rounded to exactly zero or one before the inverse-normal transformation. In those cases, check the fitted distribution and compute tail probabilities with adequate numerical precision. An outcome assigned zero probability is especially consequential: the fitted model cannot generate that observation under the stated distribution.
A common mistake is to claim that randomized residuals must be exactly independent standard normal values whenever the GLMM is correct. That statement overlooks estimated parameters, estimated random effects, and dependence induced by shared cluster effects. The correction is to treat normality as a calibration diagnostic under a specified conditional or marginal distribution, interpret it alongside dependence and structure, and avoid reading a single randomization as definitive evidence.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A fitted discrete CDF gives \(F(y^-)=0.25\) and \(F(y)=0.55\). Is \(U=0.60\) a valid randomized value for this observation? If \(U=0.40\), what is the residual's sign?
- For a Bernoulli outcome with fitted success probability \(0.10\), an observed zero has what probability interval? If \(U=0.50\), is its residual positive, negative, or zero?
- Why can marginal randomized quantile residuals have standard normal marginal distributions yet remain dependent within clusters?
- What does a large positive residual suggest, and why does it not by itself establish model failure?
Solutions. (1) No: \(0.60\) is outside \((0.25,0.55)\). The value \(0.40\) is inside the interval and gives \(\Phi^{-1}(0.40)\approx-0.253\), a negative residual. (2) The zero interval is \((0,0.90)\), since \(P(Y=0)=0.90\). With \(U=0.50\), the residual is \(\Phi^{-1}(0.50)=0\). (3) A marginal CDF can calibrate each response separately, while a shared random effect can associate responses within a cluster; standard normal marginal distributions do not guarantee joint independence. (4) It places the randomized probability in the upper tail relative to the fitted distribution. One such value can arise by chance or from the random draw, so it motivates broader checks rather than proving misspecification.