Residual Diagnostics For Negative Binomial GLMs

Why Diagnose a Negative Binomial Model?

A negative binomial GLM can provide a useful comparison of expression across conditions, but its fitted coefficients and tests are only as credible as the model used to obtain them. Residual diagnostics ask a focused question: does the model describe the observed counts for this gene well enough to support the intended inference? A poorly fitting model may miss a covariate pattern, have an unsuitable dispersion, or be unduly influenced by one sample.

Diagnostics are performed gene by gene. They complement, rather than replace, the design and test choices covered in the earlier material on count models and differential-expression tests. They also do not provide a mechanical rule that labels a gene as valid or invalid. In RNA-seq experiments, the number of samples is often modest, counts are discrete, and parameters are estimated from the same observations being checked. Those features make residuals informative but imperfect evidence.

Model and Residual Definitions

For one gene, let $Y_j$ be the count in sample $j$, and let $\mu_j$ be its fitted mean. A common negative binomial parameterization uses dispersion $\alpha \geq 0$, with conditional mean $\mu_j$ and variance $\mu_j+\alpha\mu_j^2$. A GLM specifies the mean through a design matrix and a log link, commonly with a library-size offset:

Y_j \sim \operatorname{NB}(\mu_j,\alpha), \qquad \operatorname{Var}(Y_j\mid X)=\mu_j+\alpha\mu_j^2,\qquad \log(\mu_j)=\log(s_j)+x_j^\mathsf{T}\beta .

Here $s_j$ is the sample's size factor, $x_j$ is its design row, and $\beta$ is the coefficient vector for this gene. The usual working assumptions include conditionally independent counts across samples, a correctly specified mean structure and offset, and a dispersion that adequately describes the conditional variability. These are assumptions to examine, not facts guaranteed by fitting software.

Definition: The Pearson residual for sample $j$ is the observed-minus-fitted count divided by its model-based standard deviation: $$ r_{Pj}=\frac{Y_j-\widehat{\mu}_j} {\sqrt{\widehat{\mu}_j+\widehat{\alpha}\widehat{\mu}_j^2}}. $$ The deviance residual is the signed square root of that observation's contribution to the likelihood deviance. Its sign is positive when $Y_j>\widehat{\mu}_j$ and negative when $Y_j<\widehat{\mu}_j$.

For the negative binomial parameterization above, define $d_j=r_{Dj}^2$, the deviance contribution. For $\widehat{\alpha}>0$, it is

d_j=2\left[ Y_j\log\left(\frac{Y_j}{\widehat{\mu}_j}\right) -\left(Y_j+\frac{1}{\widehat{\alpha}}\right) \log\left(\frac{1+\widehat{\alpha}Y_j}{1+\widehat{\alpha}\widehat{\mu}_j}\right) \right],\qquad r_{Dj}=\operatorname{sign}(Y_j-\widehat{\mu}_j)\sqrt{d_j}.

The term $Y_j\log(Y_j/\widehat{\mu}_j)$ is defined as zero when $Y_j=0$, by its limiting value. The Poisson limit as $\alpha$ tends to zero gives $d_j=2[Y_j\log(Y_j/\widehat{\mu}_j)-(Y_j-\widehat{\mu}_j)]$. Pearson residuals scale the discrepancy by the fitted standard deviation; deviance residuals express it through the likelihood loss relative to a model that fits that observation exactly.

Why the Deviance Formula Has This Form

The likelihood deviance compares the fitted model's log likelihood with that of the saturated model, which assigns each observation its own fitted mean equal to its observed count. For a single observation, the deviance contribution is twice the saturated log likelihood minus the fitted log likelihood. In the negative binomial log likelihood, terms involving only the observed count and dispersion cancel in this subtraction. The terms that remain give:

d_j=2\left[ Y_j\log\left(\frac{Y_j}{\widehat{\mu}_j}\right) -\left(Y_j+\frac{1}{\widehat{\alpha}}\right) \log\left(\frac{1+\widehat{\alpha}Y_j}{1+\widehat{\alpha}\widehat{\mu}_j}\right) \right].

This is nonnegative because the saturated model cannot fit worse than the constrained fitted model. Taking the square root retains the size of the contribution, while attaching the sign of $Y_j-\widehat{\mu}_j$ shows whether the model overpredicted or underpredicted. This likelihood-based construction explains why a deviance residual need not equal a Pearson residual, even though both summarize the same observation.

Worked Example: Calculating Both Residuals

Suppose a hypothetical gene has observed count $Y_j=12$, fitted mean $\widehat{\mu}_j=8$, and fitted dispersion $\widehat{\alpha}=0.2$. These fitted values are taken as given for this calculation.

The fitted variance is $8+0.2(8^2)=20.8$, so the Pearson residual is

r_{Pj}=\frac{12-8}{\sqrt{20.8}}=\frac{4}{4.5607}=0.8771.

For the deviance contribution, $1/\widehat{\alpha}=5$, and the two logarithmic terms are $12\log(12/8)=4.8656$ and $17\log(3.4/2.6)=4.5605$. Thus

d_j=2(4.8656-4.5605)=0.6102,\qquad r_{Dj}=+\sqrt{0.6102}=0.7811.

Both residuals are positive because the observed count exceeds its fitted mean. Their magnitudes are less than one, so this observation is not a large discrepancy on either scale. They are not identical because one uses a variance standardization and the other a likelihood comparison.

Patterns, Leverage, and Influence

A residual is most useful when viewed in context. Plot residuals against fitted means, important covariates, or sample order. For a reasonably specified mean and variance model, residuals should not show a strong systematic pattern. They need not look exactly normal: counts are discrete, sample sizes may be small, and the fitted values and dispersion are estimated.

Residuals that trend with a predictor can indicate that the mean structure is missing a covariate, interaction, or other systematic feature.
Residual spread that grows or shrinks across fitted means can signal that the variance model or dispersion is inadequate. Check whether the pattern is systematic rather than driven by one or two observations.
Residuals that cluster by sample order or batch can suggest an unmodeled sample-level structure. Consider whether a relevant design variable was omitted; do not infer a cause from the plot alone.
A large residual paired with high leverage deserves an influence check. Leverage reflects how unusual a sample's design values are for estimating the coefficients; it is not itself proof of a problem.

A first-order GLM screening approximation to Cook distance for observation $j$ is

D_j \approx \frac{r_{Pj}^2h_j}{p(1-h_j)},

where $h_j$ is the observation's leverage and $p$ is the number of fitted coefficients. This approximation is a prompt for investigation, not an exact influence calculation for every negative binomial fit. The most direct sensitivity check is to refit the model without the sample and compare the relevant fitted coefficients or test results. A sample should not be removed simply because it makes a residual plot look better.

Worked Example: A Residual Pattern and an Influence Screen

In a hypothetical two-condition experiment with three samples per condition, suppose one gene has counts $8,10,12$ in condition A and $9,11,30$ in condition B. With an intercept for each condition, equal size factors, and a fitted dispersion of $0.1$, the fitted mean in B is its average count: $(9+11+30)/3=16.6667$. This mean structure is consistent with the stated design; the model allows each condition its own mean.

For the sample with count 30, the fitted variance is $16.6667+0.1(16.6667^2)=44.4444$, with standard deviation $6.6667$. Its Pearson residual is

r_P=\frac{30-16.6667}{6.6667}=2.000.

The signed deviance residual is also positive. Using the deviance formula, $Y/\widehat{\mu}=1.8$ and $(1+\alpha Y)/(1+\alpha\widehat{\mu})=4/2.6667=1.5$, so $d=2[30\log(1.8)-40\log(1.5)]=2.830$, giving $r_D=1.682$. This observation has a larger discrepancy than its peers and merits inspection, but a residual alone does not establish that it is erroneous.

Within condition B the three design rows have equal fitted means and equal weights, so the high-count sample has leverage $h=1/3$; there are $p=2$ condition coefficients. The approximate Cook distance is $D\approx 2^2(1/3)/[2(1-1/3)]=1.000$. This is a screening value, not a universal cutoff. Refit without this sample and check whether the condition estimate changes materially before drawing a conclusion.

For a separate arithmetic check of the same approximation, suppose a diagnostic calculation gives $r_P=10/\sqrt{14.4}=2.6352$, $h=0.55$, and $p=2$. Then

D\approx \frac{(10/\sqrt{14.4})^2(0.55)}{2(1-0.55)} =\frac{(100/14.4)(0.55)}{0.90} =4.244.

A relatively large value calls for a sensitivity check; it is not a reason, by itself, to discard the sample.

Worked Example: Recognizing a Mean-Structure Pattern

Consider a hypothetical gene with five samples, equal size factors, and an intercept-only negative binomial model with fitted mean $\widehat{\mu}=5$ and dispersion $\widehat{\alpha}=0.1$. The observed counts in a known processing order are $1,3,5,7,9$. Their mean is 5, so the stated intercept-only fitted mean is consistent with these data.

The fitted variance is $5+0.1(25)=7.5$, with standard deviation $\sqrt{7.5}=2.7386$. The Pearson residuals, in sample order, are

\left( \frac{1-5}{\sqrt{7.5}}, \frac{3-5}{\sqrt{7.5}}, \frac{5-5}{\sqrt{7.5}}, \frac{7-5}{\sqrt{7.5}}, \frac{9-5}{\sqrt{7.5}} \right) =(-1.461,-0.7303,0,0.7303,1.461).

None is especially large in isolation, but the steadily increasing pattern suggests that the order may be associated with the counts. If sample order represents a meaningful batch or covariate, an intercept-only mean may be inadequate. This pattern does not prove a batch effect, and the residuals are not five independent tests; it identifies a feature to investigate using the experimental design.

Global Summaries and Their Limits

The residual deviance is the sum of the deviance contributions, $\sum_j r_{Dj}^2$, and the Pearson statistic is $\sum_j r_{Pj}^2$. They provide overall summaries of discrepancy for one gene. Large values relative to a suitable reference can indicate lack of fit, but a simple chi-square comparison is not automatically reliable for sparse counts, small sample sizes, or estimated dispersion. Estimation also constrains the residuals, so their reference distribution is not generally that of independent, exactly standardized errors.

A useful diagnostic combines these summaries with plots and design knowledge. Pearson residuals are straightforward to interpret in units of model-based standard deviations and are often useful for spotting observations far from their fitted means. Deviance residuals quantify likelihood discrepancy and can be more natural when comparing the fitted model with a saturated one. Neither scale should be treated as normally distributed by default. In particular, a normal quantile plot based on only a few samples has limited ability to diagnose a gene's fit.

A pattern across fitted means may indicate variance misspecification, but deciding whether the dispersion-mean relationship itself is appropriate belongs to a separate inspection of dispersion estimates and trends. Residual diagnostics can reveal that the current fit deserves scrutiny; they do not identify a unique replacement model. Likewise, a single unusual observation is distinct from evidence that the count distribution as a whole is inadequate.

Common Mistake

Wrong reasoning: “One sample has a large residual, so it must be an outlier and should be removed before testing.” This confuses a diagnostic signal with a decision rule. A large residual can arise from a real biological response, a missing covariate, an incorrect mean or variance model, or a measurement problem. Removing the sample without checking these possibilities can bias the comparison.

Correction: Check the sample's design, measurements, and residual pattern; assess leverage; and compare the fitted result with and without the observation. Document any justified exclusion, and report how sensitive the conclusion is to it. A residual is evidence to investigate, not a probability that the model or sample is wrong.

Key takeaway: Use Pearson and deviance residuals to locate mismatches between observed counts and a gene's fitted negative binomial model. Look for patterns as well as large individual values, account for leverage, and treat influence measures as prompts for sensitivity analysis rather than automatic exclusion thresholds.

Check Your Understanding

Work through the questions, then compare your reasoning with the solutions below.

For $Y=7$, $\widehat{\mu}=5$, and $\widehat{\alpha}=0.2$, calculate the Pearson residual.
For $Y=0$, $\widehat{\mu}=5$, and $\widehat{\alpha}=0.2$, calculate the signed deviance residual.
What does a steadily increasing sequence of residuals against sample order suggest, and what does it not prove?
Given $r_P=10/\sqrt{14.4}$, $h=0.55$, and $p=2$, calculate the approximate Cook distance.

Solutions. (1) The fitted variance is $5+0.2(25)=10$, so $r_P=(7-5)/\sqrt{10}=0.6325$. (2) With $Y=0$, $d=2[0-5\log(1/2)]=10\log(2)=6.931$. The observed count is below its fitted mean, so $r_D=-\sqrt{6.931}=-2.633$. (3) A trend suggests that sample order, or a covariate associated with it, may be missing from the mean structure. It does not establish a cause or prove that any individual sample is faulty. (4) $D\approx (100/14.4)(0.55)/[2(0.45)]=4.244$; use this as a reason to check influence and sensitivity, not as an automatic deletion rule.