Residual Diagnostics for Pseudobulk Models

Why Examine Residuals in Pseudobulk Models?

A gene-level pseudobulk model relates a sample’s count for a gene to its condition and other sample-level predictors, often while accounting for library size and repeated observations from donors. A fitted model can produce an estimated effect even when its assumptions are imperfect. Residual diagnostics ask a different question: does the model reproduce the observed counts without leaving a systematic pattern unexplained?

For a single gene, residuals are summaries of the discrepancy between each observed count and what the fitted model predicts. A residual that is positive or negative is not, by itself, evidence of a problem. The concern is a repeated pattern—for example, residuals mostly positive in one condition, or a trend with fitted values or sample order. Such patterns can indicate that the mean structure, dependence structure, or assumed count distribution does not adequately represent the data. Residuals flag a mismatch; they generally do not identify its cause on their own.

The focus here is on deviance residuals and randomized quantile residuals. Both compare observations with fitted count distributions, but they put those discrepancies on different scales. The introductory material on pseudobulk count models provides the setting: observations are gene counts for biological samples, and a model specifies their distribution conditional on predictors, offsets, and any modeled donor effects.

Model and Residual Definitions

For a fixed-effect count model, let \(Y_i\) be the observed count for sample \(i\), and let \(\widehat{\mu}_i\) be its fitted mean. A negative binomial model is commonly written with variance \(\mu_i+\alpha\mu_i^2\), where \(\alpha\) is a dispersion parameter. A Poisson model is the special case with variance \(\mu_i\). The model also specifies how the mean depends on predictors and any offset, such as a library-size adjustment.

The interpretations below require that the chosen count distribution and mean structure adequately describe the observations, with the specified offset and dependence structure. In an ordinary generalized linear model, observations are assumed independent conditional on their predictors. For repeated pseudobulk samples, independence may instead be assumed conditional on modeled donor effects. Diagnostics are approximate in practice because model parameters are estimated from the same data.

The saturated model gives each observation its own fitted mean, so it can reproduce that observation exactly. The deviance measures the log-likelihood loss when the fitted model replaces those separate means. In a Poisson model, the residual has the explicit form below, using the convention that \(y_i\log(y_i/\widehat{\mu}_i)=0\) when \(y_i=0\).

For a negative binomial model, use the same likelihood-ratio definition with the negative binomial log-likelihood; the Poisson shortcut above does not apply. Deviance residuals are useful for spotting whether observations tend to fall above or below their predictions. Their scale, however, does not make them exactly standard normal, especially for small or discrete counts.

Randomized quantile residuals instead use the entire fitted cumulative distribution. Let \(F_i(y)=P(Y_i\leq y)\) be the fitted cumulative distribution for observation \(i\). Draw \(V_i\) independently from a uniform distribution on \((0,1)\), and define:

Why Randomization Gives a Normal Reference Scale

If the model’s distribution is correct and its parameters are known, the randomized probability \(U_i\) is uniform on \((0,1)\). To see why, the event \(Y_i=y\) has probability \(P(Y_i=y)\). Conditional on that event, \(U_i\) is uniform across the interval from \(F_i(y-1)\) to \(F_i(y)\), whose length is exactly \(P(Y_i=y)\). Thus, across all possible counts, those intervals fill \((0,1)\) in proportion to their lengths, giving a uniform distribution. Applying the standard normal quantile function to a uniform variable gives a standard normal variable.

This is a reference result, not a guarantee that residuals calculated from fitted models will be exactly independent standard normal variables. Estimated parameters make the reference approximate, and repeated samples can remain dependent if the donor structure is not handled correctly. A normal quantile plot can help assess whether the residual distribution is broadly compatible with the reference, but plots against predictors and sample structure are needed to look for systematic patterns.

Looking for Structure in Residuals

Plot residuals against fitted means and against predictors included in the model, such as condition or time. Also consider sample order and donor grouping when they are relevant to the design. A cloud without a clear trend is broadly reassuring, while runs of residuals on one side of zero, curvature, or repeated donor-specific patterns deserve investigation. A normal quantile plot is most directly interpretable for randomized quantile residuals; it does not replace plots that preserve the design structure.

The following example illustrates a condition-related pattern. It deliberately fits a model with one common mean to hypothetical counts from two conditions. The resulting residual pattern is evidence that the fitted model has left a systematic condition difference unexplained; it is not a formal test of a condition effect.

Repeated Measures, Comparisons, and Limitations

For models with donor random effects, the distribution used to calculate a residual matters. A conditional residual evaluates an observation conditional on a specified donor effect, often an estimated effect. A marginal predictive distribution averages over the random-effect distribution. These answer different diagnostic questions: the conditional view examines discrepancies after accounting for donor effects, while the marginal view asks whether predictions for new observations from the population of donors are adequate. Treating conditional residuals as independent standard normal values can be misleading when the donor effects were estimated from the same data.

Deviance and randomized quantile residuals also differ from the Pearson residual, which divides the difference between the observed and fitted counts by a model-based standard deviation. Pearson residuals can be convenient for a quick standardized discrepancy, but they do not directly use the likelihood-ratio scale or transform the fitted discrete distribution to a continuous reference scale. No one residual plot is a complete diagnostic procedure: use residuals that match the fitted model and inspect them in ways that preserve its design.

These diagnostics have limits. A small number of pseudobulk samples makes patterns difficult to judge, and a visually striking residual may arise by chance. Parameter estimation, repeated observations, and randomization affect the reference distribution. Residuals alone cannot validate a model, prove a particular alternative, or establish a biological explanation. In particular, do not infer a specific variance or zero-generating mechanism from a residual pattern alone; those questions require targeted checks of the distributional assumptions.

Common Mistakes and Practical Checks

• Reading a residual as a count difference. A deviance residual is a signed likelihood-based quantity, not \(y_i-\widehat{\mu}_i\). Report the observed and fitted counts separately when the count-scale discrepancy matters.
• Expecting every residual to be standard normal. The standard normal reference for randomized quantile residuals assumes the fitted distribution is correct and known. In a fitted, repeated-measures model, it is approximate and may need to account for dependence.
• Ignoring which distribution generated the residual. The Poisson deviance formula is not the negative binomial formula. Use the likelihood for the family actually fitted.
• Treating one unusual value as proof of model failure. Residual diagnostics seek repeated structure. Check the full set of samples and relevant predictors before drawing conclusions.