Posterior Variances and Moderated Degrees of Freedom

Why Combine Residual Information Across Genes?

In a gene-wise linear model, the residual sum of squares estimates that gene’s variance, but the estimate can be unstable when the residual degrees of freedom are small. An empirical Bayes prior supplies additional information: it represents a distribution of gene-specific variances centered on a common scale. The update in this tutorial combines the prior information with a gene’s residual information, producing a moderated variance estimate and its posterior degrees of freedom.

The update depends on a conjugate model. Its interpretation is conditional on the prior parameters being specified. In an empirical Bayes analysis those parameters are estimated from the collection of genes and then treated as fixed for this calculation; how to estimate them is a separate problem.

Definition: For gene \(g\), let \(Q_g\) be the residual sum of squares from its fitted linear model, and let \(d_g\) be its residual degrees of freedom. Write the gene-specific error variance as \(\sigma_g^2\). The scaled inverse-chi-squared prior has degrees of freedom \(d_0>0\) and scale \(s_0^2>0\), written \(\sigma_g^2\sim \text{Inv-}\chi^2(d_0,s_0^2)\). Under this parameterization, \(d_0\sigma_g^2/s_0^2\sim\chi^2_{d_0}\). The posterior degrees of freedom are \(d_0+d_g\), and the posterior scale estimate is \(s_{g,\mathrm{post}}^2=(d_0s_0^2+Q_g)/(d_0+d_g)\).

Here \(d_0s_0^2\) acts like prior residual sum of squares, while \(d_0\) acts like prior degrees of freedom. The posterior scale estimate is therefore a weighted average when \(d_g>0\). Since the ordinary residual variance estimate is \(s_g^2=Q_g/d_g\), the same expression can be written

s_{g,\mathrm{post}}^2 = \frac{d_0}{d_0+d_g}s_0^2 + \frac{d_g}{d_0+d_g}s_g^2.

This weighted-average form makes the moderation intuitive: a gene with more residual degrees of freedom receives more weight on its own residual variance; a gene with fewer receives more weight on the prior scale. The residual-sum-of-squares form is more general, including the zero-degree-of-freedom case.

Deriving the Posterior Update

Assume that, conditional on the gene’s variance and the fitted design, the expression errors in the linear model are independent normal errors with mean zero and common variance \(\sigma_g^2\). The model is of full rank for its fitted coefficients. Under these conditions, the residual sum of squares satisfies \(Q_g/\sigma_g^2\sim\chi^2_{d_g}\). The residual component is independent of the estimated coefficients under the normal linear model.

The prior is independent across genes conditional on its common hyperparameters, and follows the scaled inverse-chi-squared distribution defined above. For the update, \(d_0\) and \(s_0^2\) are held fixed. These assumptions give a residual likelihood proportional to \((\sigma_g^2)^{-d_g/2}\exp[-Q_g/(2\sigma_g^2)]\), and a prior density proportional to \((\sigma_g^2)^{-d_0/2-1}\exp[-d_0s_0^2/(2\sigma_g^2)]\).

Multiplying the likelihood and prior gives the posterior kernel:

\begin{aligned} \pi(\sigma_g^2\mid Q_g) &\propto (\sigma_g^2)^{-(d_0+d_g)/2-1} \exp\left\{-\frac{d_0s_0^2+Q_g}{2\sigma_g^2}\right\}. \end{aligned}

This has the scaled inverse-chi-squared form. Matching its power of \(\sigma_g^2\) and its exponential term identifies the posterior degrees of freedom and scale:

d_{g,\mathrm{post}}=d_0+d_g, \qquad s_{g,\mathrm{post}}^2 = \frac{d_0s_0^2+Q_g}{d_0+d_g}.

The numerator is the sum of prior and observed residual sums of squares; the denominator is the sum of their degrees of freedom. This proves the update directly, rather than relying only on the weighted-average interpretation.

Important distinction: \(s_{g,\mathrm{post}}^2\) is the scale parameter of the posterior scaled inverse-chi-squared distribution and the customary moderated variance estimate. It is not generally the posterior mean or posterior mode of \(\sigma_g^2\). If \(d_{g,\mathrm{post}}>2\), the posterior mean is \((d_0s_0^2+Q_g)/(d_{g,\mathrm{post}}-2)\); the posterior mode is \((d_0s_0^2+Q_g)/(d_{g,\mathrm{post}}+2)\). These summaries answer different questions.

Worked Updates

Worked Example: A Gene with Few Residual Degrees of Freedom

Suppose a hypothetical gene has \(Q_g=24\) expression-units-squared and \(d_g=4\). The empirical Bayes prior has \(d_0=6\) and \(s_0^2=2\) expression-units-squared. Assume the gene-wise normal linear model and prior above apply, and treat the prior parameters as fixed for this update.

The ordinary residual variance estimate is \(s_g^2=Q_g/d_g=24/4=6\) expression-units-squared. The posterior degrees of freedom and variance scale are

d_{g,\mathrm{post}}=6+4=10, \qquad s_{g,\mathrm{post}}^2 = \frac{6(2)+24}{6+4} = \frac{36}{10} = 3.6\ \text{expression-units-squared}.

Equivalently, the weighted-average expression gives \((6/10)(2)+(4/10)(6)=1.2+2.4=3.6\). The gene’s estimate is pulled below its noisy ordinary estimate of 6, toward the prior scale of 2. The posterior distribution has 10 degrees of freedom and scale 3.6 expression-units-squared. Its posterior mean is \(36/(10-2)=4.5\) expression-units-squared, not 3.6; that difference illustrates why the reported estimate must be identified as a posterior scale estimate.

Worked Example: More Residual Information Means Less Shrinkage

Consider a second hypothetical gene under the same prior, \(d_0=6\) and \(s_0^2=2\) expression-units-squared. Suppose its residual sum of squares is \(Q_g=80\) expression-units-squared with \(d_g=20\). These values are consistent with an ordinary estimate of \(80/20=4\) expression-units-squared. The normal residual model and fixed-prior assumptions are as in the preceding example.

The update is

d_{g,\mathrm{post}}=6+20=26, \qquad s_{g,\mathrm{post}}^2 = \frac{6(2)+80}{26} = \frac{92}{26} \approx 3.538\ \text{expression-units-squared}.

The ordinary estimate, 4, is already closer to the prior scale 2 than the estimate in the first example, and it is based on more residual degrees of freedom. Accordingly, its posterior scale, about 3.538, remains closer to 4 than the first gene’s posterior scale did to its ordinary estimate of 6. The posterior degrees of freedom are 26, reflecting the 20 residual degrees of freedom plus 6 prior degrees of freedom.

Worked Example: A Valid Model with Zero Residual Degrees of Freedom

Suppose a hypothetical gene is fit by a saturated linear model: the model uses all available residual degrees of freedom, so \(d_g=0\), and its residual sum of squares is \(Q_g=0\). Let the prior have \(d_0=6\) and \(s_0^2=2\) expression-units-squared. This is a valid zero-residual-degree-of-freedom case, not a negative-df model.

The update using the residual sum of squares is well defined:

d_{g,\mathrm{post}}=6+0=6, \qquad s_{g,\mathrm{post}}^2 = \frac{6(2)+0}{6+0} = 2\ \text{expression-units-squared}.

The posterior equals the prior because the residual component supplies no variance information. The ordinary estimate \(s_g^2=Q_g/d_g=0/0\) is undefined, so the weighted-average expression involving \(s_g^2\) cannot be used. The residual-sum-of-squares formula remains valid and makes clear why the prior is unchanged.

Interpreting Degrees of Freedom and Checking the Update

The posterior degrees of freedom \(d_0+d_g\) quantify the information entering the variance update in the scaled inverse-chi-squared model. They are not the number of biological samples, nor are they necessarily an integer: empirical Bayes estimation can produce a non-integer \(d_0\), and the chi-squared distribution is defined for positive real degrees of freedom. When \(d_g\) varies across genes because designs or missing observations differ, the posterior degrees of freedom vary by gene as well.

Several checks help detect mistakes or limitations:

Check the residual degrees of freedom. Negative \(d_g\) is invalid for this residual model and usually signals an incorrect design or degrees-of-freedom calculation. Zero \(d_g\), however, can arise in a valid saturated model. With \(Q_g=0\), the posterior remains the prior and the ordinary estimate is undefined.
Check units and positivity. \(Q_g\) and \(s_0^2\) must both be in squared expression units, and both \(d_0\) and \(s_0^2\) must be positive for the stated proper prior. The resulting scale must also be positive.
Check the weighted-average identity when \(d_g>0\). The result should lie between \(s_g^2\) and \(s_0^2\). If it does not, recheck the numerator, degrees of freedom, and units.
Interpret the hyperparameters conditionally. The derivation treats estimated prior parameters as fixed. If the prior is a poor description of the gene-specific variances, or its estimated scale is unrepresentative, the update can pull estimates in an unhelpful direction.

A common mistake is to say that a zero residual degree of freedom proves the model was specified incorrectly and therefore invalidates the update. That reasoning confuses the ordinary estimate with the residual-sum-of-squares update. In a valid saturated model, \(d_g=0\) and \(Q_g=0\) provide no residual information: substituting these values into the posterior formula returns the prior unchanged. Negative residual degrees of freedom are invalid, but zero is a distinct case.

How This Update Differs from Alternatives

Approach	Variance information used	Consequence
Ordinary gene-wise estimate	\(Q_g/d_g\) for that gene	Uses no cross-gene prior; can be unstable with small \(d_g\), and is undefined when \(d_g=0\).
Common pooled variance	A single pooled residual estimate	Uses one variance value for every gene, rather than retaining gene-specific residual information.
Empirical Bayes update	\(Q_g\), \(d_g\), and prior \(d_0,s_0^2\)	Shrinks each gene’s variance scale toward the prior while retaining gene-specific information; posterior df equal \(d_0+d_g\).

Moderation is not a claim that all genes have exactly the same variance. It is a compromise between a gene’s residual evidence and a prior distribution over gene variances. Its limitations follow from that structure: poor prior hyperparameters can distort estimates, and treating hyperparameters as fixed does not propagate their estimation uncertainty. The update also does not repair a misspecified gene-wise mean model or dependence and distributional problems that violate the assumed residual model.

The posterior variance scale and posterior degrees of freedom are ingredients for later inference. In particular, the moderated \(t\)-statistic uses the variance estimate together with the relevant coefficient standard error and posterior degrees of freedom; the statistic itself is not derived here.

Check Your Understanding

Answer each question, then compare with the solutions below.

For \(d_0=4\), \(s_0^2=3\), \(Q_g=18\), and \(d_g=6\), calculate the posterior scale and posterior degrees of freedom.
For a gene with \(Q_g=12\) and \(d_g=3\), explain why the posterior scale is between the ordinary estimate and the prior scale when \(d_0>0\).
A valid saturated model has \(d_g=0\) and \(Q_g=0\). What is the posterior update, and why is \(Q_g/d_g\) not an alternative calculation?
Does \(s_{g,\mathrm{post}}^2\) generally equal the posterior mean of \(\sigma_g^2\)? State the posterior mean when it exists.
What does a negative residual degree of freedom imply, and how does that differ from zero residual degrees of freedom?

Solutions. (1) \(d_{g,\mathrm{post}}=4+6=10\), and \(s_{g,\mathrm{post}}^2=[4(3)+18]/10=30/10=3\), in the squared units of \(s_0^2\) and \(Q_g\). (2) The ordinary estimate is \(12/3=4\); the posterior scale is \([d_0s_0^2+12]/(d_0+3)\), a weighted average of 3 and 4 with positive weights summing to 1, so it lies between them. (3) The posterior df remain \(d_0\) and the posterior scale remains \(s_0^2\), because the residual contributes neither sum of squares nor degrees of freedom. \(12\) is not the value in this case; generally \(Q_g/d_g=0/0\) is undefined when both are zero. (4) No. If \(d_{g,\mathrm{post}}>2\), the posterior mean is \((d_0s_0^2+Q_g)/(d_{g,\mathrm{post}}-2)\). (5) Negative residual df are invalid for the stated residual model. Zero df can occur in a valid saturated model and, with \(Q_g=0\), leaves the prior unchanged.