Residual and Predictive Checks for Count PK/PD Models

Why Count-Model Checks Need More Than One View

A count model can reproduce the overall average response yet miss important features of the data. For example, it may predict too few zero responses at low concentration, underpredict large counts late in an observation period, or show systematic errors at particular concentrations. A single summary of fit can hide these patterns. Count-model diagnostics therefore examine the full predicted distribution and ask whether its predictions agree with observed counts across relevant conditions.

The checks in this tutorial apply to fitted Poisson, negative-binomial, zero-inflated, hurdle, and count models with subject-level random effects. They do not determine which model is universally best. Instead, they help identify where a model's predicted distribution is inconsistent with the observed data and whether that inconsistency depends on concentration, time, or repeated observations.

Definition: For observation $i$, let $Y_i$ be the observed count and let $F_i(k)=P(Y_i\leq k\mid \text{model and covariates for }i)$ be its fitted conditional cumulative distribution function. Its probability mass function is $f_i(k)=P(Y_i=k\mid \text{model and covariates for }i)$. The conditioning must match the intended prediction: for example, it may condition on a subject's estimated random effect or integrate over the random-effects distribution.

Before interpreting diagnostics, check that the model's response support and prediction target are appropriate. Counts are nonnegative integers; exposure duration or observation windows must be represented consistently; and concentration and observation time must be evaluated at the values relevant to each record. For repeated observations, the model should represent the dependence induced by subject heterogeneity. Treating correlated records as independent can make discrepancies appear more precise than they are.

Randomized Quantile Residuals

Ordinary residuals such as observed count minus fitted mean are difficult to compare across records when the outcome is discrete and the fitted means differ. A count of two, for example, may be unsurprising under one fitted distribution and unusual under another. Randomized quantile residuals use each record's full fitted distribution, then map the resulting probability scale to a common normal scale.

Definition: Let $V_i$ be independently drawn from a uniform distribution on $(0,1)$, independently of the data. The randomized probability integral transform and randomized quantile residual are $$ U_i=F_i(Y_i-1)+V_i f_i(Y_i), \qquad r_i=\Phi^{-1}(U_i), $$ where $F_i(-1)=0$ for a count outcome and $\Phi$ is the standard normal cumulative distribution function. If the fitted conditional distribution is the true distribution and its parameters are known, then $U_i$ is uniform on $(0,1)$, and $r_i$ has a standard normal distribution.

The randomization is needed because a discrete cumulative distribution function jumps at each possible count. For a realized count $Y_i=k$, the residual's probability value is drawn between $F_i(k-1)$ and $F_i(k)$, rather than being fixed at one of the jump endpoints.

Why the Randomized Transform Works

The justification follows from the probability mass assigned to each count. Conditional on $Y_i=k$, $U_i$ is uniform on the interval from $F_i(k-1)$ to $F_i(k)$. That interval has length $$ F_i(k)-F_i(k-1)=f_i(k). $$ For any subinterval inside it, the probability that $U_i$ falls in that subinterval equals its length: the probability of $Y_i=k$ is $f_i(k)$, and the conditional uniform density on the interval is $1/f_i(k)$. Their product is one. As $k$ ranges over the nonnegative integers, these adjacent intervals cover $(0,1)$ without gaps or overlap, so $U_i$ is uniform. Applying the inverse standard normal CDF to a uniform variable gives a standard normal variable.

This result is exact for a known, correct distribution. In practice, model parameters are estimated, and residuals may share information through fitted parameters or subject random effects. Thus a normal-looking residual histogram is not proof of model adequacy, and standard-normal reference lines are diagnostic guides rather than exact guarantees.

Worked Example: Computing a Randomized Residual

Suppose a hypothetical count observation has fitted Poisson mean $1.2$, and the observed count is $Y=1$. For this illustration, take the randomization draw to be $V=0.5$. The Poisson probabilities give $$ F(0)=e^{-1.2}=0.3012,\qquad P(Y=1)=1.2e^{-1.2}=0.3614. $$

The randomized probability is $$ U=F(0)+0.5P(Y=1)=0.3012+0.5(0.3614)=0.4819. $$ Therefore $$ r=\Phi^{-1}(0.4819)\approx -0.0454. $$

This residual is close to zero, so this particular count is near the center of its fitted distribution on the normal scale. It does not establish that the model fits other observations or that the model is adequate overall. A residual plot against concentration and observation time can reveal systematic patterns that a single residual cannot.

Useful displays include a residual histogram with a standard normal reference curve and residual plots against fitted mean, concentration, observation time, and subject. Look for skewness, heavy tails, clusters, and trends. A residual pattern that changes with concentration may indicate an inadequate concentration effect; a pattern that changes over time may indicate that the time course or exposure window is not represented adequately. With repeated data, plots should also be interpreted in light of subject-level dependence.

Rootograms and Zero-Frequency Checks

A rootogram compares observed and model-expected frequencies at possible count values. For each count $k$, the observed frequency is the number of records with $Y_i=k$; the expected frequency is the sum of their fitted probabilities, $\sum_i P(Y_i=k)$. Taking square roots makes differences more visible when frequencies vary in size. A hanging rootogram displays the observed and expected square-root frequencies from a shared baseline, making overprediction and underprediction easier to see.

Worked Example: Reading a Rootogram and Checking Zeros

In a hypothetical set of 100 records, suppose every fitted count distribution is Poisson with mean $1.2$. The observed frequencies are 40 zeros, 27 ones, 18 twos, and 15 counts of three or more. The fitted probabilities are $P(0)=0.3012$, $P(1)=0.3614$, and $P(2)=0.2169$. The probability of three or more is $1-F(2)=1-0.8795=0.1205$.

The expected frequencies are therefore $$ E_0=100(0.3012)=30.12,\quad E_1=100(0.3614)=36.14,\quad E_2=100(0.2169)=21.69,\quad E_{3+}=100(0.1205)=12.05. $$ They sum to 100, as they should. For the zero category, the square-root difference (observed minus expected) is $$ \sqrt{40}-\sqrt{30.12}=6.325-5.488=0.837. $$ For counts of one, it is $\sqrt{27}-\sqrt{36.14}=5.196-6.012=-0.816$. Thus the display would show more zeros and fewer ones than predicted, with a smaller shortfall at two and a modest excess in the $3+$ group. This pattern merits investigation; it does not by itself identify the reason for the mismatch.

The zero count can also be summarized directly. Under independent records and a fixed, correct model, the zero indicator for record $i$ has probability $p_{0i}=P(Y_i=0)$. The expected number of zeros is $\sum_i p_{0i}$, and its variance is $\sum_i p_{0i}(1-p_{0i})$. In this example, treating the fitted mean as fixed gives $$ SE_0=\sqrt{100(0.3012)(1-0.3012)}=\sqrt{21.05}=4.588, \qquad z_0=\frac{40-30.12}{4.588}=2.153. $$

This standardized discrepancy flags more zeros than expected under those simplifying assumptions. It is not automatically a calibrated hypothesis test: parameter estimation, repeated-measure dependence, or uncertainty in random effects can change the variability. A predictive simulation that refits the model and preserves the data structure is generally more appropriate for formal calibration.

A rootogram is useful because it can reveal problems beyond an excess of zeros, including deficits at intermediate counts or a poorly represented upper tail. A zero-frequency check focuses attention on a clinically or scientifically important feature, but it cannot diagnose the rest of the distribution. Sparse count categories also make rootograms unstable; combining a tail category such as $3+$ can help display it, but the combined category should be interpreted as a tail summary, not as a check of each tail count separately.

Predictive Intervals by Concentration and Time

A predictive interval asks whether an observed count is plausible for a new record under the fitted distribution. For a discrete outcome, define $q_a$ as the smallest nonnegative integer $k$ for which $F(k)\geq a$. An equal-tailed 90% predictive interval is $[q_{0.05},q_{0.95}]$. Because the distribution is discrete, its actual coverage can exceed 90%; the interval is not generally symmetric around the fitted mean.

Assessing intervals only in the full data set can conceal conditional misfit. Instead, compare observed counts with intervals in strata defined by concentration and observation time, or plot prediction coverage against those variables. Strata must be chosen with enough records to make the comparison interpretable. If predictions concern a new subject, integrate over the random-effects distribution; if they concern another record for an already observed subject, condition on the subject information as appropriate.

Worked Example: Stratified Predictive Intervals

Consider a hypothetical Poisson model evaluated in four concentration-by-time strata. For illustration, assume six independent records per stratum and use the fitted conditional means shown below. These intervals describe individual counts, not uncertainty in the fitted mean.

Stratum	Fitted mean	90% predictive interval	Hypothetical observed counts	Within interval
Low concentration, early	0.5	[0, 2]	0, 0, 1, 1, 2, 3	5 of 6
High concentration, early	2.0	[0, 5]	0, 1, 2, 3, 4, 6	5 of 6
Low concentration, late	1.2	[0, 3]	0, 1, 1, 2, 3, 4	5 of 6
High concentration, late	3.0	[1, 6]	1, 2, 3, 4, 6, 8	5 of 6

For example, with mean $0.5$, the Poisson CDF at 2 is $e^{-0.5}(1+0.5+0.125)=0.9856$, while $F(1)=0.9098$. Thus $q_{0.95}=2$, and $q_{0.05}=0$, giving $[0,2]$. With mean $2.0$, $F(4)=0.9473$ and $F(5)=0.9834$, so $q_{0.95}=5$; since $F(0)=0.1353$, $q_{0.05}=0$, giving $[0,5]$. The same quantile rule gives $[0,3]$ for mean $1.2$ and $[1,6]$ for mean $3.0$.

Each stratum has five of six observations inside its interval, but six records are far too few to establish the intended long-run coverage. The counts outside the intervals (3, 6, 4, and 8) are worth examining in context, particularly the high-concentration, late observation. A larger discrepancy concentrated in one stratum could indicate conditional misfit that overall coverage would miss.

For repeated counts, individual intervals alone do not assess whether the model captures the dependence between records from the same subject. A model can give reasonable record-level intervals while predicting the wrong within-subject clustering. In that setting, examine simulated subject trajectories or summaries of within-subject counts, using predictions that preserve the distinction between an existing subject and a new subject.

Choosing and Interpreting the Checks

Check	What it highlights	Main limitation
Randomized quantile residuals	Record-level distributional patterns on a common scale	Randomization adds visual noise; fitted residuals are not exactly independent standard normal values
Rootogram	Observed-versus-expected frequency differences across counts	Sparse categories and grouped tails can hide detailed discrepancies
Zero-frequency check	Whether the overall or stratum-specific number of zeros is surprising	It checks only one feature of the distribution
Stratified predictive intervals	Whether prediction coverage changes with concentration or time	Coverage is imprecise in small strata and depends on the prediction target

These approaches complement one another rather than compete. Residuals provide a record-level view; rootograms show the shape of count-frequency discrepancies; zero checks isolate the zero probability; and stratified intervals evaluate prediction in clinically relevant conditions. A global fit statistic or a visually satisfactory residual histogram cannot replace checks of concentration, time, and subject structure. Conversely, a single unusual count or a small stratum with one interval miss is not enough to conclude that the model fails.

Common Pitfall: Treating a Diagnostic as a Verdict

A common mistake is to see several observed zeros above the expected number and conclude that the data require a zero-inflated model. That conclusion is not justified by the zero count alone. The excess could reflect a poorly modeled concentration or time effect, unaccounted subject heterogeneity, incorrect exposure windows, or chance variation. The correction is to inspect several diagnostics, check the model's conditioning and prediction target, and assess whether the discrepancy persists across relevant strata. These checks identify patterns that need explanation; they do not, by themselves, establish a particular replacement model.

Another pitfall is interpreting a 90% predictive interval as a range containing the true model parameter with 90% probability. It is instead a range intended to contain a future count with the stated coverage under the model and prediction conditions. Parameter uncertainty and estimated random effects affect how that coverage should be assessed. When the model is fitted to the same data being checked, simulation-based reference distributions can account for fitting and dependence more realistically than treating all diagnostic quantities as fixed.

Key takeaway: Evaluate a count model's entire conditional distribution, not only its mean. Use randomized quantile residuals, observed-versus-expected count frequencies, a targeted zero check, and predictive assessments across concentration and time; interpret each in light of parameter estimation and repeated-subject structure.

Check Your Understanding

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

Why does a randomized quantile residual use a uniform draw between $F_i(Y_i-1)$ and $F_i(Y_i)$, rather than simply using $F_i(Y_i)$?
For a record with fitted $P(Y=0)=0.25$, $P(Y=1)=0.40$, observed $Y=1$, and $V=0.5$, calculate $U$. Is the residual above or below zero?
A rootogram shows many more zeros than expected and fewer counts of one, but the expected and observed total counts are similar. What can be concluded, and what cannot be concluded?
Why should predictive interval coverage be examined by concentration and observation time rather than only as one overall percentage?

Solutions. (1) The fitted CDF has jumps for discrete counts. Randomizing within the jump assigned to the observed count spreads probability uniformly across $(0,1)$ when the distribution is correct; using $F_i(Y_i)$ alone produces a discrete, non-uniform set of values. (2) $F(0)=0.25$, so $U=0.25+0.5(0.40)=0.45$. Since $0.45<0.5$, $\Phi^{-1}(0.45)$ is below zero. (3) The fitted distribution may misallocate probability among count categories despite a similar overall count total. The pattern flags a distributional discrepancy but does not establish its cause or prove that a particular alternative model is needed. (4) Overall coverage can conceal overprediction in one range and underprediction in another; stratifying reveals whether predictive performance changes with exposure conditions.