Randomized Quantile Residuals for Discrete Responses

Why Discrete Responses Need an Interval

For a continuous response, a fitted cumulative distribution function can locate an observation on the probability scale. For a discrete response, the CDF jumps at each possible outcome: an observed value corresponds not to one probability but to an interval of probabilities. Randomized quantile residuals use a random point inside that interval, then convert it to the normal scale. This provides a common way to calculate residuals for different discrete response distributions.

The key calculation is the size and location of the CDF jump at the observed value. Use the fitted distribution for observation $i$, not a single distribution shared across observations unless the model specifies one. The fitted CDF is denoted $F_i$, and its probability mass function is $f_i$.

Definition: For an observed discrete response $y_i$, define the lower and upper ends of its fitted CDF interval by $F_i(y_i-)=P(Y_i<y_i)$ and $F_i(y_i)=P(Y_i\leq y_i)$. Draw $V_i\sim\operatorname{Uniform}(0,1)$, independently of the response randomization for other observations, and set $$ U_i=F_i(y_i-)+V_i\{F_i(y_i)-F_i(y_i-)\}. $$ The randomized quantile residual is $r_i=\Phi^{-1}(U_i)$, where $\Phi$ is the standard normal CDF. Since $F_i(y_i)-F_i(y_i-)=f_i(y_i)$, this can also be written $U_i=F_i(y_i-)+V_i f_i(y_i)$.

The notation $F_i(y_i-)$ means the probability strictly below the observed response; it is not the probability at the observation. The upper endpoint includes the probability at the observation. Thus the interval has width exactly equal to the fitted probability mass $f_i(y_i)$. If the response can take only nonnegative integer values, for example, $F_i(y_i-)=P(Y_i\leq y_i-1)$.

Result: The interval spans exactly the observed outcome's probability mass. For any discrete response and observed value $y_i$, $$ F_i(y_i)-F_i(y_i-)=P(Y_i=y_i)=f_i(y_i). $$ Proof: The event $\{Y_i\leq y_i\}$ is the disjoint union of $\{Y_i<y_i\}$ and $\{Y_i=y_i\}$. Additivity of probability therefore gives $P(Y_i\leq y_i)=P(Y_i<y_i)+P(Y_i=y_i)$. Subtracting $P(Y_i<y_i)$ from both sides gives the stated identity. Consequently $U_i$ lies between the two fitted cumulative probabilities, with the random fraction $V_i$ determining its position within the jump. $\square$

Calculation Procedure

For a particular observation, identify its fitted distribution and parameter values first. Then calculate the probability strictly below the observed value, the probability at that value, and the randomized point within the resulting interval. Apply the standard normal inverse CDF only after obtaining $U_i$. This order makes it less likely to substitute a mass for a cumulative probability or to use the wrong endpoint.

Find the fitted mass at the observation.
Calculate $f_i(y_i)$ under the fitted model for observation $i$.

Find the lower cumulative probability.
Calculate $F_i(y_i-)=P(Y_i<y_i)$, using a strict inequality.

Randomize within the jump.
For the specified draw $V_i$, calculate $U_i=F_i(y_i-)+V_i f_i(y_i)$.

Transform to the normal scale.
Calculate $r_i=\Phi^{-1}(U_i)$, retaining enough precision in $U_i$ for the final calculation.

Bernoulli Responses

For a Bernoulli response, $Y_i=1$ denotes success and $Y_i=0$ denotes failure. Let the fitted success probability be $\hat p_i$, so the fitted failure probability is $1-\hat p_i$. The only two possible outcomes make the intervals especially direct: failure occupies the CDF interval from 0 to $1-\hat p_i$, while success occupies the interval from $1-\hat p_i$ to 1.

Result: Bernoulli intervals. If $Y_i\sim\operatorname{Bernoulli}(\hat p_i)$, then $$ \begin{array}{c|c|c} y_i & F_i(y_i-) & F_i(y_i)\\ 0 & 0 & 1-\hat p_i\\ 1 & 1-\hat p_i & 1 \end{array} $$ Thus $U_i=V_i(1-\hat p_i)$ when $y_i=0$, and $U_i=1-\hat p_i+V_i\hat p_i$ when $y_i=1$.

These formulas follow directly from the two fitted masses: $P(Y_i=0)=1-\hat p_i$ and $P(Y_i=1)=\hat p_i$. For a failure, there is no probability strictly below zero. For a success, the probability strictly below one is the failure probability. The interval endpoints therefore follow from the definition of the CDF.

Worked Example: Bernoulli Success

Suppose the fitted success probability is $\hat p_i=0.30$, the observed response is $y_i=1$, and the random draw is $V_i=0.50$. The probability below the observation is $P(Y_i<1)=P(Y_i=0)=0.70$, and the mass at the observation is $P(Y_i=1)=0.30$.

The randomized probability is $$ U_i=0.70+(0.50)(0.30)=0.85. $$ Therefore $$ r_i=\Phi^{-1}(0.85)\approx 1.0364. $$ The residual is positive because the randomized probability is above one half. The chosen random draw matters: another draw inside the same fitted interval would give another residual for this same observed response.

Grouped-Binomial Responses

A grouped-binomial response records $y_i$ successes among $m_i$ trials. Under the fitted model, suppose $Y_i\sim\operatorname{Binomial}(m_i,\hat p_i)$. For $k=0,\ldots,m_i$, its fitted mass is $$ P(Y_i=k)=\binom{m_i}{k}\hat p_i^k(1-\hat p_i)^{m_i-k}. $$ The lower endpoint for an observed count $y_i$ is the sum of the masses at counts strictly less than $y_i$; the interval width is the mass at $y_i$. In particular, $F_i(0-)=0$ and $F_i(m_i)=1$, so the endpoint counts require no special correction.

Worked Example: Grouped-Binomial Count

Suppose $Y_i\sim\operatorname{Binomial}(5,0.40)$, the observed number of successes is $y_i=2$, and $V_i=0.30$. We calculate the masses at counts zero, one, and two because the lower endpoint is $P(Y_i<2)$.

For zero successes, $$ P(Y_i=0)=\binom{5}{0}(0.40)^0(0.60)^5=0.07776. $$ For one success, $$ P(Y_i=1)=\binom{5}{1}(0.40)(0.60)^4=5(0.40)(0.1296)=0.2592. $$ For two successes, $$ P(Y_i=2)=\binom{5}{2}(0.40)^2(0.60)^3=10(0.16)(0.216)=0.3456. $$

The lower cumulative probability is $F_i(2-)=P(Y_i=0)+P(Y_i=1)=0.07776+0.2592=0.33696$. Adding the specified fraction of the mass at two gives $$ U_i=0.33696+(0.30)(0.3456)=0.33696+0.10368=0.44064. $$ Hence $r_i=\Phi^{-1}(0.44064)\approx -0.1494$. The lower endpoint is not $F_i(2)$: that upper endpoint includes the full mass at two, whereas the randomized point includes only the fraction selected by $V_i$.

Poisson Responses

For a Poisson response with fitted mean $\hat\lambda_i$, the fitted mass at a nonnegative integer $k$ is $e^{-\hat\lambda_i}\hat\lambda_i^k/k!$. For an observed count $y_i$, sum the masses from zero through $y_i-1$ for the lower endpoint, and use the mass at $y_i$ as the interval width. The cumulative probability at the observation is the lower endpoint plus that full mass.

Worked Example: Poisson Count With Exact CDF Endpoints

Suppose $Y_i\sim\operatorname{Poisson}(2)$, the observed count is $y_i=2$, and $V_i=0.50$. The probabilities at zero and one are $e^{-2}$ and $2e^{-2}$, respectively. Therefore the lower endpoint is $$ F_i(2-)=P(Y_i\leq 1)=e^{-2}+2e^{-2}=3e^{-2}\approx 0.406005850. $$ The mass at two is $$ P(Y_i=2)=\frac{2^2e^{-2}}{2!}=2e^{-2}\approx 0.270670566. $$

The upper endpoint, calculated from the exact expression rather than rounded intermediate values, is $$ F_i(2)=F_i(2-)+P(Y_i=2)=5e^{-2}\approx 0.676676416. $$ Randomizing halfway through the jump gives $$ U_i=3e^{-2}+(0.50)(2e^{-2})=4e^{-2}\approx 0.541341133. $$ Thus $r_i=\Phi^{-1}(0.541341133)\approx 0.1039$. The observation is not assigned the upper endpoint automatically; the random draw places it halfway through its fitted probability interval in this example.

Practical Interpretation and Common Pitfalls

The randomized probability $U_i$ is a construction based on an observed response and its fitted discrete distribution. The normal transformation expresses that probability as a standard-normal quantile; it does not change the fitted model or make the response itself continuous. For a fixed observation, changing $V_i$ changes the position within the same CDF jump, so a reproducible calculation should record the randomization or its seed.

Using the upper CDF endpoint as the residual input. The correction is to use $F_i(y_i-)+V_i f_i(y_i)$, not simply $F_i(y_i)$.
Putting the observed mass below the observation. The lower endpoint is $P(Y_i<y_i)$; add the observed mass only as part of the randomized jump.
Rounding before adding. Keep exact expressions or adequate precision for the lower endpoint and mass, then calculate $U_i$ and its normal quantile. Separately rounded terms can produce an incorrectly rounded result.
Assuming the residual is uniquely determined by the response. For a discrete outcome the specified random draw is part of the residual calculation; different draws can produce different values within the same interval.

For an interior count, a useful check is that the randomized probability lies between the lower and upper endpoints. Since $0<V_i<1$ and the observed fitted mass is positive, $$ F_i(y_i-) < U_i < F_i(y_i). $$ A value outside these bounds signals an endpoint, mass, or arithmetic error. In numerical work, choose random draws strictly between zero and one so the inverse normal CDF does not receive zero or one.

Key takeaway: For a discrete response, identify the probability strictly below the observation and the probability at the observation. Randomize within that CDF jump, then apply the standard normal inverse CDF. Bernoulli, grouped-binomial, and Poisson calculations differ in their fitted masses, but use this same interval construction.

Check Your Understanding

For each question, identify the relevant fitted CDF interval before transforming to the normal scale.

A Bernoulli observation is zero with fitted success probability $0.25$. What are the lower and upper interval endpoints, and what expression gives $U_i$ for a draw $V_i$?
For a grouped-binomial observation $y_i=3$, explain why the lower endpoint is $P(Y_i\leq 2)$, not $P(Y_i\leq 3)$.
A Poisson response has fitted mean $1.5$ and observed count zero. Write the lower endpoint, upper endpoint, and randomized probability in terms of $V_i$.
Why can rounding the lower endpoint and observed mass separately change the reported randomized residual?
If two random draws are used for the same discrete observation, what remains fixed and what can change in the resulting residuals?