A Two-Group Negative Binomial GLM

From a Count Matrix to a Two-Group Model

In RNA-seq, each entry in the gene-by-sample count matrix is a nonnegative integer: the number of reads assigned to a gene in a sample. A two-condition experiment asks whether the expected count for a gene differs between, for example, control and treated samples. The count is the response, but samples can have different sequencing depths. A model therefore needs both a probability model for counts and a way to account for sample-specific exposure.

A generalized linear model (GLM) connects a response distribution to explanatory variables through a link function. Here the response is modeled with a negative binomial distribution, which allows variability greater than the mean, and the log link keeps expected counts positive. The normalization factor enters as an offset: it adjusts the expected count for sample depth without adding another freely estimated condition effect.

Definition: For gene $g$ in sample $j$, let $Y_{gj}$ be the observed read count, $s_j>0$ the sample’s normalization factor, and $x_j$ the condition indicator, with $x_j=0$ for control and $x_j=1$ for treated. The two-group negative binomial GLM is $$ Y_{gj}\mid x_j,s_j \sim \operatorname{NB}(\mu_{gj},\phi_g),\qquad \log(\mu_{gj})=\log(s_j)+\beta_{0g}+\beta_{1g}x_j. $$ The negative binomial parameterization used here has $\mathbb{E}(Y_{gj})=\mu_{gj}$ and $\operatorname{Var}(Y_{gj})=\mu_{gj}+\phi_g\mu_{gj}^{2}$, where $\phi_g\geq 0$ is the dispersion for gene $g$.

The coefficients and dispersion are gene-specific: the model is fitted separately in its biological meaning for each gene, even though analysis software may process many genes together. The intercept $\beta_{0g}$ describes the log expected count for the control condition at normalization factor $s_j=1$. The condition coefficient $\beta_{1g}$ describes the change in the linear predictor when the indicator changes from control to treated. The parameter $\phi_g$ controls the extra variability beyond a Poisson model.

The offset $\log(s_j)$ is known for the model fit, not estimated as a coefficient alongside the condition term. It represents a sample’s relative sequencing exposure after normalization. Size factors are generally scaled so their typical value is near one; consequently, the intercept refers to the control mean at the reference exposure $s_j=1$, rather than an unadjusted count for a particular library.

Assumptions and Model Specification

The model is useful only when its design and probability assumptions are plausible. First, the response must be a count, and the negative binomial mean–variance relationship must be a reasonable description of its variability. Second, samples are treated as independent biological replicates conditional on the modeled condition and normalization factors. Technical details or shared experimental structure that induce dependence are not represented by this simple two-group model.

The condition indicator must match the experiment: every sample is assigned to one of the two conditions, and both conditions need observations to distinguish a condition effect from the intercept. The design also assumes that the condition is the only modeled systematic difference. If a batch or other measured factor is associated with condition and affects counts, omitting it can confound the estimated condition difference; a covariate-extended model addresses that situation separately.

The offset assumes the supplied normalization factors adequately adjust for differences in effective library size. It does not correct arbitrary sample-specific biases or make the underlying biological populations comparable by itself. Normalization choices and their assumptions matter, particularly when broad composition differences affect many genes.

Worked Example: Assigning the Response, Offset, and Design Term

Suppose a hypothetical experiment has four samples. Their conditions are control, control, treated, treated; their normalization factors are $0.8, 1.0, 1.25, 0.9$. For one gene, the observed counts are $7, 9, 16, 11$. Identify the model inputs and write the mean for each sample.

The response vector for this gene is $Y_g=(7,9,16,11)$. The design indicators are $x=(0,0,1,1)$, and the offsets are $\log(0.8),\log(1.0),\log(1.25),\log(0.9)$. The first two samples are controls; the last two are treated.

Using the general model, the four expected counts are $\mu_{g1}=0.8\exp(\beta_{0g})$, $\mu_{g2}=1.0\exp(\beta_{0g})$, $\mu_{g3}=1.25\exp(\beta_{0g}+\beta_{1g})$, and $\mu_{g4}=0.9\exp(\beta_{0g}+\beta_{1g})$. The observed counts are not set equal to these means: they are random outcomes around the model means.

The design has observations in both conditions, so the condition indicator is not identical to the intercept column. The example states no basis for concluding that the negative binomial variance assumption is adequate; that requires checking the data and experimental context rather than reading it off the design.

How the Offset and Condition Effect Work

Exponentiating the linear predictor makes the roles of the terms explicit. For a control sample, the indicator is zero, so $\mu_{gj}=s_j\exp(\beta_{0g})$. For a treated sample, it is one, so $\mu_{gj}=s_j\exp(\beta_{0g}+\beta_{1g})$. Thus the normalization factor scales the expected count for either condition, while the condition term changes the expected count relative to the control at the same normalization factor.

A short derivation shows why the offset has coefficient one. Suppose the expected count is proportional to the normalization factor, with a condition-specific rate $r_{gc}$ for gene $g$ and condition $c$. Then $\mu_{gj}=s_jr_{gc}$. Taking logarithms gives $\log(\mu_{gj})=\log(s_j)+\log(r_{gc})$. In the two-group model, the log rate is represented by $\beta_{0g}+\beta_{1g}x_j$. Since the exposure contribution is already $\log(s_j)$, its multiplier is fixed at one; estimating that multiplier would no longer encode the intended proportional adjustment for exposure.

Worked Example: Calculating Expected Counts From Parameters

For a hypothetical gene, suppose $\beta_{0g}=\log(10)$, $\beta_{1g}=\log(2)$, and the normalization factors for one control and one treated sample are $0.5$ and $2.0$, respectively. Calculate the model means.

For the control sample, $x_j=0$. Therefore, $$ \mu_{g,\mathrm{control}} =0.5\exp\{\log(10)+\log(2)\times 0\} =0.5\times 10 =5. $$

For the treated sample, $x_j=1$. Therefore, $$ \mu_{g,\mathrm{treated}} =2.0\exp\{\log(10)+\log(2)\} =2.0\times 10\times 2 =40. $$

These are expected counts under the specified model, not predictions that the observed counts must equal 5 and 40. The treated sample’s larger normalization factor contributes to its larger expected count, as does the condition term. Because the factors differ, comparing these two means alone does not isolate the condition term; the model includes the offset to account for that exposure difference.

Design Matrix and Identifiability

The two-group design can be written as a matrix with one row per sample and two columns: an intercept column of ones and the condition indicator. For the sample order control, control, treated, treated, the matrix is

X= \begin{pmatrix} 1&0\\ 1&0\\ 1&1\\ 1&1 \end{pmatrix}.

The first column represents the intercept and the second represents condition. Both columns are needed to distinguish the baseline level from the condition change. The offset is not a column of this design matrix with a freely estimated coefficient; it is added directly to the linear predictor.

Worked Example: Detecting a Design That Cannot Separate Conditions

Consider a hypothetical four-sample experiment in which every sample is treated. The proposed indicators are $x=(1,1,1,1)$. Can this design estimate a separate control intercept and treated-condition effect?

The design matrix would be $$ X= \begin{pmatrix} 1&1\\ 1&1\\ 1&1\\ 1&1 \end{pmatrix}. $$ The two columns are identical. For each sample, the linear predictor contains $\beta_{0g}+\beta_{1g}$, so the data provide no distinction between the two coefficients. There are infinitely many pairs with the same sum; for example, increasing $\beta_{0g}$ by $a$ and decreasing $\beta_{1g}$ by $a$ leaves every fitted linear predictor unchanged.

Therefore, this experiment has no control observations and cannot estimate a treated-versus-control difference using this parameterization. More reads or a different count distribution cannot repair the missing condition comparison; the design needs observations from both groups.

Fitting, Diagnostics, and Alternatives

In practice, fitting the model means supplying the gene’s counts, the normalization factors as offsets, and the two-column design. A negative binomial GLM procedure estimates the gene-specific coefficients under the chosen dispersion specification. This tutorial focuses on what is being modeled, not on the likelihood calculations, dispersion-estimation procedures, or tests used to assess a condition coefficient.

Several checks help identify a poor specification. Verify that each count is matched to the correct sample, condition, and offset; a sample-order mismatch can silently attach the wrong exposure or group label. Confirm that both conditions are represented and that the design columns are not redundant. Compare observed counts with fitted means across samples: persistent condition-specific patterns or very large discrepancies can indicate a missing covariate, an unusual sample, or a mean structure that is too simple. Examine whether variability is substantially greater than expected under a Poisson model; the negative binomial dispersion is intended to accommodate extra variation, but it does not make every outlier or dependence structure harmless.

The main simple alternatives clarify why this model is used. Comparing raw counts ignores differences in library exposure. Comparing normalized counts can be useful descriptively, but does not by itself specify a count distribution or model the mean–variance relationship. A Poisson GLM uses the same log-link and offset structure but imposes $\operatorname{Var}(Y_{gj})=\mu_{gj}$; when biological variation produces overdispersion, that variance can be too small. The negative binomial model relaxes this restriction through $\phi_g$, while still relying on a particular mean–variance form and on suitable normalization.

The two-group model is deliberately limited. It has no batch, pairing, interaction, or additional covariate terms, and it assumes a single condition effect for all samples within each group. If the experimental design contains important structure not captured by those two groups, this model may attribute that structure to condition or leave systematic variation unexplained. A more elaborate design must reflect the actual experiment rather than being added merely to improve apparent fit.

Common Mistakes and Key Takeaways

Treating the offset as an estimated condition effect. The offset is the known $\log(s_j)$ term with coefficient fixed at one. It adjusts expected counts for relative exposure.
Calling an observed count the model mean. A count is one realized outcome; $\mu_{gj}$ is its expected value under the model.
Using the Poisson variance automatically. A Poisson model requires variance equal to the mean. The negative binomial permits extra variation through $\phi_g$, but its assumptions still need consideration.
Trying to estimate a contrast without both groups. If all samples have the same condition indicator, the intercept and condition columns cannot be separated.
Assuming normalization removes all confounding. Size factors address relative sequencing exposure, not every systematic difference between samples or conditions.

Key takeaway: For each gene, the two-group negative binomial GLM models the observed count using a negative binomial response, a log mean, a sample-specific log-offset, and an intercept-plus-condition design. A valid interpretation depends on appropriate normalization, a usable two-group design, and assumptions about variability and sample independence.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions below.

For a sample with $x_j=0$, what is the model mean in terms of $s_j$ and $\beta_{0g}$?
A treated sample has $s_j=1.5$, $\beta_{0g}=\log(4)$, and $\beta_{1g}=\log(1.5)$. Calculate its expected count.
Why is the offset coefficient fixed at one rather than estimated freely?
What design problem occurs if all samples are controls?
State one reason a negative binomial model may be preferred to a Poisson model for RNA-seq counts, and one limitation that remains.

Solutions. (1) With $x_j=0$, $\mu_{gj}=s_j\exp(\beta_{0g})$. (2) Since $x_j=1$, $\mu_{gj}=1.5\exp\{\log(4)+\log(1.5)\}=1.5\times4\times1.5=9$ expected counts. (3) The offset represents the specified proportional exposure adjustment, $\mu\propto s_j$; taking logs gives $\log(s_j)$ with multiplier one. (4) The condition indicator is constant and the condition effect cannot be separated from the intercept; there is no treated-versus-control comparison. (5) The negative binomial allows variance $\mu_{gj}+\phi_g\mu_{gj}^2$, which can represent overdispersion; it still depends on suitable normalization and does not automatically handle confounding or dependence.