Tutorials › High-Dimensional Genomic Data Analysis › Diagnosing Variance Trends and Hyperparameter Fits

Linear Models for Expression Microarrays · Tutorial 409 of 795

Diagnosing Variance Trends and Hyperparameter Fits

Use residual-variance plots and prior-fit diagnostics to judge whether expression-dependent variance moderation is credible and where it may be failing.

Intermediate 10 min read

What You'll Learn

  • Distinguish observed gene-wise residual variances from the fitted prior variance trend.
  • Diagnose whether variance changes systematically with average expression.
  • Check whether a fitted empirical Bayes prior represents the distribution of gene-wise variances.
  • Recognize hypervariable genes that may influence prior estimates.
  • Identify signs of inadequate moderation and understand the limits of diagnostic plots.
  • Compare a trend prior with a constant prior and ordinary gene-wise inference.

Why Diagnose the Variance Fit?

In expression microarray analyses, a gene-wise linear model estimates a residual variance from the samples available for that gene. When the number of samples is small, this estimate can be noisy. Empirical Bayes moderation stabilizes it by combining gene-specific information with a prior distribution estimated across genes. If residual variability also changes with average expression, the prior may need to vary with expression. The quality of the resulting inference therefore depends on two related questions: does the fitted trend describe the systematic variance pattern, and does the fitted prior describe the remaining cross-gene variability?

Diagnostics address those questions; they do not establish that a model is correct. A smooth-looking curve can conceal a poor distributional fit, and individual genes can depart substantially from a well-fitted prior. The aim is to recognize systematic mismatch, unusual influence on the hyperparameter fit, and circumstances where moderation is too strong or too weak to be credible.

Definition: For gene \(g\), let \(s_g^2\) be the residual variance from its fitted linear model, \(d_g\) its residual degrees of freedom, and \(A_g\) its average expression on the analysis scale. A fitted empirical Bayes prior supplies a prior variance scale \(s_0^2(A_g)\), which may be constant or may depend on \(A_g\), and a prior degrees of freedom \(d_0\). The moderated variance estimate is $$ s_{g,\mathrm{mod}}^2 = \frac{d_0s_0^2(A_g)+d_gs_g^2}{d_0+d_g}. $$ This is a weighted combination of the prior scale and observed residual variance; it is not the arithmetic mean of all genes' observed variances.

Model and Assumptions Behind the Diagnostics

For each gene, suppose the expression measurements satisfy a linear model \(y_g=X\beta_g+\varepsilon_g\), where \(X\) is the design matrix, \(\beta_g\) is the gene-specific coefficient vector, and the residual errors have mean zero and variance \(\sigma_g^2\). The model used for inference assumes an appropriate error distribution (commonly normal errors for these linear-model calculations), a correctly specified design, and independent residual information after accounting for the design. If samples are paired or otherwise dependent, that structure must be represented in the model rather than ignored. Residual degrees of freedom must be positive for the usual gene-wise residual variance estimate to be available.

Across genes, empirical Bayes estimation further assumes that gene-specific variances are sufficiently comparable to estimate a shared prior, possibly conditional on average expression. A standard prior representation is scaled inverse-chi-square: conditional on \(A_g\), the prior variance has scale \(s_0^2(A_g)\) and degrees of freedom \(d_0\). In particular, \(d_0s_0^2(A_g)/\sigma_g^2\) has a chi-square distribution with \(d_0\) degrees of freedom. This is a modeling device, not a claim that all genes have the same true variance. The prior describes a distribution of variances across genes.

A residual-variance plot typically places \(A_g\) on the horizontal axis and \(\log(s_g^2)\) on the vertical axis. The log scale helps display variance patterns across orders of magnitude. The trend to assess is the fitted prior scale \(s_0^2(A_g)\), not a curve that must pass through every observed \(s_g^2\). Individual residual variances are noisy estimates, particularly when \(d_g\) is small, so scatter around the trend is expected.

What to Look for in a Variance-Trend Plot

First check whether the fitted trend follows the broad direction of the cloud. A decline in residual variance as average expression rises, for example, may indicate that low-expression genes have greater variability on the analysis scale. A trend that is nearly flat while the cloud has a sustained slope, or that bends in the opposite direction over a substantial expression range, suggests a mismatch. Also inspect the edges: sparse data at very low or high expression can make a smooth fit unstable, so apparent curvature there may be driven by few genes.

Next consider whether deviations are random scatter or structured departures. A band of points consistently above the curve over an expression region can indicate that the trend is too low there. A fan-shaped spread may indicate that the prior distribution's dispersion does not represent the observed spread, even if its central trend is reasonable. Extreme residual variances should be checked for data or model problems, but should not automatically be deleted: a real hypervariable gene is possible.

Residual degrees of freedom matter when judging the cloud. With few residual degrees of freedom, \(s_g^2\) has high sampling variability even if the gene's variance follows the prior. Thus, a wide cloud alone does not prove that the prior is poorly fitted. If \(d_g\) differs across genes, comparisons should account for that difference. For example, one can inspect log residual variances after allowing for the expected sampling variation at each \(d_g\), rather than treating every vertical deviation as equally surprising.

Worked Example: Detecting a Missed Expression Trend

Suppose a hypothetical diagnostic groups genes into three average-expression ranges. The median residual variances are \(0.100\), \(0.060\), and \(0.040\), while a fitted constant prior scale is \(0.065\) in all three ranges. The units are squared log-expression units.

In the lowest-expression range, the observed median is \(0.100/0.065=1.54\) times the fitted scale. In the middle range, it is \(0.060/0.065=0.923\) times the scale; in the highest range, it is \(0.040/0.065=0.615\) times the scale. The systematic decrease across the ranges is not the pattern expected from random scatter around a constant scale.

This summary suggests that a constant prior may miss expression-dependent variance. It does not, by itself, prove that a particular smooth trend is correct: the individual gene variances are noisy, and the ranges could differ in residual degrees of freedom or data quality. The next diagnostic would examine the full scatter and compare it with a fitted expression-dependent prior scale.

Checking the Prior Fit and Its Influence

A trend can match the central pattern while the distribution around it remains poorly represented. In the scaled inverse-chi-square model, the observed residual variance reflects both the gene's underlying variance and sampling variation. Consequently, the prior scale should not be judged by asking whether \(s_g^2=s_0^2(A_g)\) for each gene. A useful fit check instead asks whether the observed distribution of residual variances, allowing for \(d_g\), is broadly compatible with the fitted prior.

One way to make that comparison is to use simulated or calculated prior-predictive quantiles: the range of residual variances expected under the fitted prior and the gene's residual degrees of freedom. If far more genes than expected lie above the upper quantiles, the fit may understate the high-variance tail; an excess below the lower quantiles may indicate that the prior scale or dispersion is off in the other direction. Such checks should be interpreted across groups of genes, not as a rule that every gene must lie inside a chosen band.

Hyperparameter estimates can also be affected by a small number of extreme genes. A robust empirical Bayes fit is designed to reduce the influence of hypervariable genes on the common prior fit, while still allowing those genes to retain their own variance estimates. Compare fitted trend and prior parameters with and without a suspicious subset, or compare robust and non-robust fits, as a diagnostic of influence. Large changes signal that the shared fit is sensitive to a few genes; they do not prove those genes are erroneous.

Worked Example: An Influential Hypervariable Gene

Consider a hypothetical expression range with \(d_g=4\) for each gene. A fit including all genes gives \(d_0=8\) and \(s_0^2(A_g)=0.040\). One gene has observed residual variance \(s_g^2=0.400\), much larger than the prior scale.

Its moderated variance under that fit is $$ s_{g,\mathrm{mod}}^2 =\frac{8(0.040)+4(0.400)}{8+4} =\frac{0.320+1.600}{12} =0.160. $$ The observed variance is pulled toward the prior scale, but remains above it.

Suppose, for illustration, that refitting with a robust procedure gives \(d_0=2\) and \(s_0^2(A_g)=0.025\), while this gene's \(s_g^2\) and \(d_g\) remain unchanged. Then $$ s_{g,\mathrm{mod}}^2 =\frac{2(0.025)+4(0.400)}{2+4} =\frac{0.050+1.600}{6} =0.275. $$ The larger estimate reflects less prior weight in this illustrative refit. The change from \(0.160\) to \(0.275\) shows that the fitted prior can materially affect moderation for this gene. It is a diagnostic illustration, not a universal direction or size of change: actual robust fitting estimates the prior from all genes under its specified procedure.

A Short Justification for the Moderated Variance

The weighted formula follows from combining the prior variance information with the gene's residual sum of squares. Under the model, the residual sum of squares is \(d_gs_g^2\). The scaled inverse-chi-square prior contributes \(d_0s_0^2(A_g)\) in the same variance-times-degrees-of-freedom form. Multiplying the prior density by the normal-model likelihood adds these contributions and adds their degrees of freedom. Thus, the posterior scaled inverse-chi-square distribution has scale $$ \frac{d_0s_0^2(A_g)+d_gs_g^2}{d_0+d_g}. $$ This is the moderated variance scale used in the moderated linear-model calculation. The derivation explains the diagnostic stakes: if either the trend scale or the prior degrees of freedom is poorly estimated, the variance shrinkage can be poorly calibrated.

Worked Example: Recognizing Weak Moderation

Suppose a hypothetical gene has \(s_g^2=0.090\) squared log-expression units and \(d_g=4\). At its average expression, the fitted prior scale is \(s_0^2(A_g)=0.050\). If the fitted prior degrees of freedom are \(d_0=2\), then $$ s_{g,\mathrm{mod}}^2 =\frac{2(0.050)+4(0.090)}{2+4} =\frac{0.100+0.360}{6} =0.0767. $$ The prior receives \(2/(2+4)=1/3\) of the weight, so the estimate remains relatively close to the observed \(0.090\).

If instead the fitted prior degrees of freedom were \(d_0=20\), with the same prior scale and observed variance, then $$ s_{g,\mathrm{mod}}^2 =\frac{20(0.050)+4(0.090)}{20+4} =\frac{1.000+0.360}{24} =0.0567. $$ Now the prior receives \(20/24=0.833\) of the weight. This much stronger moderation is justified only if the cross-gene prior fit supports that concentration. A very large \(d_0\) paired with visible, systematic departures from the fitted trend is a warning that the prior may be overconfident or misspecified.

Diagnostics, Alternatives, and Limitations

A practical diagnostic sequence is to inspect the expression-variance scatter, assess the fitted scale trend across the expression range, check the distribution of deviations after accounting for residual degrees of freedom, and examine whether a few extreme genes materially alter the prior fit. Then verify that conspicuous patterns are not explained by a design error, unmodeled sample structure, normalization problem, or an inappropriate analysis scale. A variance trend cannot repair a misspecified gene-wise model.

ApproachWhat it assumes or usesDiagnostic implication
Ordinary gene-wise varianceUses each gene's residual variance without cross-gene shrinkageA useful comparison, but noisy estimates can be unstable when residual degrees of freedom are small
Constant empirical Bayes priorUses a common prior scale across expression levelsCan be efficient if no meaningful expression trend remains; a visible systematic slope is a warning
Expression-dependent priorAllows the prior scale to vary with average expressionCan follow a genuine mean-variance pattern, but requires enough genes across the expression range to estimate it

These alternatives do not make diagnostics optional. A constant prior may be preferable when the data support no trend, while a flexible trend can overfit sparse regions. Ordinary inference avoids prior-fit dependence but may be imprecise with little sample information. The appropriate comparison is based on the observed trend, the cross-gene fit, and the design—not on a preference for the largest or smallest moderated variance.

A common mistake is to interpret every point far from the trend as evidence that the trend is wrong, or to delete every high-variance gene until the plot looks neat. The correction is to distinguish a systematic pattern across genes from expected gene-wise sampling noise and from genuine outliers. Diagnose the prior at the distribution level, account for residual degrees of freedom, and investigate extreme genes without assuming in advance that they are errors.

Key takeaway: A credible moderation fit should describe the broad relation between average expression and residual variance and should yield a plausible distribution of gene-wise variances around that trend. The plot is a model check, not proof of validity; its conclusions are limited by the design, data quality, and number of genes supporting each part of the fit.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

  1. Why should a residual-variance plot be interpreted with residual degrees of freedom in mind?
  2. A fitted prior scale is constant, but residual variances decline consistently across expression ranges. What is the main diagnostic concern?
  3. For \(d_0=6\), \(s_0^2(A_g)=0.030\), \(d_g=3\), and \(s_g^2=0.090\), calculate the moderated variance.
  4. If a few hypervariable genes substantially change the fitted prior degrees of freedom, what does that indicate, and what does it not establish?
  5. Why is it incorrect to require every gene's observed residual variance to equal the fitted prior scale?

Solutions. (1) Small residual degrees of freedom produce noisy residual variance estimates, so a broad cloud may arise from sampling variability rather than prior misspecification. (2) The constant prior may miss an expression-dependent variance trend; the scatter and degrees of freedom should be checked before concluding that a trend is needed. (3) \(s_{g,\mathrm{mod}}^2=[6(0.030)+3(0.090)]/(6+3)=(0.180+0.270)/9=0.0500\) squared log-expression units. (4) The shared prior fit is sensitive to those genes. This does not establish that the genes are erroneous; they may be genuinely hypervariable. (5) The prior describes a distribution of true gene variances, while each observed residual variance also contains sampling variation. Individual deviations are therefore expected.

Clinical Biostats