Estimating the Variance Prior Across Genes

Why Estimate a Variance Prior Across Genes?

In a gene-wise linear model, each gene has its own residual variance, estimated from a limited number of residual degrees of freedom. With few samples, these estimates can be noisy: an unusually small or large sample variance may reflect sampling fluctuation rather than a genuinely unusual gene. Empirical Bayes moderation addresses this by combining each gene's variance estimate with a prior distribution estimated across genes. The earlier material on scaled inverse-chi-squared priors and moderated t-statistics gives the role of that prior; here the focus is how its parameters are estimated.

The prior has two parameters. Its scale, $s_0^2$, describes the central size of gene-wise variances, while its degrees of freedom, $d_0$, controls how widely those variances are spread around the prior. Estimating them requires separating real variation between genes from the sampling variation in each gene's residual variance estimate. Log-variance moments make that separation tractable.

Model and Log-Variance Moments

Definition: For gene $g$, let $s_g^2$ be its residual variance estimate and $d_g$ its residual degrees of freedom. Under the gene-wise normal linear model, conditional on the true residual variance $\sigma_g^2$, $$ \frac{d_g s_g^2}{\sigma_g^2}\sim\chi^2_{d_g}. $$ Across genes, assume the true variances are independent draws from a shared scaled inverse-chi-squared prior: $$ \sigma_g^2\sim\frac{d_0s_0^2}{\chi^2_{d_0}}. $$ Here $s_0^2>0$ is the prior scale and $d_0>0$ is the prior degrees of freedom. This model assumes that, after accounting for the design and residual degrees of freedom, the genes share these prior parameters.

The sampling model says that $s_g^2$ is a scaled chi-squared estimate of $\sigma_g^2$. It is unbiased conditional on $\sigma_g^2$, but its logarithm is not: the expected log of a chi-squared variable differs from the log of its expectation. The correction uses the digamma function $\psi(x)$, the derivative of the log gamma function, and its derivative $\psi_1(x)$, the trigamma function.

If $U\sim\chi^2_d$, then $E[\log U]=\psi(d/2)+\log 2$ and $\operatorname{Var}(\log U)=\psi_1(d/2)$. Applying these facts to the sampling model gives

E[\log s_g^2\mid \sigma_g^2] =\log\sigma_g^2+\psi(d_g/2)-\log(d_g/2), \qquad \operatorname{Var}(\log s_g^2\mid \sigma_g^2) =\psi_1(d_g/2).

Thus define the corrected log variance

z_g=\log s_g^2-\psi(d_g/2)+\log(d_g/2).

Its conditional expectation is $\log\sigma_g^2$. It removes the degrees-of-freedom-dependent bias in the log estimate, but it does not remove the random sampling error. Across genes, that remaining error contributes $\psi_1(d_g/2)$ to the variance of $z_g$.

The prior also has convenient log moments. Since $\sigma_g^2=d_0s_0^2/U$ for $U\sim\chi^2_{d_0}$,

E[\log\sigma_g^2] =\log s_0^2+\log(d_0/2)-\psi(d_0/2), \qquad \operatorname{Var}(\log\sigma_g^2)=\psi_1(d_0/2).

Combining the prior and sampling variation gives $E[z_g]=\log s_0^2+\log(d_0/2)-\psi(d_0/2)$ and $\operatorname{Var}(z_g)=\psi_1(d_0/2)+\psi_1(d_g/2)$. The first term in the variance is between-gene variation in the true variances; the second is sampling noise.

Estimating the Prior Degrees of Freedom

Let $G$ be the number of genes used to estimate the prior, and let $S_z^2$ be the usual sample variance of their corrected log variances $z_g$, using denominator $G-1$. Under a shared prior, the $z_g$ have a common mean even when their residual degrees of freedom differ. The expected sample variance is therefore the prior log-variance contribution plus the average sampling contribution. This yields the method-of-moments estimate

\widehat{v}_0 = S_z^2-\frac{1}{G}\sum_{g=1}^{G}\psi_1(d_g/2).

If $\widehat{v}_0>0$, estimate $d_0$ by solving

\psi_1(\widehat{d}_0/2)=\widehat{v}_0, \qquad\text{or equivalently}\qquad \widehat{d}_0=2\psi_1^{-1}(\widehat{v}_0).

The trigamma function is strictly decreasing on positive arguments, from infinity toward zero, so a positive corrected moment has a unique positive solution. Numerical root-finding is generally used. A larger corrected log-variance dispersion implies a smaller $\widehat{d}_0$, meaning a more diffuse prior. A smaller dispersion implies a larger $\widehat{d}_0$, meaning the gene variances are more concentrated around the prior scale.

This moment calculation can be justified directly. Write $z_g=\log\sigma_g^2+\epsilon_g$, where $E[\epsilon_g\mid\sigma_g^2]=0$, and $\operatorname{Var}(\epsilon_g)=\psi_1(d_g/2)$. The genes are independent and their corrected log variances have the same mean. For independent observations with a common mean but possibly different variances, the expected usual sample variance is the average of their variances. Consequently,

E[S_z^2] =\psi_1(d_0/2)+\frac{1}{G}\sum_{g=1}^{G}\psi_1(d_g/2).

Subtracting the known sampling contribution estimates the prior component. This argument depends on the shared prior mean and independent gene-wise sampling errors; it does not require all $d_g$ to be identical.

Worked Example: Correcting Log-Variance Dispersion

Suppose four hypothetical genes have residual degrees of freedom $d_g=6$. After applying the correction for log bias, their corrected log variances are $-1, 0, 0,$ and $1$. The values are on the log of the model's residual-variance scale.

Their mean is zero, and their sample variance is $$ S_z^2=\frac{(-1-0)^2+(0-0)^2+(0-0)^2+(1-0)^2}{4-1} =\frac{2}{3}=0.666667. $$

For $d_g=6$, the sampling-variance contribution is $\psi_1(3)\approx0.394934$. Therefore, $$ \widehat{v}_0=\frac{2}{3}-0.394934\approx0.271733\approx0.2717. $$

Solving $\psi_1(\widehat{d}_0/2)=0.271733$ numerically gives approximately $\widehat{d}_0=8.3$. This positive estimate indicates detectable between-gene dispersion beyond the sampling noise. The value is an estimated prior degrees of freedom, not a residual degrees of freedom for any one gene.

Estimating the Prior Scale

Once $\widehat{d}_0$ is available, estimate the prior scale from the mean corrected log variance. Let $\bar{z}=G^{-1}\sum_g z_g$. Equating this observed mean to the prior expectation gives

\log\widehat{s}_0^2 =\bar{z}-\log(\widehat{d}_0/2)+\psi(\widehat{d}_0/2), \qquad \widehat{s}_0^2 =\exp\!\left\{\bar{z}-\log(\widehat{d}_0/2)+\psi(\widehat{d}_0/2)\right\}.

The correction involving the digamma function matters: the mean of $\log\sigma_g^2$ is not simply $\log s_0^2$. The resulting estimate $s_0^2$ is a scale parameter of the prior, not generally the arithmetic mean of the true gene variances. The arithmetic prior mean exists only when $d_0>2$, in which case it is $d_0s_0^2/(d_0-2)$.

Worked Example: Recovering the Prior Scale

Suppose a hypothetical set of genes has mean corrected log variance $\bar{z}=0.200$, and the degrees-of-freedom calculation gives $\widehat{d}_0=8.32$. Then $\widehat{d}_0/2=4.16$, $\log(4.16)\approx1.4255$, and $\psi(4.16)\approx1.3005$. Substitution gives $$ \log\widehat{s}_0^2 =0.200-1.4255+1.3005=0.0750, \qquad \widehat{s}_0^2=\exp(0.0750)\approx1.078. $$

On the scale used for these hypothetical residual variances, the estimated prior scale is about $1.078$ variance units. The log-variance scale parameter is estimated from the mean only after accounting for the fitted prior degrees of freedom.

Boundary Cases, Diagnostics, and Alternatives

The corrected moment can be zero or negative. For example, finite-gene sampling may produce a raw variance $S_z^2$ smaller than the estimated average sampling contribution. A negative value cannot equal $\psi_1(d_0/2)$, since trigamma is positive. In this boundary case the data do not support a positive moment estimate of between-gene log-variance dispersion; practical fitting procedures typically use a boundary or limiting estimate rather than interpret the negative number as a valid variance. An extremely large $d_0$ corresponds to a prior concentrated near one variance value.

Worked Example: A Nonpositive Corrected Moment

Suppose three hypothetical genes have residual degrees of freedom $4,6,$ and $10$, and their corrected log variances have sample variance $S_z^2=0.300$. The average sampling contribution is $$ \frac{\psi_1(2)+\psi_1(3)+\psi_1(5)}{3} \approx\frac{0.644934+0.394934+0.221323}{3} \approx0.420397. $$

Thus $\widehat{v}_0=0.300-0.420397=-0.120397$. This is not a permissible estimate of $\psi_1(d_0/2)$. It signals that the observed dispersion is no greater than the estimated sampling contribution, not that the prior variance of log variances is negative. A boundary treatment is needed, and the fitted prior should be interpreted cautiously.

Useful diagnostics begin with the assumptions behind the calculation. Check that residual degrees of freedom and variance estimates are correctly paired gene by gene, and that the residual variance estimates are positive so their logarithms exist. Inspect the distribution of corrected log variances for extreme observations or strong skew. A small number of genes with unusually large variance can inflate $S_z^2$ and reduce the fitted $d_0$. Also assess whether a shared prior is plausible: systematic differences in variance across genes, if left unexplained, can be mistaken for random prior dispersion. A mean-variance trend is one such concern; it requires a model that accounts for expression level rather than treating all genes as exchangeable under one prior.

Moment estimation is not the only option. Pooling all residual sums of squares estimates a common variance and effectively assumes little or no genuine gene-to-gene variance spread; it can be misleading when that spread is real. Using each $s_g^2$ without moderation preserves gene-specific estimates but leaves small-sample instability untouched. Log-moment estimation provides a practical compromise: it estimates a shared prior from many genes and retains gene-specific variance information in moderated inference. More elaborate likelihood-based fitting can use the full distributional model, but it does not remove the need to assess the shared-prior assumptions. Moment estimates are summaries of the gene collection, not guarantees that every gene follows the prior closely.

Common pitfall: Subtracting the sampling contribution from a variance computed with a different denominator can change the estimate. Use the usual sample variance with denominator $G-1$ in the stated moment formula, and carry unrounded values through the subtraction. A nonpositive corrected moment is a boundary warning, not a valid negative prior variance.

Key takeaway: Correct each gene's log residual variance for its degrees-of-freedom-dependent bias, subtract average log-sampling variance from the across-gene sample variance, and invert the trigamma function to estimate $d_0$. Then use the mean corrected log variance to estimate $s_0^2$.

Check Your Understanding

Try each question, then compare your reasoning with the brief solutions.

Why is $\log s_g^2$ corrected by $-\psi(d_g/2)+\log(d_g/2)$?
If the corrected across-gene dispersion increases while sampling degrees of freedom stay fixed, what happens to the estimated prior degrees of freedom, and why?
For $G=5$, suppose $S_z^2=0.80$ and the average sampling contribution is $0.30$. What is the moment estimate of the prior log-variance dispersion?
What does a negative corrected moment mean, and what does it not mean?
Why does estimating $s_0^2$ require $\widehat{d}_0$ as well as $\bar{z}$?

Solutions. 1. The conditional expected log sample variance differs from $\log\sigma_g^2$ by $\psi(d_g/2)-\log(d_g/2)$; subtracting this term removes that bias. 2. The estimated $d_0$ decreases because $\psi_1$ is decreasing, so greater between-gene dispersion corresponds to a more diffuse prior. 3. The estimate is $0.80-0.30=0.50$. 4. It means the observed dispersion does not exceed the estimated sampling contribution; it is not a valid negative prior variance. 5. The expected corrected log variance is $\log s_0^2+\log(d_0/2)-\psi(d_0/2)$, so the scale correction depends on the prior degrees of freedom.