Why Use Simulation-Based Residuals?
For a continuous Gaussian response, a residual can often be scaled by its standard deviation and compared with a normal distribution. That approach does not transfer directly to non-Gaussian PD outcomes. A binary response has only two possible values, an ordinal response has a finite set of categories, and a count response consists of integers. Their raw differences from predictions are therefore discrete, may have unequal variances, and need not be symmetric. A residual of one unit, for example, cannot have the same interpretation across these response types.
Simulation-based residuals address this by asking where the observed outcome falls in the outcome distribution generated by the fitted model. The distribution incorporates the response type and its predicted probabilities, rather than assuming that the difference between observation and fitted mean has a common Gaussian scale. The resulting residual is useful for checking whether the model's implied outcome distributions are broadly compatible with observed data.
The randomization spreads the probability mass associated with an observed category across the interval of cumulative probability assigned to that category. For a discrete response, \(F_i(Y_i^-)\) is the probability of outcomes strictly below the observation, while \(F_i(Y_i)-F_i(Y_i^-)\) is the probability of the observed value itself. Thus \(U_i\) lies between the lower and upper cumulative-probability limits for the observation. The normal quantile transformation expresses that location on a familiar scale; it does not make the response itself continuous or Gaussian.
Why the Transformation Works
The key result is about the distribution of \(U_i\) when the fitted conditional distribution is the true distribution for \(Y_i\), with parameters treated as known. For a discrete response with possible ordered values \(y_1,y_2,\ldots\), the cumulative probability intervals are \([F_i(y_j^-),F_i(y_j)]\). The length of the interval for \(y_j\) is exactly \(P(Y_i=y_j)\). Conditional on observing \(y_j\), randomization makes \(U_i\) uniform over that interval. Since the interval is selected with probability equal to its length, the combined distribution is uniform on \((0,1)\).
In the second sum there can be at most one interval containing \(u\). Its contribution is \(u-F_i(y_j^-)\), because the probability of observing \(y_j\) equals the interval length. The first sum contributes the lengths of all completed intervals, which add to \(F_i(y_j^-)\) for the interval containing \(u\). Together they give \(u\). Therefore \(U_i\sim\mathrm{Uniform}(0,1)\), and applying the increasing quantile function \(\Phi^{-1}\) gives \(R_i\sim N(0,1)\).
This result describes the marginal distribution of each residual under a correctly specified model. It does not guarantee that residuals from repeated observations on the same subject are independent. Nor does it ensure exact standard-normal residuals when parameters have been estimated from the same data. Those distinctions matter when interpreting residual plots and formal summaries.
Worked Examples
Worked Example: A Binary PD Response
Suppose a hypothetical model predicts that a binary response indicating whether a specified PD event occurs has probability \(p=0.30\) for one observation. The observed response is \(Y=1\). Compute a randomized residual using the hypothetical random draw \(V=0.50\).
For a Bernoulli outcome, \(P(Y=0)=0.70\) and \(P(Y=1)=0.30\). For the observed value \(Y=1\), \(F(1^-)=P(Y<1)=0.70\), and \(F(1)=1.00\). Substitution gives
The residual is \(R=\Phi^{-1}(0.85)\approx 1.036\). On the standard-normal residual scale, this realized randomization places the observed event above the median of its fitted outcome distribution. It does not mean the event was 1.036 response units above its prediction: the response remains binary. A different random draw would produce a different residual for the same observation, somewhere between \(\Phi^{-1}(0.70)\) and the upper limit approached as \(U\) approaches 1.
Worked Example: An Ordinal PD Response
Consider a hypothetical ordinal symptom response coded 0, 1, or 2, with fitted category probabilities \(P(Y=0)=0.20\), \(P(Y=1)=0.50\), and \(P(Y=2)=0.30\). The observed category is \(Y=1\). Using \(V=0.50\), calculate its residual.
The cumulative probability strictly below category 1 is \(F(1^-)=P(Y=0)=0.20\). The cumulative probability through category 1 is \(F(1)=0.20+0.50=0.70\). Therefore,
The residual is close to zero because this randomized location falls near the center of the model's probability interval for category 1. Its sign and size depend on the fitted probabilities and the randomization, not on treating the categories as equally spaced numerical measurements. The category codes establish order here; the calculation uses the model's category probabilities.
Worked Example: A Count PD Response
Suppose a hypothetical count-response model gives a Poisson distribution with fitted mean \(\lambda=2.5\) events for one observation. The observed count is \(Y=4\). Calculate the residual using \(V=0.50\).
For a Poisson distribution, \(P(Y=k)=e^{-\lambda}\lambda^k/k!\). Summing the probabilities through 3 gives \(F(4^-)=P(Y\leq3)\approx0.7576\). The probability at 4 is \(P(Y=4)=e^{-2.5}(2.5^4/4!)\approx0.1336\), so \(F(4)\approx0.8912\). Thus,
The positive residual indicates that the chosen randomized location for count 4 is above the median of the fitted count distribution. It is not a count difference of 0.931 events. If the variance of the count response is materially larger than its mean, a Poisson distribution may be inadequate; the residual calculation cannot repair a misspecified outcome distribution.
Using Simulations When the CDF Is Not Convenient
If the fitted outcome distribution is available analytically, calculate its cumulative probabilities directly as in the examples. In more complex models, it may be easier to generate \(M\) outcomes \(Y_i^{(1)},\ldots,Y_i^{(M)}\) from the model for the same design and covariates as observation \(i\). These simulations should represent the same predictive target as the residual: for example, either a distribution conditional on specified subject effects or one that integrates over the relevant between-subject variability. Mixing those targets can make a residual look unusual simply because the wrong distribution was used.
A simulation-based approximation can use the simulated outcomes to estimate the lower and upper cumulative-probability limits for the observed value, then randomize within the estimated jump. One rank-based implementation randomizes the observed value's position among the \(M+1\) exchangeable values consisting of the observation and its \(M\) simulations. Let \(L\) be the number of simulations strictly below the observation, \(E\) the number equal to it, and \(V\sim\mathrm{Uniform}(0,1)\). A jittered rank approximation is
This approximation relies on the observed value and simulations being exchangeable under the model being checked. With finite \(M\), the attainable ranks are limited and the estimated residual can vary with the simulation draw. More simulations improve resolution but do not remove model-parameter estimation error or dependence among observations. The calculation also requires avoiding exact endpoints \(0\) and \(1\), which correspond to infinite normal quantiles; continuous jitter gives endpoints probability zero in the ideal calculation.
State what is held fixed and whether relevant random effects or parameter uncertainty are included.
Use the model CDF directly, or estimate its location from simulations with ties handled by randomization.
Apply the standard-normal quantile function and assess distributional calibration, while retaining the response type and dependence structure in mind.
Diagnostics, Alternatives, and Limitations
If the model is adequate for the chosen predictive target and the observations are independent under that target, the residuals should have an approximately standard-normal marginal distribution. A histogram or quantile-quantile plot can reveal broad departures such as skewness, heavy tails, or an excess of extreme residuals. A few residuals can also be compared with their fitted probabilities or expected outcome distributions to identify observations for closer review. These summaries are diagnostic clues, not automatic tests that identify a unique model flaw.
Simulation-based residuals should be interpreted with several qualifications:
- Randomization adds Monte Carlo variation. Repeating the calculation with a different seed can change residuals, especially when the observed category has a narrow probability interval. Assess whether conclusions persist across randomizations or use enough simulations for adequate resolution.
- Estimated parameters affect calibration. The standard-normal result is exact when the conditional distribution is correct and its parameters are known. Fitting parameters to the same observations generally makes the residual distribution only approximately standard normal.
- Dependence is not removed by marginal transformation. Each residual can have a standard-normal marginal distribution while residuals within a subject remain correlated. Do not infer independence from a normal-looking histogram.
- Very coarse outcomes provide limited information per observation. A binary residual is randomized across one of two probability intervals; an apparently unusual value should be interpreted in light of the underlying event probability and the number of observations.
- A poor distributional fit can be hidden by a convenient summary. Overall residual normality does not establish that probabilities are calibrated across all relevant covariate or prediction ranges. Examine patterns in context and compare the simulated outcome distribution with the observed outcome features of interest.
Pearson residuals instead divide the difference between an observation and its fitted mean by a model-based standard deviation; they offer a direct standardized discrepancy but can remain skewed or discrete for non-Gaussian outcomes. Deviance residuals summarize a contribution to a likelihood-based discrepancy and can be useful for model comparisons, but their scale is not generally the same normal-quantile calibration described here. Predictive checks compare observed summaries or distributions with replicated data and retain the response scale; they are often easier to communicate but do not assign each observation a single transformed residual. These approaches answer related, not identical, diagnostic questions.
A common mistake is to say, “The residual is standard normal by construction, so a normal-looking residual plot proves the model fits.” The correction is that the transformation has an exact standard-normal marginal distribution only when the outcome distribution used in the transformation is correct and known. With fitted parameters, finite simulations, dependence, or a misspecified conditioning target, that guarantee does not hold. Even a good marginal residual distribution cannot rule out structured lack of fit elsewhere.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- A binary outcome has fitted event probability \(0.20\) and the observed outcome is 0. What are \(F(0^-)\) and \(F(0)\)?
- Why is randomization needed when an observed count has positive probability under a discrete model?
- For an ordinal outcome, why should category probabilities be used rather than subtracting the category code from its fitted mean?
- What does a standard-normal-looking residual histogram fail to establish?
Solutions. (1) For outcome 0, \(F(0^-)=0\) and \(F(0)=P(Y=0)=0.80\); hence the randomized \(U\) is between 0 and 0.80. (2) Without randomization, the discrete CDF has jumps, so its value at the observation is not continuously uniform; spreading each mass across its CDF interval produces a uniform value under the correct model. (3) Category codes indicate order, but their numerical distances need not represent equal changes in response; the fitted probability distribution respects the ordinal model. (4) It does not establish conditional calibration across covariates, correctness of the predictive distribution, independence of residuals, or adequate parameter estimation; marginal normality alone is not proof of model adequacy.