Residual Diagnostics for Binary, Ordinal, and Count PD

Why Discrete PD Outcomes Need Specialized Residuals

A pharmacodynamic endpoint may be a binary response, an ordinal severity grade, or a count of events. For these outcomes, the difference between the observed value and its fitted mean is not distributed like a continuous, symmetric error. For example, a binary observation can only be 0 or 1, and a count cannot be negative. A residual defined as “observed minus predicted” can still be useful, but its scale and spread depend on the predicted probability or count. Discrete outcomes therefore need diagnostics that respect their full predicted distributions.

The earlier material on observed-predicted calibration and residual patterns provides useful starting points: compare predictions with outcomes on the response scale, and look for structure against time, exposure, and covariates. Here the focus is on diagnostics that use the probability assigned to the observed category, rather than treating the outcome as continuous. These checks assess calibration, dispersion, and remaining patterns; they do not identify a unique cause of misfit.

Definition: Let $Y_i$ be a discrete PD observation for occasion $i$, and let $F_i(y)=P(Y_i\le y\mid \text{model and observed predictors})$ be its fitted conditional cumulative distribution function. Write $F_i(y-)=P(Y_i<y\mid \text{model and observed predictors})$. A randomized quantile residual is $$ U_i=F_i(Y_i-)+V_i\{F_i(Y_i)-F_i(Y_i-)\},\qquad r_i=\Phi^{-1}(U_i), $$ where $V_i\sim\mathrm{Uniform}(0,1)$, independently randomized for each observation, and $\Phi^{-1}$ is the standard normal quantile function. The randomization spreads the probability mass at the observed outcome across its CDF jump.

Why the Randomization Works—and What It Assumes

If the model distribution is correct and its parameters are known, the randomized value $U_i$ is uniform on $(0,1)$. To see why, consider any possible outcome $y$ with probability $p_y=F_i(y)-F_i(y-)$. Conditional on $Y_i=y$, $U_i$ is uniform on the interval from $F_i(y-)$ to $F_i(y)$, whose width is $p_y$. Its unconditional density on that interval is therefore $p_y/p_y=1$. The intervals for all possible outcomes partition $(0,1)$, so the density is one throughout. Applying the standard normal quantile function then gives a standard normal marginal distribution for $r_i$.

This result relies on the fitted distribution being the correct distribution for the observation, including its conditioning structure. For repeated observations, the assumed dependence among observations from the same subject also matters. If the model includes subject-level random effects, predictions and simulations must handle those effects consistently—for example, by conditioning on estimated subject effects for a conditional check, or integrating over them for a population-level check. Drawing each occasion independently when the model implies within-subject dependence can create misleading diagnostics.

In practice, model parameters are estimated from the same data, so fitted residuals are not exactly independent standard normal variables even when the model is adequate. They may also change with the random number used for $V_i$. A residual plot or histogram is therefore evidence to investigate, not an exact normality test. Repeating the randomization can show whether a visual feature is stable or is mainly due to the added random variation.

Calculating Residuals for Different Outcome Types

The same definition applies to binary, ordinal, and count outcomes, but their fitted CDFs are constructed differently. For a binary response $Y_i\in\{0,1\}$ with $P(Y_i=1)=p_i$, the probability mass at 0 is $1-p_i$ and at 1 is $p_i$. For an ordinal response, the model assigns a probability to each ordered category; the cumulative probabilities are sums of those category probabilities. For a count response, the CDF comes from the specified count distribution, such as a Poisson or negative binomial distribution.

For a discrete observation, the residual does not simply measure distance from the fitted mean. It locates the observed outcome within the probability mass the model assigns to it. If an outcome is rare under the model, its CDF interval is narrow; if common, its interval is wider. The random draw chooses a location within that interval. Consequently, individual randomized residuals should not be overinterpreted as fixed scores.

Worked Example: A Binary PD Response

Suppose a fitted model gives a hypothetical probability $p_i=0.30$ of a binary PD response at one occasion, and the observed response is $Y_i=1$. For the event $Y_i=1$, $F_i(1-)=P(Y_i=0)=0.70$, while $F_i(1)=1$. Use a hypothetical random draw $V_i=0.40$.

The randomized probability is $$ U_i=0.70+0.40(1-0.70)=0.70+0.12=0.82. $$ Thus $$ r_i=\Phi^{-1}(0.82)\approx 0.915. $$

This positive residual says that this randomized location within the model’s probability interval is above the standard normal median. It does not mean the response is “0.915 units” above its prediction, nor does one such value establish model misfit. Across observations, calibration would be questionable if responses of 1 occurred systematically more often than the model’s predicted probabilities, or if residual patterns persisted across time or exposure.

Worked Example: An Ordinal PD Grade

Suppose a hypothetical ordinal endpoint has four ordered grades, with fitted probabilities $0.10, 0.20, 0.50,$ and $0.20$. The cumulative probabilities at grades 1, 2, and 3 are $0.10, 0.30,$ and $0.80$. For an observation in grade 3, the CDF interval is $(0.30,0.80]$. With $V_i=0.60$, $$ U_i=0.30+0.60(0.80-0.30)=0.60, $$ so $$ r_i=\Phi^{-1}(0.60)\approx 0.253. $$

The fitted probabilities also give a direct calibration summary for a hypothetical group of 20 comparable observations. Their expected grade counts are $20(0.10)=2$, $20(0.20)=4$, $20(0.50)=10$, and $20(0.20)=4$. If the observed counts are 2, 5, 9, and 4, the observed-minus-expected counts are $0,1,-1,0$. This small discrepancy is descriptive, not by itself a formal test. Larger or recurring discrepancies—especially if concentrated at severe grades or at particular exposures—would prompt investigation of the category probabilities or covariate structure.

Worked Example: A Poisson Count PD Outcome

Suppose a hypothetical count response follows a Poisson distribution with fitted mean $\lambda_i=1.5$, and the observed count is 3. The Poisson probability of count $k$ is $P(Y_i=k)=e^{-\lambda_i}\lambda_i^k/k!$. The probabilities through count 2 are approximately $P(0)=0.2231302$, $P(1)=0.3346952$, and $P(2)=0.2510214$, giving $$ F_i(2)=P(Y_i\le 2)=0.2231302+0.3346952+0.2510214=0.8088468. $$

For count 3, $$ P(Y_i=3)=e^{-1.5}\frac{1.5^3}{3!}\approx 0.1255107, \qquad F_i(3)=P(Y_i\le 3)=0.9343575\approx 0.9344. $$ The cumulative probability is computed from the unrounded probabilities; adding displayed, rounded components can give a slightly different last digit. Taking $V_i=0.50$, the randomized value is $$ U_i=F_i(2)+0.50\{F_i(3)-F_i(2)\} =0.8088468+0.50(0.1255107) \approx 0.8716022. $$ Therefore $r_i=\Phi^{-1}(0.8716022)\approx 1.13$. This is one randomized residual for a count of 3 under a Poisson mean of 1.5. A broader diagnostic should ask whether counts across the data are more variable than the model predicts, whether zeros or large counts are systematically misrepresented, and whether discrepancies vary with time or exposure.

Simulation-Based Checks

Simulation-based diagnostics compare observed features of the data with features expected under the fitted model. They are especially useful when discrete residuals have awkward finite-sample distributions, parameters have been estimated, or the data are repeated within subjects. A simulation should reproduce the design and dependence structure relevant to the question: observed covariates and sampling times are generally held fixed, while outcomes are generated from the fitted model. If uncertainty from parameter estimation is important, each simulated data set can be refitted before its diagnostic summaries are calculated.

One approach is to simulate $B$ replicated outcomes $Y_{ib}^{\mathrm{rep}}$ for each observation, then compare the observed outcome with its simulated distribution. A randomized rank can be formed from the number of replicates below the observation plus a random allocation among ties. Under an adequate model, these ranks should not show systematic concentration near the low or high ends. Another approach is to calculate a discrepancy—such as the proportion of zeros, a Pearson-type dispersion summary, or the difference between early- and late-time residual means—for the observed data and for every simulated data set. The observed discrepancy is concerning when it is unusual relative to the simulated reference distribution.

Worked Example: Testing for a Remaining Time Pattern by Simulation

Consider a hypothetical repeated binary PD data set. After fitting a model, define $T$ as the mean randomized residual at late observation times minus the mean at early times. Suppose the observed value is $T_{\mathrm{obs}}=0.82$. Simulate 199 data sets under the fitted model, preserving each subject’s observation times and the model’s within-subject structure. For each replicate, refit the model and calculate the same $T$. Suppose 4 of the 199 simulated data sets have $|T|\ge 0.82$.

A simulation-based two-sided tail estimate, using the usual plus-one adjustment, is $$ \hat{p}_{\mathrm{sim}}=\frac{1+\#\{|T_b|\ge |T_{\mathrm{obs}}|\}}{B+1} =\frac{1+4}{199+1}=0.025. $$ This indicates that a discrepancy at least as large in absolute value as the observed early-late difference was uncommon in these fitted-model simulations. It is a diagnostic warning about a possible remaining time pattern, not proof that the model is false or that time itself causes the pattern. The model, simulation scheme, and scientific context should be reviewed before changing the structure.

Reading Calibration, Dispersion, and Residual Patterns

Use several complementary summaries rather than relying on a single plot or statistic. On the response scale, compare observed proportions, category frequencies, or count summaries with their model-based expectations, preferably across clinically relevant exposure, time, baseline, or covariate strata. Randomized quantile residual histograms and normal quantile plots can reveal skewness, heavy tails, or clustering, but apparent departures should be checked across repeated randomizations or against simulations from the fitted model.

For counts, a Pearson-type dispersion summary is $$ D_P=\sum_i\frac{(Y_i-\mu_i)^2}{\operatorname{Var}(Y_i\mid\text{model})}, $$ where $\mu_i$ and the conditional variance come from the fitted count distribution. For a Poisson model, both equal $\lambda_i$. Under idealized independent, correctly specified settings, $D_P$ is often compared with the number of observations minus the number of estimated parameters, but that reference can be unreliable for repeated data, estimated random effects, or complex designs. Simulation of the full fitting procedure is safer than applying that simple reference mechanically.

Overdispersion means that variability exceeds what the fitted distribution predicts; underdispersion means it is smaller. Excess zeros, too many high counts, or a mismatch concentrated in particular ordinal categories can also indicate distributional misfit. A pattern in residuals against exposure may suggest that predicted outcome probabilities or means do not change adequately with exposure. A trend over time can point to an omitted time structure or a misspecified dynamic component. These are clues, not diagnoses: a pattern may also arise from dependence, an incorrect baseline, omitted covariates, or an unsuitable outcome distribution.

Practical diagnostic principle: A satisfactory diagnostic picture is not “all residuals are normal.” It is the absence of systematic, reproducible discrepancies from the model’s predicted outcome distributions, assessed on the response scale and through residual or simulation summaries that respect the design and dependence structure.

Common Pitfalls and Alternatives

Treating one randomized residual as fixed. The value depends partly on $V_i$. Repeat the randomization or use simulation-based summaries before treating a single unusual value as evidence.
Using continuous-residual intuition without adjustment. Discrete outcomes have ties and probability jumps. The randomized transform addresses those jumps, but fitted residuals are not exactly independent standard normal observations after parameter estimation.
Simulating occasions independently for repeated subjects. This destroys modeled dependence and can make predictive variability appear wrong. Simulate whole subject records using the model’s dependence structure.
Calling a residual pattern a causal effect. A trend with exposure or time flags a possible mismatch; it does not establish the cause or justify a specific structural change on its own.

Randomized quantile residuals provide a compact common scale for different discrete outcome types. Response-scale calibration plots retain direct clinical interpretation, while simulation-based checks can account for discreteness, fitting, and repeated sampling structure. None is universally best: use response-scale summaries to see what the model predicts, residual plots to locate patterns, and simulations to judge whether those patterns are unusual under the full fitted procedure. Each approach remains limited by the adequacy of the model and the realism of the simulation.

Check Your Understanding

Answer each question, then compare with the solutions below.

A binary model predicts response probability 0.20. For an observed response of 0, state the CDF interval used to construct its randomized quantile residual.
Why can two randomized quantile residuals for the same observation differ when the fitted model and observed outcome are unchanged?
For a Poisson model with mean $\lambda$, what feature of a count residual diagnostic would suggest overdispersion?
Why should a simulation-based check preserve the repeated-measures structure when the PD data contain multiple observations per subject?

Solutions. (1) For $Y=0$, $F(0-)=0$ and $F(0)=P(Y=0)=1-0.20=0.80$, so the randomized value is $U=0+V(0.80)$, with $V\sim\mathrm{Uniform}(0,1)$. (2) The residual depends on the random draw $V$, which selects a location within the CDF jump for the observed outcome. (3) Overdispersion is suggested when the observed count variability, summarized using squared deviations scaled by the model variances or assessed by simulation, is systematically larger than the Poisson model predicts. (4) The model predicts a joint pattern of outcomes within a subject; independent occasion-by-occasion simulation would remove that dependence and give an incorrect reference for variability and residual patterns.