Why Robust Dispersion Estimation Matters
RNA-seq genes can differ greatly in how much their counts vary beyond what is expected from their means. A small number of genes may have unusually large gene-wise dispersion estimates. These may reflect genuine biology, an omitted covariate, a problematic sample, or uncertainty in estimating dispersion from limited replication. Whatever their source, such estimates can pull an ordinary global dispersion–mean trend upward. That can affect the trend used to stabilize estimates for many genes, not just the unusual ones.
The starting model is a negative binomial generalized linear model for the raw integer counts. For gene \(g\) in sample \(j\), write \(K_{gj}\) for the count, \(s_j\) for a sample-specific size factor, and \(\mathbf{x}_j\) for the design vector. The size factor enters the modeled mean, for example as an offset:
Here \(\alpha_g\) is the gene-specific dispersion. Under this parameterization, the variance is \(\mu_{gj}+\alpha_g\mu_{gj}^{2}\). Size-factor normalization therefore adjusts the expected raw count; it does not make non-integer divided counts into negative binomial observations. As established in the earlier material on dispersion–mean trends and shrinkage, the trend describes how typical dispersion varies with expression level, and gene-wise estimates can be stabilized toward that trend.
What Robust Fitting Changes
Let \(\widehat{\alpha}_g\) be a gene-wise dispersion estimate and let \(\bar{\mu}_g\) denote a chosen expression summary for gene \(g\). The trend is a smooth function \(f\) relating expected dispersion to expression. It is often convenient to fit on the log scale: \(z_g=\log(\widehat{\alpha}_g)\) is related to \(x_g=\log(\bar{\mu}_g)\). This avoids negative fitted dispersions and makes multiplicative differences easier to represent. The particular trend form and the definition of the expression summary depend on the estimation procedure.
Ordinary least squares penalizes a residual by its square, so a gene far above the trend can exert substantial influence. A robust fit replaces the squared-residual penalty with a loss that grows more slowly for large residuals, or uses another resistant summary such as a median within expression ranges. For example, the Huber loss is quadratic for small residuals and linear for large ones:
A robust smooth trend can be estimated by minimizing a sum of such losses between observed log dispersions and the trend, possibly with precision weights and a smoothness penalty. The tuning constant \(c\), weights, and smoothing method are choices of the fitting procedure; there is no single universal robust estimator. The key idea is that one very large residual has less effect than it would under squared-error fitting.
Worked Examples
Worked Example: How Dispersion Changes Count Variability
Suppose two hypothetical genes have the same expected raw count, \(\mu=20\), but their dispersion values are \(0.10\) and \(0.60\). Under the negative binomial variance formula, compare their variances.
For dispersion \(0.10\):
For dispersion \(0.60\):
At the same mean count of 20, the second gene has much greater modeled variability. If its estimate is one of a few such large values, including it without protection in trend fitting can make the estimated typical dispersion too high. The calculation does not show that the second gene is faulty: it only quantifies the model-implied variability at those assumed means and dispersions.
Worked Example: A Median Resists One Extreme Estimate
Consider five hypothetical genes in a narrow expression range, with gene-wise dispersion estimates \(0.08, 0.09, 0.10, 0.11,\) and \(0.60\). For illustration, summarize the local trend by the arithmetic mean or by the median.
The arithmetic mean is
The ordered middle value, and therefore the median, is \(0.10\). In this small illustration, the extreme value raises the mean from \(0.095\), the mean of the first four values, to \(0.196\), while the median remains at \(0.10\). A median-based summary is not automatically a complete smooth-trend estimator, but it demonstrates why resistant summaries can protect the estimated typical dispersion from a minority of large values.
Worked Example: A High-Leverage Gene Can Distort a Trend
Consider a simplified linear fit of log dispersion \(y\) on log mean \(x\). Four hypothetical genes follow \(y=-2-0.2x\) at \(x=0,1,2,3\), giving \(y=-2,-2.2,-2.4,-2.6\). Add one gene at \(x=10\) with \(y=2\). This gene is extreme in both expression and log dispersion.
For all five points, \(\bar{x}=3.2\) and \(\bar{y}=-1.44\). Also, \(\sum x_gy_g=5.2\) and \(\sum x_g^2=114\). Thus
The least-squares slope is \(28.24/62.8=0.4497\), and the intercept is \(-1.44-(0.4497)(3.2)=-2.879\). Without the fifth gene, the four points lie exactly on a line with slope \(-0.2\). The unusual gene has changed the fitted slope from negative to positive.
This example is also a warning: robust loss functions limit the influence of large vertical residuals, but a gene far out in the expression range can have high leverage. A robust procedure based only on downweighting residuals may not fully protect against high-leverage points. Plotting the data, checking the expression range, and assessing whether a trend is supported at its extremes remain important.
Justification: Why a Resistant Summary Helps
The median’s resistance can be seen directly in the five-value example. If the largest value \(0.60\) is replaced by any value \(M\ge 0.11\), the ordered values still have \(0.10\) in the middle. The median therefore stays \(0.10\), however large \(M\) becomes. In contrast, the arithmetic mean becomes
which increases without bound as \(M\) increases. This is the basic reason resistant summaries can reduce the influence of a few extreme estimates. Smooth robust regression extends the same principle across expression levels, although its behavior depends on the fitting method, the distribution of genes, and the locations of extreme points.
How to Diagnose Distortion and Choose an Approach
A dispersion–mean plot is a practical diagnostic: examine gene-wise estimates alongside the fitted trend, preferably on scales that make multiplicative dispersion differences visible. A small number of points far above the trend suggests vertical outliers. Points at the edge of the expression range deserve separate attention because they can have high leverage even if their residuals are not initially large. Also check whether the apparent extremes occur in a narrow region, across the full range, or in a broad group of genes. A broad shift is not the same problem as a few isolated departures.
Robust trend fitting does not replace checking the data and model. Unusual dispersion can arise from genuine between-sample biological variation, unmodeled structure such as a batch effect, or a small number of influential observations. The design and sample information should be considered before interpreting an extreme estimate as an isolated gene property. Robust fitting limits an extreme gene’s impact on the common trend; it does not identify the cause of its variability or guarantee that the negative binomial model is adequate.
Different approaches answer different questions. Ordinary least squares uses all residual magnitudes directly and can be efficient when its working assumptions are suitable, but it is sensitive to large residuals. Robust regression or median-based summaries limit the contribution of isolated extremes, while requiring choices about loss, tuning, smoothing, and possibly weights. Trimming or winsorizing can also reduce the effect of extremes, but discards or alters information and depends on a cutoff. Finally, retaining an unusually large gene-wise estimate rather than pulling it strongly toward the trend can avoid underrepresenting that gene’s variability; it does not by itself improve the global trend.
| Approach | Strength | Main limitation |
|---|---|---|
| Ordinary least-squares trend | Uses all genes in a straightforward fit | A few large residuals or high-leverage genes can dominate |
| Robust trend fit | Limits the influence of a minority of extreme residuals | Depends on choices; does not automatically solve high leverage or broad contamination |
| Trimming or winsorizing | Can cap or remove extreme contributions | Cutoff-dependent and may discard genuine structure |
| Keep an extreme gene-wise estimate | Avoids forcing that gene toward a potentially unsuitable trend | Does not correct distortion of the shared trend or explain the unusual variability |
Common Pitfalls
- Modeling divided normalized counts as negative binomial observations. Dividing raw counts by size factors generally produces non-integer values. In the count GLM, model raw counts and incorporate size factors in the expected means, commonly through offsets.
- Calling every large dispersion estimate an error. A large estimate flags unusual variability relative to a trend; it may be biologically meaningful or caused by features of the experiment that need investigation.
- Assuming robust fitting removes all influence. Robust residual losses are principally designed to limit large residual contributions. High-leverage points and a broad shift affecting many genes can still distort the fit.
- Confusing robustness with shrinkage. Robust fitting protects the trend estimate from a minority of extremes. Shrinkage concerns how gene-wise estimates are combined with a trend; the two steps have related goals but are not interchangeable.
Check Your Understanding
Try each question, then compare your reasoning with the solutions.
- At mean count \(\mu=10\), calculate the negative binomial variance for \(\alpha=0.20\). What does the result describe?
- For dispersion estimates \(0.05, 0.06, 0.07, 0.08,\) and \(0.40\), calculate the mean and median. Which is less affected by the largest value?
- Why can a robust fit still be distorted by a gene at an extreme expression level?
- Does a robust trend fit prove that a gene with unusually large dispersion is erroneous? Explain.
Solutions. (1) \(\operatorname{Var}(K)=10+0.20(10^2)=10+20=30\). This is the model-implied variance of the raw count at the specified mean and dispersion. (2) The mean is \((0.05+0.06+0.07+0.08+0.40)/5=0.66/5=0.132\); the median is \(0.07\). The median is less affected by the single largest value. (3) An extreme expression value can have high leverage in a regression and influence the trend’s shape; residual-based robustness alone may not eliminate that influence. (4) No. Robust fitting limits the gene’s influence on the shared trend. The unusual dispersion may reflect biology, model misspecification, or data issues, and requires contextual investigation.