Robust Empirical Bayes for Hypervariable Genes

Why Robust Variance Moderation Is Needed

In a microarray analysis, each gene has its own fitted linear model and residual variance estimate. When the number of samples is small, a gene's residual variance can be highly unstable. Empirical Bayes moderation stabilizes these estimates by combining gene-specific information with information shared across genes. But an unusually large residual variance can distort the fitted shared prior if the prior is estimated using ordinary, nonrobust summaries.

A gene with unusually high residual variability is often called hypervariable. Such a gene should not automatically be treated as having the same variance distribution as typical genes. Robust empirical Bayes fitting limits the influence of outlying variances when estimating the shared prior, then gives identified outlier genes less prior weight. The gene's observed variance is not discarded: it continues to contribute to that gene's moderated estimate.

Model and assumptions: For gene \(g\), suppose the fitted linear model has residual variance estimate \(s_g^2\) and residual degrees of freedom \(d_g\). Conditional on the true residual variance \(\sigma_g^2\), assume the model errors are independent and normally distributed with common variance within that gene's samples, so that \(d_gs_g^2/\sigma_g^2\sim\chi^2_{d_g}\). Across genes, assume the variances are sufficiently comparable, after accounting for any mean-variance trend, to estimate a shared prior distribution. The design and residual degrees of freedom must be correctly specified; robust fitting does not repair a misspecified gene-wise model.

The Shared Variance Prior and Its Posterior

Use the scaled-inverse-chi-squared parameterization in which a prior for gene \(g\)'s true variance has scale \(s_0^2\) and prior degrees of freedom \(d_{0g}\). In the ordinary fit, prior degrees of freedom are common across genes; in a robust fit, they can be gene-specific. The prior scale describes the typical variance level, while the prior degrees of freedom describe how strongly that prior information is weighted.

Definition: Under the scaled-inverse-chi-squared prior, \(d_{0g}s_0^2/\sigma_g^2\sim\chi^2_{d_{0g}}\). Given the gene's observed residual variance \(s_g^2\), the posterior degrees of freedom are \(d_{0g}+d_g\), and the posterior scale is \[ \tilde{s}_g^2=\frac{d_{0g}s_0^2+d_gs_g^2}{d_{0g}+d_g}. \] This scale is not the posterior mean. When \(d_{0g}+d_g>2\), the posterior mean is \[ E(\sigma_g^2\mid s_g^2)=\frac{d_{0g}s_0^2+d_gs_g^2}{d_{0g}+d_g-2}. \]

The distinction matters: the posterior scale is a weighted average of the prior scale and observed variance, but the mean of an inverse-chi-squared distribution has a denominator two smaller than its degrees of freedom. When \(d_{0g}+d_g\leq 2\), that posterior mean is not finite, even though the posterior scale remains defined.

The update follows directly from multiplying the variance likelihood by the prior. Ignoring factors that do not depend on \(\sigma_g^2\), the likelihood from the normal linear model is proportional to \((\sigma_g^2)^{-d_g/2}\exp[-d_gs_g^2/(2\sigma_g^2)]\). The prior density is proportional to \((\sigma_g^2)^{-d_{0g}/2-1}\exp[-d_{0g}s_0^2/(2\sigma_g^2)]\). Their product is proportional to \[ (\sigma_g^2)^{-(d_{0g}+d_g)/2-1} \exp\left[-\frac{d_{0g}s_0^2+d_gs_g^2}{2\sigma_g^2}\right], \] which is the scaled-inverse-chi-squared posterior with the degrees of freedom and scale given above. This derivation also shows why reducing \(d_{0g}\) reduces the prior's contribution relative to the observed variance.

How Robust Fitting Limits Outlier Influence

The empirical Bayes fit estimates the prior from the collection of gene-wise residual variances, taking account of the residual degrees of freedom. Under a common scaled-inverse-chi-squared prior, the ratio \(s_g^2/s_0^2\) has an \(F\) distribution with numerator degrees of freedom \(d_g\) and denominator degrees of freedom \(d_0\). This gives a distributional basis for estimating the prior parameters across genes.

In a nonrobust fit, exceptionally large variances can strongly affect the estimated spread of this distribution. Robust fitting limits that influence by fitting to a Winsorized distribution: values in the chosen tails are replaced, for the purpose of estimating the prior, by the corresponding tail cutoffs. The precise tail settings are analysis choices and should be reported. The fitted prior then supplies gene-specific prior degrees of freedom: a gene judged to be an outlier receives a smaller \(d_{0g}\), so its unusually high observed variance is pulled less toward the shared prior scale.

Crucially, Winsorization is used to protect the prior fit; it does not replace the gene's actual \(s_g^2\) in the posterior update. A hypervariable gene still contributes its full observed residual variance. Its moderation is weaker because the prior contributes fewer degrees of freedom, not because the observed variance has been erased or capped.

Worked Example: A Hypervariable Gene Receives Less Shrinkage

Suppose a robust prior fit gives a typical prior scale \(s_0^2=1.10\) expression-units squared. Consider a gene with \(s_g^2=9.00\) expression-units squared and \(d_g=4\). Assume the robust fit assigns it \(d_{0g}=2\), while a typical gene would receive \(d_0=8\). These values are hypothetical fitted results, used here to illustrate the update.

For the hypervariable gene, the posterior scale is \[ \tilde{s}_g^2=\frac{(2)(1.10)+(4)(9.00)}{2+4} =\frac{2.20+36.00}{6} =\frac{38.20}{6}=6.367 \] expression-units squared. Its posterior degrees of freedom are \(6\), so its posterior mean exists and is \[ E(\sigma_g^2\mid s_g^2)=\frac{38.20}{6-2} =\frac{38.20}{4}=9.550 \] expression-units squared.

If the same observation were assigned the typical prior degrees of freedom \(d_0=8\), its posterior scale would instead be \[ \frac{(8)(1.10)+(4)(9.00)}{8+4} =\frac{8.80+36.00}{12} =3.733 \] expression-units squared, and its posterior mean would be \(44.80/(12-2)=4.480\) expression-units squared. The robust assignment retains more of the gene's high variance. The posterior scale values, \(6.367\) and \(3.733\), are not posterior means; the respective means are \(9.550\) and \(4.480\).

Worked Calculations and Interpretation

Worked Example: Separating Posterior Scale From Posterior Mean

For another hypothetical gene, let \(s_0^2=0.800\) expression-units squared, \(d_{0g}=3\), \(s_g^2=1.40\) expression-units squared, and \(d_g=5\). The robust prior degrees of freedom are lower than those of a typical gene, but are still positive.

The posterior scale is \[ \tilde{s}_g^2=\frac{(3)(0.800)+(5)(1.40)}{3+5} =\frac{2.40+7.00}{8} =\frac{9.40}{8}=1.175 \] expression-units squared. The posterior degrees of freedom are \(8\). Since \(8>2\), the posterior mean is \[ E(\sigma_g^2\mid s_g^2)=\frac{9.40}{8-2} =\frac{9.40}{6}=1.567 \] expression-units squared. Reporting \(1.175\) as the posterior mean would be incorrect: it is the posterior scale parameter.

Worked Example: Comparing the Amount of Shrinkage

Suppose two hypothetical genes each have \(s_g^2=8.00\) expression-units squared and \(d_g=4\), and the prior scale is \(s_0^2=1.20\) expression-units squared. The robust fit assigns \(d_{0g}=2\) to one high-variance outlier and \(d_{0g}=10\) to the other gene.

The fraction of the posterior scale's numerator contributed by the observed variance is \(d_g/(d_g+d_{0g})\). For the outlier this fraction is \(4/(4+2)=0.667\), whereas for the other gene it is \(4/(4+10)=0.286\). Thus, the outlier's estimate is less dominated by the prior.

The respective posterior scales are \[ \frac{(2)(1.20)+(4)(8.00)}{6} =\frac{34.40}{6}=5.733 \quad\text{and}\quad \frac{(10)(1.20)+(4)(8.00)}{14} =\frac{44.00}{14}=3.143 \] expression-units squared. The first is closer to the observed variance of \(8.00\), as expected with less shrinkage. This comparison concerns the posterior scale; if posterior means are requested, divide the same numerators by \(6-2=4\) and \(14-2=12\), respectively.

Assumptions, Diagnostics, and Alternatives

Robust empirical Bayes is most useful when most genes share a reasonably stable variance distribution but a minority have genuine or unexplained extreme variability. It depends on correctly estimated residual variances and degrees of freedom, a suitable model for the mean-variance relationship, and enough genes to estimate the prior distribution. The earlier material on mean-variance trends is relevant: robust fitting across genes does not substitute for modeling systematic changes in variability with expression level.

Inspect the distribution of residual variances, preferably relative to the fitted mean-variance trend and their gene-specific degrees of freedom. Look for a small extreme upper tail, check whether robust and nonrobust prior estimates differ materially, and examine whether the genes receiving low prior degrees of freedom have plausible model fits. A broad skew or several variance subpopulations may indicate that one shared prior is inadequate rather than merely contaminated by a few outliers. A poor design, dependence not represented in the model, or a systematic variance trend also cannot be fixed just by robustifying the prior fit.

Approach	How the prior is fitted	Implication for an extreme gene
Ordinary empirical Bayes	Uses all genes without robust protection against outlying variances	Extreme values can affect the common prior fit; prior degrees of freedom are typically shared
Robust empirical Bayes	Limits tail influence when estimating the shared prior and permits gene-specific prior degrees of freedom	Hypervariable genes can receive lower prior degrees of freedom and less shrinkage
No variance moderation	Uses each gene's residual variance without cross-gene borrowing	Avoids prior assumptions but retains the instability of small-sample variance estimates

These approaches answer different needs. Robust fitting is not a guarantee that every gene-specific variance is correct, and it does not alter the gene-wise coefficient estimates. Its purpose is to make variance moderation less sensitive to a minority of outliers while preserving useful borrowing for the majority. It also does not itself adjust for testing many genes; that is a separate inferential issue.

Common Mistakes

Calling the posterior scale the posterior mean. The scale divides the weighted numerator by \(d_{0g}+d_g\). The mean divides it by \(d_{0g}+d_g-2\), and exists only when the total degrees of freedom exceed two.
Assuming robust fitting caps a gene's observed variance. Tail handling protects estimation of the prior; the actual \(s_g^2\) remains in that gene's posterior update.
Assuming every gene should receive the same shrinkage. Robust fitting can lower prior degrees of freedom for hypervariable genes, reducing the prior's weight for those genes.
Treating an outlier flag as proof of biological hypervariability. It indicates unusual residual variability under the fitted model. Technical artifacts, model misspecification, or unmodeled groups can also produce it.

Key takeaway: Robust empirical Bayes protects the shared variance-prior fit from extreme genes, then allows identified hypervariable genes to receive smaller prior degrees of freedom. Their observed variances are retained and shrink less. Always distinguish the posterior scale from the posterior mean.

Check Your Understanding

Answer each question, then compare with the solutions below.

A gene has \(d_{0g}=4\), \(s_0^2=0.90\), \(d_g=6\), and \(s_g^2=1.50\). Find its posterior scale and posterior mean, with units of expression-units squared.
For the same \(d_g\), which receives less shrinkage toward the prior scale: a gene with \(d_{0g}=2\) or one with \(d_{0g}=8\)? Explain using the observed-variance weight.
Does robust empirical Bayes replace an extreme gene's observed variance with a Winsorized value in the posterior update?
Name one reason an extreme residual variance might not indicate genuine biological hypervariability.

Solutions. (1) The numerator is \((4)(0.90)+(6)(1.50)=3.60+9.00=12.60\). The posterior scale is \(12.60/(4+6)=1.260\) expression-units squared. The total degrees of freedom are \(10\), so the posterior mean is \(12.60/(10-2)=1.575\) expression-units squared. (2) The \(d_{0g}=2\) gene receives less shrinkage: its observed-variance weight is \(6/(6+2)=0.750\), compared with \(6/(6+8)=0.429\). (3) No. Winsorization limits tail influence in fitting the prior; the gene's actual observed variance remains in its posterior update. (4) A technical artifact, an incorrect design, or an unmodeled source of variation can produce an extreme residual variance.