Estimating a Smooth Mean-Variance Trend

Why Estimate a Mean–Variance Trend?

Genes with higher average expression often have larger observed variability simply because they are more highly expressed. Ranking genes by variance alone can therefore favor abundant genes, even when their variability is ordinary for that expression level. A smooth mean–variability trend provides a reference: it estimates the typical variability associated with each gene’s average expression, so that unusually variable genes can be identified by their departures from that reference.

The goal is descriptive feature ranking, not a hypothesis test. A large positive departure says that a gene is more variable across the observed cells than the fitted trend predicts at its mean expression. It does not, by itself, identify the source of that variability or show that the gene is differentially expressed.

Gene-Level Summaries and the Trend

Let \(Y_{ig}\) denote the normalized expression of gene \(g\) in cell \(i\), for \(n\) cells. Assume that the same set of cells and the same expression scale are used to summarize the genes being compared. Define the gene’s average expression as \(\bar{Y}_g\), and its sample variance across cells as \(s_g^2\). Following the mean-adjusted dispersion convention, define \(D_g=s_g^2/\bar{Y}_g\) when \(\bar{Y}_g>0\). Both summaries are calculated across cells for each gene.

A common choice is to work on logarithmic scales. Set \(x_g=\log(\bar{Y}_g)\) and \(z_g=\log(D_g)\) for positive means and dispersions, then fit a smooth function \(f\) so that \(z_g\) is approximately \(f(x_g)\) plus a residual. The fitted dispersion on the original scale is \(\widehat{D}_g=\exp\{\widehat{f}(x_g)\}\). A log residual is

Thus, the residual has a relative interpretation: \(\exp(r_g)=D_g/\widehat{D}_g\) is the observed dispersion divided by the fitted dispersion. For example, \(r_g=\log(1.5)\) means the observed dispersion is 1.5 times the fitted value. Log residuals treat proportional differences symmetrically: being twice the trend corresponds to \(\log(2)\), while being half the trend corresponds to \(\log(0.5)=-\log(2)\).

This transformation requires positive values. If some estimated dispersions are zero, a defined alternative is \(z_g=\log(D_g+\delta)\), where \(\delta>0\) is a small, fixed offset chosen on the dispersion scale before fitting. The resulting residual compares \(D_g+\delta\) with the fitted value on that shifted scale; it is not exactly a log ratio of the unshifted dispersions. The offset and its effect on rankings should be reported and checked. Genes with zero mean cannot be placed on the log-mean scale and should not be included in this fit.

Fitting a Smooth Curve

A local polynomial smoother, such as LOESS, estimates the trend from nearby genes rather than imposing one global straight line. For a target location \(x_0\), it assigns each gene a weight that typically decreases with its distance from \(x_0\), then fits a low-degree polynomial to the weighted observations \((x_g,z_g)\). A local linear fit can adapt its slope across the expression range. The neighborhood width, often called the span or bandwidth, controls smoothness: a wider neighborhood gives a steadier but less flexible trend, while a narrow one follows local fluctuations more closely.

Here is the local linear calculation. Write \(u_g=x_g-x_0\), let \(w_g\) be the local weight for gene \(g\), and approximate the trend near \(x_0\) by \(a+b u_g\). The weighted least-squares criterion is \(\sum_g w_g(z_g-a-bu_g)^2\). Define \(S_k=\sum_g w_g u_g^k\) for \(k=0,1,2\), \(T_0=\sum_g w_g z_g\), and \(T_1=\sum_g w_g u_gz_g\). Setting the derivatives of the criterion with respect to \(a\) and \(b\) to zero gives the two normal equations \(S_0a+S_1b=T_0\) and \(S_1a+S_2b=T_1\). Solving them gives the fitted value at the target, which is the intercept \(a\):

This calculation is defined when the local observations provide enough variation in \(u_g\) for \(S_0S_2-S_1^2>0\). At the edges of the observed mean range, there may be fewer useful neighbors on one side, making the estimate less stable. The smoother estimates a typical location on the transformed scale; it does not guarantee that every gene’s residual has constant variance or that residuals follow a normal distribution.

Choosing and Checking the Smoother

The trend is meaningful only if the gene summaries are comparable. Before fitting, the expression values should be on a scale that makes cell-to-cell comparisons appropriate. Normalization should have addressed systematic library-size differences, and filtering should remove genes whose low mean makes their summaries too unstable for the intended comparison. Remaining batch, quality, or composition effects can still create variability that the curve will treat as part of the observed pattern.

The smoother also relies on a practical assumption: across the range being analyzed, the typical transformed dispersion changes smoothly enough for a local approximation to be useful. It does not require a particular parametric mean–variance formula, but it can be unreliable where few genes provide information, especially at the extremes of mean expression. Genes are not independent biological replicates for this purpose; they are points used to estimate a descriptive cross-gene trend. Residual size is therefore not, on its own, a calibrated measure of statistical significance.

Useful diagnostics include plotting log dispersion against log mean with the fitted curve overlaid, and plotting residuals against log mean. Look for systematic curvature in the residual plot, a trend that follows isolated extreme genes, and abrupt changes caused by sparse support at either endpoint. Compare reasonable smoother spans: if small span changes substantially reorder candidates, the ranking is sensitive to the fitting choice. A trend that looks smooth overall can still leave structured residuals, which may indicate an unsuitable expression scale, unremoved technical effects, or groups of genes with different variability patterns.

Weights can be used when there is a defensible estimate of the precision of each gene’s dispersion. In practice, that precision can be difficult to estimate, especially for sparse genes. Unweighted smoothing gives each included gene a direct role; precision weighting can reduce the influence of noisy estimates, but incorrect weights can distort the trend. The weighting rule should be specified and checked rather than chosen solely to obtain a preferred set of candidates.

Alternatives, Limitations, and Common Mistakes

A parametric mean–variance model, such as a negative-binomial model for counts, specifies a particular functional form for how variability changes with the mean. It can be more efficient and interpretable when its assumptions fit the data, but misspecification can leave systematic departures. A smooth trend is more flexible and can capture shapes not specified in advance, at the cost of depending on smoothing choices and adequate gene coverage across the mean range. Ranking genes by raw variance is simpler, but does not adjust for the systematic mean–variability relationship.

The trend residual is a relative, conditional description, not a biological label. Excess dispersion may reflect genuine cell-state heterogeneity, technical artifacts, residual normalization effects, or a mixture of these. Conversely, a gene with a small residual can still have biologically important expression differences among subsets of cells. The trend does not establish causality, identify cell types, or replace downstream modeling of a specific scientific question.

• Mistake: treating a large raw dispersion as evidence of unusually high variability. Correction: compare it with the fitted dispersion at that gene’s mean; a high value may be ordinary for a highly expressed gene.
• Mistake: interpreting a positive residual as a significance test. Correction: it is a departure from a descriptive fitted curve. A calibrated test would require an appropriate sampling model and uncertainty assessment.
• Mistake: choosing a very flexible smoother because it follows every point. Correction: an overly narrow neighborhood can fit noisy gene-level estimates and absorb the very departures of interest. Check residual structure and ranking sensitivity across reasonable spans.
• Mistake: calling a highly variable gene differentially expressed. Correction: variability across cells and a difference in average expression between specified groups are distinct questions.