Diagnosing Dispersion and Mean-Variance Fit

Why Diagnose Dispersion?

A pseudobulk model must represent more than the expected count for each gene in each sample. It must also represent how much counts vary around those expectations. If the assumed variance is too small, standard errors can be underestimated and gene-level tests can be poorly calibrated; if it is too large, estimates may be less precise than necessary. Dispersion diagnostics ask whether the model’s variance assumptions are plausible for the observed pseudobulk data, both for individual genes and across genes with different expression levels.

The focus here is the fit of the variance model, not a decision rule for declaring genes differentially expressed. A conspicuous residual pattern is a reason to investigate the model and data; it does not, by itself, establish a biological effect, prove that a sample is erroneous, or justify removing a gene.

The Mean–Variance Model

Let \(Y_{gj}\) be the pseudobulk count for gene \(g\) in sample \(j\), and let \(\mu_{gj}\) be its fitted mean. A common negative-binomial model uses the NB2 variance parameterization:

The mean model must account for the study design and relevant sample-level variables. For example, a log-link model may include a library-size offset \(\log L_j\), a condition effect, and other design terms. A dispersion diagnostic is meaningful only relative to the fitted means and the covariance structure the model specifies. A misspecified mean, omitted batch term, or unmodeled repeated-measures correlation can appear as extra residual variation.

In a mixed model, the conditional variance formula does not by itself describe all marginal variability across samples. Random effects induce covariance and can add marginal variation. Therefore, residuals must be interpreted against the fitted model’s covariance structure: simple Pearson residuals are most direct for independent observations, while mixed-model diagnostics should account for the covariance implied by the random effects. The earlier material on residual diagnostics for pseudobulk models provides the general role of residual checks; here the emphasis is how their variation relates to dispersion and expression level.

Gene-Wise Residual Variation

For independent observations under the NB2 model, a Pearson residual for a gene and sample is the observed-minus-fitted count divided by its model-based standard deviation:

If the fitted variance is adequate, these standardized residuals should not systematically be much more variable than the model predicts. One compact summary is the Pearson residual sum of squares divided by residual degrees of freedom:

Here \(df_g\) is the residual degrees of freedom used for that gene’s fitted model. In a simple model with \(n\) independent observations and \(p\) estimated mean parameters, a rough reference is \(df_g=n-p\). In more complex models, residual degrees of freedom and the appropriate residual calculation depend on the fitting method. Under a well-calibrated model and adequate information, \(Q_g\) is often around 1, but it is not an exact universal cutoff: small samples, estimated parameters, leverage, and dependence affect its reference distribution.

For a mixed model, a covariance-aware version uses residuals transformed so that their fitted covariance is approximately the identity matrix. If \(\mathbf{e}_g\) is the residual vector and \(V_g\) is its fitted covariance matrix, a whitening matrix \(W_g\) satisfies \(W_gV_gW_g^{T}=I\); the transformed residuals are \(W_g\mathbf{e}_g\). Their squared magnitude can then be compared with an appropriate residual degrees-of-freedom reference. Treating correlated observations as independent can make a residual-variation summary misleading.

Assessing the Mean–Dispersion Trend

Across genes, dispersion estimates are often noisy, especially when there are few biological replicates. A mean–dispersion plot places a measure of average expression on one axis and gene-wise dispersion estimates on the other. A fitted trend summarizes how dispersion typically changes with expression level. The trend may be approximately flat or may decrease as mean expression rises, depending on the data and estimation procedure; it should be estimated rather than assumed to have a particular shape.

The diagnostic question is whether individual gene estimates show large, persistent departures from the trend, or whether residual variation is systematically related to mean expression. Such patterns can indicate that a common variance trend is inadequate, that a subset of genes has unusual variability, or that an aspect of the mean or covariance model is not represented. Because individual dispersion estimates are noisy, many workflows stabilize or shrink them toward a trend. A trend or stabilized estimate is a modeling aid, not a claim that every gene has exactly the same dispersion.

A Moment-Based Justification

The form of the NB2 variance explains why dispersion is assessed relative to the mean rather than by raw variance alone. If observations have a common mean \(\mu\) and variance \(\mu+\alpha\mu^2\), solving for \(\alpha\) gives the population moment relation below. Replacing population moments with estimates yields a rough diagnostic estimate; in a regression, fitted means differ across samples and parameter estimation must be accounted for, so this calculation is illustrative rather than a substitute for the model’s estimator.

For example, if hypothetical counts have mean 20 and residual variance 52, the moment relation gives \(\alpha=(52-20)/20^2=32/400=0.0800\). The Poisson component accounts for variance 20, and the additional variance component is 32. This calculation also shows why raw variance cannot be compared fairly across genes with different means: even when \(\alpha\) is identical, the \(\alpha\mu^2\) contribution grows quadratically with the mean.

Practical Workflow, Alternatives, and Failure Modes

A disciplined diagnostic workflow checks the fitted assumptions before treating an unusual dispersion as a gene-specific property.

The main alternatives make different assumptions. A Poisson model sets \(\alpha_g=0\) and is unsuitable when variation materially exceeds the Poisson variance. A negative-binomial model with one common dispersion is simpler but can miss gene-to-gene variability. Gene-wise dispersions are more flexible but can be unstable with limited replication. Trend-based or shrinkage approaches borrow information across genes to stabilize estimates while allowing departures; their quality depends on the trend being a useful summary. Robust standard errors can change uncertainty estimates under some forms of misspecification, but they do not make an incorrect mean or covariance model correct.

Several patterns deserve caution. A large \(Q_g\) for one gene may be driven by a single unusual count; it does not show that all samples for that gene are more variable. Broadly elevated residual variation across many genes may indicate a variance or covariance problem, but can also reflect inadequate mean modeling. A trend that bends sharply at the sparsely populated ends of the expression range may be driven by few genes. Finally, small numbers of biological replicates limit the precision of gene-wise dispersion estimates; many cells within a donor do not replace independent biological replication.

Do not use this diagnostic as a test for excess zeros, as a sample-influence analysis, or as a rule for identifying differentially expressed genes. Those questions require distinct checks. In particular, a dispersion departure should prompt investigation, not automatic exclusion: removing a gene or sample solely because it conflicts with the fitted model can bias inference.