Why Derive the Shrunk Estimate?
In a negative binomial model for RNA-seq counts, dispersion describes extra variation beyond the mean-dependent variation in the counts. A gene with limited information may have a noisy gene-wise dispersion estimate, while estimates across genes often follow a systematic dispersion-mean trend. Shrinkage combines these sources of information: it moderates a noisy gene-wise estimate toward the trend, while allowing a more precise estimate to remain closer to its gene-specific value.
The introductory material on shrinkage of dispersion estimates establishes this general idea. Here the calculation is made explicit using a simplified normal model on the log-dispersion scale. This model makes the weighting rule and its dependence on uncertainty transparent; it is a working approximation, not a claim that every software implementation uses exactly this posterior calculation.
The Working Model and Its Assumptions
For a given gene, suppose the log-dispersion estimate has an approximately normal sampling distribution around the true log-dispersion, with sampling variance \(s_g^2\). Suppose also that true log-dispersions vary among genes around the trend, with variance \(\tau^2\):
Here \(s_g^2\) describes uncertainty in the estimate for this gene; \(\tau^2\) describes genuine across-gene spread around the trend on the log scale. They are not interchangeable. The simplified calculation treats \(s_g^2\), \(m_g\), and \(\tau^2\) as known and treats the trend and prior parameters as fixed. In practice they are estimated or approximated from the count data and the collection of genes, so the calculation is best read as an explanation of the weighting principle.
The normal sampling approximation is most plausible when the gene-wise likelihood contains enough information to give a reasonably regular estimate on the log scale. It can be poor when information is extremely limited, the likelihood is highly asymmetric, or an estimate lies at a boundary. The normal prior is a modeling choice for the distribution of true log-dispersions around the trend. Neither assumption says that the raw counts themselves are normally distributed; the counts are modeled with a negative binomial distribution, and this approximation concerns the estimated log-dispersion.
Deriving the Posterior Weights
The posterior density for \(\theta_g\) is proportional to the sampling density multiplied by the prior density. Omitting constants that do not depend on \(\theta_g\), its log density is
Collect the terms involving \(\theta_g\). The coefficient of \(\theta_g^2\) is \(-\frac{1}{2}(1/s_g^2+1/\tau^2)\), and the coefficient of \(\theta_g\) is \(y_g/s_g^2+m_g/\tau^2\). Completing the square therefore gives a normal posterior with precision equal to the sum of the sampling precision and prior precision:
Thus the posterior mean is a precision-weighted average. The weight on the gene-wise estimate is \(\tau^2/(s_g^2+\tau^2)\); the weight on the trend is \(s_g^2/(s_g^2+\tau^2)\). These weights are nonnegative and sum to one. Greater sampling uncertainty \(s_g^2\) gives the trend more weight. Greater across-gene spread \(\tau^2\) makes the prior less concentrated around the trend and gives the gene-wise estimate more weight.
The posterior variance \(v_g^*\) is no larger than either \(s_g^2\) or \(\tau^2\), because adding independent precision reduces uncertainty. The posterior mean is also the posterior mode in this normal model. Exponentiating \(\mu_g^*\) gives the posterior median of \(\alpha_g\), since the posterior for \(\theta_g=\log(\alpha_g)\) is symmetric. It is not generally the posterior mean of \(\alpha_g\): under this model that mean is \(\exp(\mu_g^*+v_g^*/2)\). This distinction matters if a report specifically requires a posterior mean on the dispersion scale.
Worked Example: A Moderately Uncertain Gene
Suppose a hypothetical fitted count model gives a gene-wise dispersion estimate \(\hat\alpha_g=\exp(-1.20)\approx0.3010\), so \(y_g=-1.20\). At this gene’s mean expression, suppose the fitted trend is \(m_g=-1.50\), the gene-wise sampling variance is \(s_g^2=0.1600\), and the across-gene variance is \(\tau^2=0.0400\). Assume the normal sampling approximation and normal prior above are adequate for this illustrative calculation.
The weight on the gene-wise estimate is \(\tau^2/(s_g^2+\tau^2)=0.0400/(0.1600+0.0400)=0.2000\). The weight on the trend is \(0.1600/0.2000=0.8000\). Therefore,
The posterior variance and standard deviation are \(v_g^*=0.1600(0.0400)/0.2000=0.03200\) and \(\sqrt{0.03200}=0.1789\) log-dispersion units. Exponentiating the posterior mean gives \(\exp(-1.440)\approx0.2369\). This is the posterior median dispersion under the model, compared with the raw estimate of about \(0.3010\) and the trend value \(\exp(-1.50)\approx0.2231\). The shrunk estimate lies between the two because it is their weighted average on the log scale.
How Precision Changes the Result
The next two hypothetical calculations hold the gene-wise log estimate and trend fixed. Only the uncertainty changes. This isolates the role of the weights and shows why shrinkage is gene-specific even when genes share the same trend.
Worked Example: A More Precise Gene-Wise Estimate
Suppose another hypothetical gene has \(y_g=-1.20\), trend \(m_g=-1.50\), and a smaller sampling variance \(s_g^2=0.0100\). Keep \(\tau^2=0.0400\). The stated normal approximation and fixed trend assumptions are the same as above; only this gene’s estimated sampling uncertainty is smaller.
The gene-wise weight is \(0.0400/(0.0100+0.0400)=0.8000\), and the trend weight is \(0.2000\). Thus
The posterior variance is \(0.0100(0.0400)/0.0500=0.008000\), with posterior standard deviation \(\sqrt{0.008000}=0.08944\) log-dispersion units. The posterior median dispersion is \(\exp(-1.260)\approx0.2837\). It remains closer to the raw estimate than the result in the first example because the gene-wise estimate is more precise.
Worked Example: A Less Precise Gene-Wise Estimate
Now suppose a hypothetical gene again has \(y_g=-1.20\) and \(m_g=-1.50\), but its sampling variance is \(s_g^2=0.3600\). Let \(\tau^2=0.0400\) as before. Under the same working assumptions, this estimate is much less precise than the prior information around the trend.
The gene-wise weight is \(0.0400/(0.3600+0.0400)=0.1000\), and the trend weight is \(0.9000\). Hence
The posterior variance is \(0.3600(0.0400)/0.4000=0.03600\), and the posterior standard deviation is \(\sqrt{0.03600}=0.1897\) log-dispersion units. Exponentiating gives a posterior median dispersion of \(\exp(-1.470)\approx0.2298\), close to the trend value of about \(0.2231\). The posterior uncertainty is still appreciable: shrinkage moves the point estimate but does not make the gene’s data more informative.
What the Calculation Does—and Does Not—Establish
The examples demonstrate three central properties. First, the estimate is shrunk on the log scale, so it is a weighted average of log-dispersions, not generally a weighted average of dispersions. Second, uncertainty determines the amount of shrinkage: noisy gene-wise estimates move further toward the trend. Third, shrinkage is not the same as replacing every gene’s dispersion by the trend. If the sampling variance is small relative to \(\tau^2\), the gene-wise estimate receives most of the weight.
This calculation relies on treating the trend and variance components as known. If the trend itself is uncertain, treating it as fixed can make the posterior variance appear too small. In applications, trend uncertainty, estimation of \(\tau^2\), and the way the likelihood departs from a normal approximation may all matter. A method may use likelihood-based or empirical-Bayes calculations that differ in details while retaining the broad principle of borrowing information across genes.
A useful diagnostic is to compare the gene-wise estimates with the fitted trend on the log scale and to examine how much uncertainty each estimate has. If almost all estimates are pulled nearly to the trend, check whether the assumed across-gene variance is too small or the gene-wise uncertainties are too large. If estimates remain extremely variable despite substantial stated uncertainty, check the model fit and how the uncertainties were obtained. Genes far from the trend deserve particular attention: an unusual value may reflect real biology, weak information, or a poor approximation. Shrinkage alone cannot decide which explanation is correct.
| Approach | What it does | Main trade-off |
|---|---|---|
| Use each raw gene-wise estimate | Retains each gene’s estimate without borrowing strength. | Avoids trend-induced bias but can leave high-variance estimates unstable. |
| Use the trend for every gene | Assigns each gene its trend value. | Very stable, but discards gene-specific evidence and across-gene variation. |
| Shrink toward the trend | Combines gene-wise evidence and trend information according to uncertainty. | Balances stability and gene-specific information, but depends on the trend and uncertainty model. |
Common Mistakes and Failure Modes
- Reversing the weights. A large \(s_g^2\) means the gene-wise estimate is uncertain, so it increases the weight on the trend, not on the gene-wise estimate.
- Mixing log-scale and dispersion-scale calculations. The displayed weighted average applies to \(y_g\) and \(m_g\). Exponentiate after calculating the log-scale estimate.
- Calling the exponentiated posterior mean the posterior mean dispersion. \(\exp(\mu_g^*)\) is the posterior median dispersion in this model; the mean is \(\exp(\mu_g^*+v_g^*/2)\).
- Assuming shrinkage proves a value is correct. Shrinkage stabilizes estimates under a model. A poor trend, an inaccurate uncertainty estimate, or a non-normal likelihood can still produce misleading results.
- Interpreting a trend-centered result as no gene-specific variation. The posterior mean can be close to the trend even when its posterior variance remains substantial.
Check Your Understanding
For each question, try the calculation before reading its solution.
- If \(s_g^2=0.0300\) and \(\tau^2=0.0700\), what are the weights on the gene-wise estimate and trend?
- For \(y_g=-1.00\), \(m_g=-1.40\), \(s_g^2=0.0900\), and \(\tau^2=0.0100\), calculate the posterior mean log-dispersion.
- In the same setting as question 2, calculate the posterior variance and standard deviation.
- Why is exponentiating the posterior mean log-dispersion not generally the posterior mean dispersion?
Solutions. (1) The gene-wise weight is \(0.0700/(0.0300+0.0700)=0.7000\); the trend weight is \(0.3000\). (2) The gene-wise weight is \(0.0100/0.1000=0.1000\), so the posterior mean is \(0.1000(-1.00)+0.9000(-1.40)=-1.360\). (3) The posterior variance is \(0.0900(0.0100)/0.1000=0.009000\), and the posterior standard deviation is \(\sqrt{0.009000}=0.09487\) log-dispersion units. (4) Exponentiation is nonlinear: \(\exp(\mu_g^*)\) is the posterior median of the dispersion under the log-normal posterior. The posterior mean is \(\exp(\mu_g^*+v_g^*/2)\), which includes the posterior variance.