Pearson Residuals from Fitted NB Models

Why Standardize a Count Around Its Fitted Mean?

A raw count is difficult to interpret without context. An observed count of 4 may be unusually high for a gene-cell pair whose fitted mean is 0.3, but unremarkable for one whose fitted mean is 12. Negative-binomial models account for this mean difference and for the extra variability commonly seen in counts. A Pearson residual expresses how far an observation lies from its fitted mean in units of the model's fitted standard deviation.

In this setting, the fitted mean comes from a gene-wise count model that accounts for sequencing depth, using a size factor or depth covariate as specified by the model. The residual therefore describes an observation relative to its model-predicted count at that cell's depth. It is dimensionless: it is not a count, a fold change, or a probability.

Definition: Let $Y_{gc}$ be the observed count for gene $g$ in cell $c$, let $\hat{\mu}_{gc}$ be its fitted negative-binomial mean, and let $\hat{\phi}_g$ be the fitted dispersion. Under the parameterization $\operatorname{Var}(Y_{gc}\mid\text{model})=\mu_{gc}+\phi_g\mu_{gc}^2$, the Pearson residual is $$ r_{gc}=\frac{Y_{gc}-\hat{\mu}_{gc}}{\sqrt{\hat{\mu}_{gc}+\hat{\phi}_g\hat{\mu}_{gc}^{\,2}}}. $$ It is the observed-minus-fitted count divided by the fitted conditional standard deviation.

Model Assumptions and the Residual Formula

The calculation rests on a fitted mean-variance model. Conditional on the model covariates, the count for gene $g$ in cell $c$ is treated as negative binomial with mean $\mu_{gc}$ and dispersion $\phi_g$. In this convention, $\phi_g\geq 0$, and the variance is $\mu_{gc}+\phi_g\mu_{gc}^2$. The first term represents Poisson sampling variation; the second allows variance to grow more quickly than the mean. When $\phi_g=0$, the variance reduces to the Poisson variance $\mu_{gc}$.

The mean is obtained from the fitted depth model. For example, a log-link model may write $\log(\mu_{gc})$ as a gene-specific intercept plus a depth term and any other included predictors. In an offset formulation, the log size factor enters with coefficient fixed at one; in a regression formulation, depth can enter through an estimated coefficient. The Pearson residual calculation uses the resulting $\hat{\mu}_{gc}$, regardless of which of these allowed mean specifications produced it.

The formula assumes that the fitted mean is positive and that the dispersion used in its denominator is appropriate for that gene. Usually the mean and dispersion are estimated from the data, so this is a fitted-model calculation rather than a comparison to known parameters. A well-specified mean-variance relationship is important: if the model systematically misses the mean or understates variability, the residuals can appear more extreme than their intended scale suggests.

Deriving the Standardization

For a random variable with mean $\mu$ and variance $V$, a Pearson residual standardizes the difference between an observation and its expected value by the square root of its variance. This puts deviations on a common scale: a difference of 2 counts has a different meaning when the standard deviation is 1 than when it is 10. For the negative-binomial model, substitute its conditional mean and variance into that general standardization.

\operatorname{Var}(Y_{gc}\mid\text{model}) =\mu_{gc}+\phi_g\mu_{gc}^2

\frac{Y_{gc}-\mu_{gc}}{\sqrt{\operatorname{Var}(Y_{gc}\mid\text{model})}} = \frac{Y_{gc}-\mu_{gc}}{\sqrt{\mu_{gc}+\phi_g\mu_{gc}^2}}.

In practice, the model parameters are estimated, so replacing $\mu_{gc}$ and $\phi_g$ by their fitted values gives $r_{gc}$ in the definition. Squaring the residual gives the observation's Pearson contribution:

r_{gc}^{\,2} = \frac{(Y_{gc}-\hat{\mu}_{gc})^2} {\hat{\mu}_{gc}+\hat{\phi}_g\hat{\mu}_{gc}^{\,2}}.

This contribution compares a squared discrepancy with its fitted variance. It is nonnegative, while the unsquared residual retains direction: positive means the observed count exceeds its fitted mean, and negative means it falls below. The square root in the denominator is essential because the numerator is measured in counts and the denominator is a standard deviation in counts; their ratio has no units.

Worked Example: A Count Above Its Depth-Adjusted Mean

Suppose a hypothetical gene-cell observation has count $Y_{gc}=8$, fitted mean $\hat{\mu}_{gc}=4$ counts, and fitted dispersion $\hat{\phi}_g=0.25$. Calculate and interpret its Pearson residual.

First calculate the fitted variance:

\widehat{\operatorname{Var}}(Y_{gc}) =4+0.25(4^2)=4+4=8\ \text{counts}^2.

The fitted standard deviation is $\sqrt{8}=2.828$ counts. Therefore,

r_{gc} =\frac{8-4}{\sqrt{4+0.25(4^2)}} =\frac{4}{2.828} =1.414.

The observation is 4 counts above its depth-model fitted mean, or about 1.414 fitted standard deviations above that mean. This is a positive, dimensionless residual. It does not say that the gene is expressed 1.414 times more, nor does it by itself provide a p-value.

How Dispersion Changes the Scale

The dispersion term matters because the same difference from the fitted mean can be less surprising for a gene with larger modeled variability. Consider two hypothetical gene-cell observations, each with count 6 and fitted mean 4. Suppose the first gene has dispersion 0.1 and the second has dispersion 1.0. The numerator is 2 counts in both cases, but the denominators differ.

Worked Example: Comparing Two Dispersions

For the first observation, the fitted variance is $4+0.1(4^2)=5.6$ counts squared, so the residual is

r_1=\frac{6-4}{\sqrt{5.6}}=\frac{2}{2.366}=0.8452.

For the second observation, the fitted variance is $4+1.0(4^2)=20$ counts squared, giving

r_2=\frac{6-4}{\sqrt{20}}=\frac{2}{4.472}=0.4472.

Both observations are 2 counts above their fitted means. The second residual is smaller because its model allows greater variability at that mean. A residual compares the discrepancy with the fitted uncertainty; it is not a direct comparison of raw count differences across genes.

The residual also adjusts for depth through the fitted mean and its corresponding variance. Cells with different depths can have different predicted counts for the same gene. Subtracting each cell's own fitted mean and dividing by its own fitted standard deviation makes the resulting deviations more comparable than raw counts, provided the model captures the depth relationship adequately.

Worked Example: A Zero Count at Two Fitted Means

Take two hypothetical observations for the same gene, with dispersion $\hat{\phi}_g=0.2$. In a lower-depth cell the count is zero and the fitted mean is 0.5; in a higher-depth cell the count is also zero but the fitted mean is 2. Calculate both residuals.

For the lower-depth cell:

r_{\mathrm{low}} =\frac{0-0.5}{\sqrt{0.5+0.2(0.5^2)}} =\frac{-0.5}{\sqrt{0.55}} =\frac{-0.5}{0.7416} =-0.6742.

For the higher-depth cell:

r_{\mathrm{high}} =\frac{0-2}{\sqrt{2+0.2(2^2)}} =\frac{-2}{\sqrt{2.8}} =\frac{-2}{1.673} =-1.195.

Both are below their fitted means, so both residuals are negative. The zero count is a larger standardized shortfall in the higher-depth cell, where the fitted model expected 2 counts rather than 0.5. The comparison is conditional on the fitted model; it does not imply that the higher-depth cell has less biological expression.

Interpreting Residuals and Recognizing Their Limits

Under a correctly specified model with known parameters, the residual has conditional mean zero and variance one. This follows directly from standardization:

E\left[\frac{Y_{gc}-\mu_{gc}}{\sqrt{V_{gc}}}\,\middle|\,\text{model}\right]=0, \qquad \operatorname{Var}\left(\frac{Y_{gc}-\mu_{gc}}{\sqrt{V_{gc}}}\,\middle|\,\text{model}\right)=1,

where $V_{gc}=\mu_{gc}+\phi_g\mu_{gc}^2$ is fixed conditional on the model parameters. This is a statement about the first two moments, not a claim that each residual follows a standard normal distribution. Negative-binomial counts are discrete and may be strongly right-skewed, especially when their means are small. Fitting parameters from the same data also affects the residuals' distribution and can induce dependence among them.

A large absolute residual signals a count that is far from its fitted mean relative to the modeled standard deviation. It is a useful scale for describing discrepancies, but it is not automatically a formal significance test. Any inferential claim would require an appropriate reference distribution and attention to how the model and parameters were estimated.

The denominator must use the fitted negative-binomial variance, not just the fitted mean. Using $\sqrt{\hat{\mu}_{gc}}$ silently substitutes a Poisson variance and can mis-scale residuals when $\hat{\phi}_g$ is appreciable. Conversely, using $\hat{\mu}_{gc}+\hat{\phi}_g\hat{\mu}_{gc}$ is also incorrect: the dispersion term contains the squared mean.

Common pitfall: “The residual is $Y-\hat{\mu}$, so it is a corrected count.” This is wrong because $Y-\hat{\mu}$ is only the numerator. The Pearson residual also divides by the fitted standard deviation, making it dimensionless. For example, a residual of 1.4 does not mean 1.4 counts above expectation; the count difference must be recovered using the fitted standard deviation.

Other failure modes include using a fitted mean from a different cell or gene, mixing a dispersion estimate with a different parameterization, and rounding intermediate quantities too early. Check that the mean and dispersion correspond to the same gene-cell model, that the variance is nonnegative, and that the residual's sign agrees with $Y_{gc}-\hat{\mu}_{gc}$. If $\hat{\mu}_{gc}>0$, the denominator is positive and the sign check is immediate.

What Pearson Residuals Do—and Do Not—Summarize

Pearson residuals retain the direction and standardized size of each gene-cell deviation from a fitted depth model. They can make observations from cells with different fitted means easier to compare on a common model-based scale. They remain conditional on the model: a poor depth specification, an unsuitable dispersion, or unmodeled structure can distort their scale and interpretation.

They are not the only possible residual measure. A raw difference $Y_{gc}-\hat{\mu}_{gc}$ is in counts and preserves the original scale, but its magnitude is difficult to compare when means and variances differ. A deviance residual is based on the likelihood contribution relative to a saturated fit and has a different construction and interpretation. Pearson residuals are specifically variance-standardized differences; replacing them with either alternative changes what the numerical value represents.

Most importantly, a Pearson residual is not a replacement observation from which the original count can be read directly. It records how the count compares with its model-predicted mean and variance. Its sign, magnitude, and units (none) should be reported with that conditional interpretation.

Check Your Understanding

Answer each question, then compare with the solutions below.

A hypothetical observation has $Y=3$, fitted mean 5, and dispersion 0.2. Calculate its Pearson residual and state its sign and interpretation.
If the fitted mean and observed count stay fixed but the dispersion increases, what happens to the absolute Pearson residual? Explain.
Why is a Pearson residual not exactly a standard-normal variable for every gene-cell observation?
A residual is reported as $-1.2$. Does this mean the count is 1.2 below the fitted mean? Explain.

Solutions. For question 1, the fitted variance is $5+0.2(5^2)=10$ counts squared, so $r=(3-5)/\sqrt{10}=-2/3.162=-0.6325$. It is negative because the count is below its fitted mean, by about 0.6325 fitted standard deviations. For question 2, increasing dispersion increases $\hat{\mu}+\hat{\phi}\hat{\mu}^2$, so the denominator grows and the absolute residual decreases, when the numerator is unchanged. For question 3, the negative-binomial count is discrete and may be skewed, particularly at low means; additionally, model parameters are estimated, not known. Mean zero and variance one do not establish a standard-normal distribution. For question 4, no. The residual is dimensionless. The count difference is $-1.2$ multiplied by the fitted standard deviation, so its size in counts depends on the fitted mean and dispersion.