Why Replication Matters for Dispersion
A count model needs to describe both the expected expression level of a gene and how much its counts vary around that expectation. In RNA-seq, the negative binomial model allows variability beyond the Poisson model, whose variance equals its mean. The gene-specific dispersion controls the extra variability. Estimating it requires information about variation that remains after the modeled means have been accounted for. Biological replicates provide that information: samples with the same modeled condition can differ, and their residual variation helps reveal whether the Poisson variance is adequate.
Replication does not make an estimate automatically precise. With few samples, the same observations must inform both the condition means and the dispersion. Fitting the means uses degrees of freedom and can absorb some observed variation. A method that treats every difference from a fitted mean as though the mean had been known in advance will overstate the information about dispersion.
A typical mean model is \(\log(\mu_{ij})=\log(s_j)+\mathbf{x}_j^\mathsf{T}\boldsymbol{\beta}_i\), where \(s_j\) is the sample’s library-size offset and \(\mathbf{x}_j\) is its design row. The gene-specific coefficient vector \(\boldsymbol{\beta}_i\) describes the mean pattern allowed by the design; \(\alpha_i\) describes residual count variability around that pattern. This tutorial treats the size factors as given offsets for the gene-wise fit.
The model assumes that, conditional on the design and offsets, sample counts for a gene are independent and follow the specified negative binomial distribution, with a common \(\alpha_i\) across the samples in that fit. It also assumes that the mean model is adequate: for example, condition effects and relevant covariates have been represented in the design. If samples are dependent, the design omits an important source of mean variation, or the negative binomial distribution is unsuitable, the fitted dispersion may absorb those problems rather than represent ordinary biological variability.
How a Gene-Wise Estimate Is Obtained
One approach estimates the mean coefficients and dispersion together by maximizing the gene’s negative binomial likelihood. For a given \(\alpha_i\), the mean coefficients determine the fitted values \(\hat{\mu}_{ij}\). The likelihood can be maximized over the coefficients for each candidate dispersion, and the resulting profile likelihood can then be maximized over \(\alpha_i\geq 0\). This is a likelihood-based estimate, not a direct calculation from a single sample variance. The mean and dispersion are coupled: changing the dispersion changes the likelihood’s assessment of the counts, while changing fitted means changes which variation remains to be explained.
A simpler way to see the role of replication is with a moment calculation. Suppose \(n\) independent replicate counts share the same mean \(\mu\) and dispersion \(\alpha\), and let \(S^2\) be their usual unbiased sample variance. Since \(E(S^2)\) equals the count variance under these conditions,
Replacing the unknown mean and variance by their observed estimates gives the plug-in moment estimate \((S^2-\bar{x})/\bar{x}^2\). A negative value is incompatible with the constraint \(\alpha\geq 0\); it indicates that these data do not provide evidence of extra-Poisson variation by this calculation, not that the true dispersion is negative. The estimate is approximate because it substitutes estimated quantities and does not account for all features of a fitted GLM.
Worked Estimates and What They Tell Us
Worked Example: A Residual Moment Estimate With Two Conditions
Consider hypothetical counts for one gene from two biological replicates in each of two conditions, with equal offsets: condition A has counts 5 and 13, and condition B has counts 13 and 29. The fitted condition means are \(\hat{\mu}_A=(5+13)/2=9\) and \(\hat{\mu}_B=(13+29)/2=21\). The design fits two condition means, one for each group.
The residual sum of squares is
Each condition has two replicates and uses one fitted mean, so each contributes \(2-1=1\) residual degree of freedom. Under the negative binomial variance model, the plug-in expected residual sum of squares is approximately
Equating this approximation to the observed residual sum gives \(160=30+522\alpha\), so \(\hat{\alpha}=130/522\approx0.2490\). This is a plug-in moment estimate, not an exact maximum-likelihood estimate. The arithmetic correctly accounts for estimating each condition mean: using two full variances per condition would incorrectly count variation already consumed in fitting those means.
Worked Example: A Single Group of Replicates
Suppose four hypothetical biological replicates for a gene have counts 2, 6, 14, and 18 under the same condition and equal offsets. Their sample mean is \(\bar{x}=(2+6+14+18)/4=10\). The squared deviations sum to \(64+16+16+64=160\), so the unbiased sample variance is \(S^2=160/(4-1)=53.33\).
The moment estimate is
The observed variance is substantially greater than the Poisson variance of 10, so this simple calculation indicates extra-Poisson variation. The estimate remains uncertain: only four counts contribute, and one unusual count can change the sample variance considerably. It should not be interpreted as a precise biological constant merely because the formula returns four decimal places.
Worked Example: Why One Sample per Condition Is Weak Support
Suppose a two-condition experiment has only one hypothetical count for this gene in each condition: 8 in condition A and 24 in condition B. With separate condition means, the design has two observations and fits two means. Thus the residual degrees of freedom are \(2-2=0\). The fitted means are 8 and 24, and both residuals are zero.
A residual moment calculation cannot estimate within-condition variation here: there are no within-condition replicate differences, and the residual degrees of freedom are zero. This does not prove that the gene’s dispersion is zero. It means these observations provide no replicate-based residual estimate for it. A likelihood optimizer might still return a numerical value under a particular parameterization or constraints, but that value should not be mistaken for dispersion information supplied by replication; with one observation per fitted mean, gene-wise estimation is generally weakly supported.
Likelihood Estimates, Assumptions, and Diagnostics
For a likelihood fit, the negative binomial probability mass function under this course’s mean-dispersion parameterization is
The gene likelihood is the product of these probabilities across samples, or equivalently the sum of their log probabilities. For each candidate \(\alpha_i\), the mean coefficients are fitted under the design and offsets. Maximizing the resulting profile likelihood over nonnegative dispersion values gives a gene-wise maximum-likelihood estimate. At \(\alpha_i=0\), the model is understood as the Poisson limit. Because the mean coefficients are estimated from the same data, the profile likelihood is the relevant likelihood for comparing dispersion values; a residual variance calculation is a useful approximation, not a replacement for that joint fitting.
The fit relies on independent biological samples and a suitable mean design. Replicates should represent the independent experimental units relevant to the scientific question. Repeated technical measurements of the same biological material do not, by themselves, provide the same information as independent biological replication. The model also assumes a single dispersion for the gene across samples in the fitted experiment. If there are very few replicates, even a correctly specified likelihood may be broad or sensitive to individual observations.
Warning signs include a dispersion estimate at the zero boundary, a very large estimate driven by one count, poor optimizer convergence, or a profile likelihood that changes little across a wide range of dispersions. These signs do not all have the same interpretation: a boundary estimate can arise when the observed variation is consistent with Poisson variation, whereas a flat profile indicates that the data do not sharply distinguish dispersion values. A poor mean fit can also distort dispersion estimation, so the condition pattern and covariates in the design must be scientifically defensible.
Comparing Approaches and Avoiding Common Errors
| Approach | What it uses | Main caution |
|---|---|---|
| Gene-wise likelihood | The negative binomial likelihood across samples, with mean coefficients fitted under the design | Can be imprecise or unstable when a gene has few informative replicates |
| Moment estimate | Observed variability compared with the negative binomial variance formula | Plug-in estimates and simple residual calculations are approximations; account for fitted means and residual degrees of freedom |
| Poisson model | The special variance relationship \(\operatorname{Var}(Y)=\mu\) | It cannot represent extra-Poisson variability when that variability is present |
The central pitfall is to equate the residual sum of squares to one variance contribution for every observation even though the means were estimated from those observations. In the two-condition example, that mistake would use \(2(9+81\alpha)+2(21+441\alpha)=60+1044\alpha\). It incorrectly doubles the expected within-condition residual variation. The correct simple group-mean calculation uses one residual degree of freedom per condition, giving \(30+522\alpha\). In a general design, the exact relationship is more involved, but the principle remains: fitting mean parameters affects the residual information available for estimating dispersion.
A second error is to treat a gene-wise estimate as equally reliable for every gene. Each gene is fitted separately here, so a gene with few samples or little count information may have a highly variable estimate. Approaches that use information across genes are another strategy, but they answer a different estimation problem; a stand-alone gene-wise estimate should not be described as though it had borrowed that information. Nor should a large dispersion automatically be attributed to biological heterogeneity: unmodeled design effects, dependence, or an outlying observation can also inflate it.
Check Your Understanding
Try each question, then compare your reasoning with the solutions below.
- For a negative binomial count with mean 12 and dispersion 0.2, what variance does the model imply?
- Three independent counts share a mean of 10 and have dispersion \(\alpha\). How many residual degrees of freedom remain after fitting their common mean, and what is the plug-in expected residual sum of squares?
- Four replicate counts have mean 8 and unbiased sample variance 20. Calculate the simple moment estimate of dispersion.
- Why is a gene-wise residual moment estimate unsupported by within-condition replication when there is one sample per condition and a separate mean is fitted for each condition?
- A fitted dispersion is close to zero. Does this establish that the gene has no biological variability? Explain.
Solutions. (1) \(\operatorname{Var}(Y)=\mu+\alpha\mu^2=12+0.2(144)=40.8\) counts squared. (2) The fitted common mean uses one parameter, leaving \(3-1=2\) residual degrees of freedom. The expected residual sum of squares is approximately \(2(10+100\alpha)=20+200\alpha\). (3) \(\hat{\alpha}_{\mathrm{mom}}=(20-8)/8^2=12/64=0.1875\). (4) The design has as many fitted condition means as observations, leaving zero residual degrees of freedom; there are no replicate differences within conditions from which to estimate residual variation. (5) No. It says the fitted data do not provide evidence of extra-Poisson variation under this model; limited replication or information can make the estimate uncertain, and a boundary estimate is not proof of zero biological variability.