Residual Diagnostics for Remaining Library-Size Effects

Why Check Residuals After Normalization?

A normalization procedure can produce values that appear comparable while still leaving systematic depth patterns. For example, cells with larger libraries might have positive residuals for many genes, suggesting that the fitted mean model has not accounted for their depth adequately. Conversely, an overall residual average near zero can conceal an increasing trend across library sizes. Diagnostics therefore examine how residuals behave in relation to library size, not only whether they are centered overall.

The introductory material on negative-binomial mean models and Pearson residuals provides the modeling context. Here the goal is narrower: to assess whether residual evidence is consistent with the fitted depth relationship. A diagnostic is evidence about model adequacy, not a proof that all biological or technical structure has been handled correctly.

Model and Residual Quantities

Let \(Y_{gi}\) be the count for gene \(g\) in cell \(i\), and let \(L_i\) denote that cell’s library size. A negative-binomial model represents the conditional count distribution with mean \(\mu_{gi}\) and dispersion \(\phi_g\):

Y_{gi}\mid L_i \sim \mathrm{NB}(\mu_{gi},\phi_g), \qquad \operatorname{Var}(Y_{gi}\mid L_i)=\mu_{gi}+\phi_g\mu_{gi}^{2}.

The fitted mean \(\hat{\mu}_{gi}\) describes the expected count under the estimated depth model, and \(\hat{\phi}_g\) is the estimated dispersion. The Pearson residual is the observed-minus-fitted count divided by the fitted conditional standard deviation.

Definition: The Pearson residual for gene \(g\) in cell \(i\) is

r_{gi} = \frac{Y_{gi}-\hat{\mu}_{gi}} {\sqrt{\hat{\mu}_{gi}+\hat{\phi}_g\hat{\mu}_{gi}^{2}}}.

A positive residual means the observed count exceeds its fitted mean; a negative residual means it falls below that mean. The residual is measured in fitted standard-deviation units, not in counts.

Under a correctly specified negative-binomial model with known parameters, the conditional expectation of \(Y_{gi}-\mu_{gi}\) is zero and its conditional variance is \(\mu_{gi}+\phi_g\mu_{gi}^{2}\). Consequently, dividing by the conditional standard deviation gives a residual with conditional mean zero and variance one. This is a useful reference, not a guarantee that fitted residuals have exactly those properties: means and dispersions are estimated, coefficients may be regularized, and sparse counts can yield discrete, asymmetric residual distributions.

Three Complementary Checks

A residual-versus-library-size plot places \(r_{gi}\) on the vertical axis and a depth measure, commonly \(\log L_i\), on the horizontal axis. A horizontal smooth or binned averages can make a broad trend easier to see than a cloud of points. For an individual gene, the plot asks whether that gene’s fitted residuals vary systematically with depth. A summary across genes, such as the median residual per cell, asks whether many genes show a shared cell-level pattern. These are different summaries and should be labelled accordingly.

Gene-wise correlations provide a numerical summary for each gene. For gene \(g\), calculate the Pearson correlation between its residuals \(r_{gi}\) and \(\log L_i\) across cells. A correlation near zero indicates little linear association in that set of cells; it does not rule out a curved relationship or a depth trend limited to part of the range. Plotting remains important. If outliers dominate, a rank correlation can be a useful supplementary summary, but it answers a different, rank-based question.

Distribution checks describe residuals across cells or genes: their center, spread, asymmetry, and tails. A histogram or quantile plot can expose unusual patterns; summaries stratified by library-size range can show whether spread or centering changes with depth. Do not demand a standard-normal histogram from sparse-count Pearson residuals. The reference mean and variance are more informative than exact normal shape, and even those can be altered by fitting and regularization.

Worked Example: Finding a Trend in Binned Residuals

Suppose a hypothetical gene has the following residuals for eight cells. The first four cells have log library sizes \(0,1,2,3\) and residuals \(-1.2,-0.8,-0.4,0.0\). The other four have log library sizes \(5,6,7,8\) and residuals \(0.1,0.5,0.9,1.3\). Compare the residual averages in the two depth ranges.

For the lower-depth cells, the average residual is

\bar r_{\mathrm{low}} = \frac{-1.2-0.8-0.4+0.0}{4} = \frac{-2.4}{4} = -0.6.

Their mean log library size is \((0+1+2+3)/4=1.5\). For the higher-depth cells,

\bar r_{\mathrm{high}} = \frac{0.1+0.5+0.9+1.3}{4} = \frac{2.8}{4} = 0.7, \qquad \overline{\log L}_{\mathrm{high}}=6.5.

The difference in average residual is \(0.7-(-0.6)=1.3\) residual units over a log-library-size difference of \(6.5-1.5=5\), or \(1.3/5=0.26\) residual units per log unit between these two group averages. The shift from negative residuals at lower depth to positive residuals at higher depth is a concerning pattern for this gene. It is a diagnostic signal, not a formal test or proof that depth alone caused the pattern; the plot and broader gene-level behavior should also be examined.

Worked Example: Calculating a Gene-Wise Correlation

For a second hypothetical gene, take six cells with centered log-library-size values \(x_i=(-2,-1,0,0,1,2)\) and residuals \(r_{gi}=(-0.6,-0.3,0,0.1,0.3,0.5)\). Both vectors have mean zero. Calculate their Pearson correlation.

The sums needed are \(\sum_i x_i^2=10\), \(\sum_i r_{gi}^2=0.80\), and \[ \sum_i x_i r_{gi} =1.2+0.3+0+0+0.3+1.0 =2.8. \] Therefore,

r(x,r_g) = \frac{\sum_i x_i r_{gi}} {\sqrt{\sum_i x_i^2\sum_i r_{gi}^2}} = \frac{2.8}{\sqrt{10(0.80)}} = \frac{2.8}{\sqrt{8}} \approx 0.990.

This is a strong positive linear association in these six hypothetical cells: residuals tend to increase with log library size. The least-squares slope for these centered values is \(2.8/10=0.28\) residual units per log unit. The small number of cells makes this an illustration, not a stable population estimate. The value of the correlation is a prompt to inspect the scatterplot and check whether a small number of observations or a nonlinear shape drives the pattern.

Worked Example: Checking Center and Spread

A third hypothetical gene has eight residuals \((-1.1,-0.8,-0.5,-0.2,0.1,0.3,0.6,1.6)\). Calculate their mean and sample variance, then describe what the distribution check does and does not show.

The residuals sum to zero, so their mean is \(0/8=0\). Their squared values sum to \[ 1.21+0.64+0.25+0.04+0.01+0.09+0.36+2.56=5.16. \] Thus the sample variance is

s_r^2 = \frac{\sum_i(r_{gi}-\bar r_g)^2}{n-1} = \frac{5.16}{7} \approx 0.737, \qquad s_r\approx \sqrt{0.737}=0.858.

The residuals are centered, and their sample spread is below one. The positive value \(1.6\) gives the distribution a longer positive tail than the negative side. These features merit comparison with other genes and with residuals across library-size ranges. They do not, on their own, establish a model failure: fitted residuals need not have sample variance exactly one, and eight values provide limited information about distribution shape.

Assumptions, Interpretation, and Failure Modes

These checks are most interpretable when the fitted mean-variance model is a reasonable description of the counts, the library-size measure is defined consistently, and there is enough spread in library size to reveal a trend. They also depend on the set of cells being compared. If cells differ in composition or other measured characteristics, a residual-depth pattern may reflect a misspecified model, a group imbalance, or another structure—not simply an incorrect normalization factor. A diagnostic plot can reveal a problem without identifying its cause.

Look for more than a nonzero overall average. A curved smooth, a change in residual spread at high depth, or a consistent shift across many genes can matter even if positive and negative residuals cancel overall. Conversely, one unusual gene or cell may dominate a correlation without representing a broad depth effect. Examine both gene-wise patterns and cross-gene summaries, and identify whether a plotted point represents a cell, a gene, or a summary over genes.

A raw-count-versus-library-size plot is an alternative descriptive view, but counts are expected to vary with their fitted means and are not standardized for gene-specific dispersion. A residual plot instead checks deviations from the fitted mean on a standardized scale. A distribution histogram can reveal skew and extreme residuals, while a correlation compresses the relationship to one number; neither replaces the other. Formal significance tests for many gene-wise correlations can also be difficult to interpret because many genes are assessed and residuals share cells and fitted quantities. For model assessment, effect size, shape, consistency, and context are central.

A residual diagnostic has limits. A flat trend does not validate the entire model, and a non-flat trend is not an automatic instruction to change the normalization. Residual patterns may arise from model misspecification, sparse observations, influential cells, or biological differences among cells. In particular, these checks alone cannot decide whether a pattern should be removed as technical or retained as meaningful structure.

Key takeaway: A convincing assessment combines residual-versus-log-library-size plots, gene-wise association summaries, and distribution checks. Seek systematic patterns across depth and genes, but treat each as diagnostic evidence rather than a standalone verdict.

Common Mistake

Wrong reasoning: “The average residual over all cells is zero, so library size has been fully accounted for.” A positive trend at high depth can cancel a negative trend at low depth, producing exactly that overall average. Correction: examine residuals conditional on depth—using plots, depth-stratified summaries, and gene-wise associations—and check whether the pattern recurs across genes. Likewise, a correlation near zero rules out only a strong linear association in the evaluated cells; it does not rule out curvature or a localized trend.

Check Your Understanding

Try each question, then compare your reasoning with the solutions.

A residual-versus-log-library-size plot has negative residuals at low depth and positive residuals at high depth, but the overall mean residual is zero. What does this suggest?
For centered values, \(\sum_i x_i^2=4\), \(\sum_i r_i^2=9\), and \(\sum_i x_i r_i=3\). Calculate the Pearson correlation.
Why is a non-normal histogram not, by itself, proof that a negative-binomial fit is inadequate for sparse counts?
A gene has correlation zero between residuals and log library size. Name one depth-related pattern this result might fail to detect.

Solutions. 1. The residuals have a depth trend despite averaging to zero overall; inspect the conditional pattern and whether it recurs across genes. 2. \(r=3/\sqrt{4(9)}=3/6=0.5\), a moderate positive linear association. 3. Pearson residuals from sparse counts are discrete, and fitted parameters can affect their distribution; exact normality is not the required diagnostic target. 4. A curved or depth-range-specific relationship can have zero overall linear correlation.