The Negative Binomial GLM Likelihood

Why the Likelihood Matters

A negative binomial generalized linear model (GLM) connects observed RNA-seq counts to experimental conditions through the expected count. The likelihood makes that connection precise: it assigns a probability to the observed counts under a proposed set of regression coefficients. Estimating the coefficients means finding values that maximize this probability, or equivalently its logarithm.

Here the focus is the gene-specific likelihood with a log link and a library-size offset. The offset adjusts the expected count for differences in sequencing depth; it is not a coefficient estimated from the gene’s counts. The dispersion is also treated as fixed while deriving and maximizing this likelihood. Estimating gene-wise dispersion is a separate task.

Definition: For gene \(g\) and sample \(j\), let \(Y_{gj}\) be the observed nonnegative integer count, and let \(X\) be the design matrix with row \(x_j^\mathsf{T}\) for sample \(j\). Conditional on the design and known positive size factor \(s_j\), assume the sample counts for this gene are independent and \(Y_{gj}\sim\operatorname{NB}(\mu_{gj},\alpha_g)\), with mean \(\mu_{gj}>0\), dispersion \(\alpha_g>0\), and variance \(\mu_{gj}+\alpha_g\mu_{gj}^2\). The log-linked mean is \(\log(\mu_{gj})=x_j^\mathsf{T}\beta_g+\log(s_j)\). Here \(\beta_g\) is the gene’s coefficient vector and \(\log(s_j)\) is a known offset.

From the Count Model to the Log-Likelihood

Use the mean-dispersion parameterization in which the negative binomial size parameter is \(r_g=1/\alpha_g\). For one count \(y\), with mean \(\mu\) and dispersion \(\alpha\), the probability mass function is

P(Y=y)= \frac{\Gamma(y+1/\alpha)}{\Gamma(1/\alpha)\,\Gamma(y+1)} \left(\frac{1}{1+\alpha\mu}\right)^{1/\alpha} \left(\frac{\alpha\mu}{1+\alpha\mu}\right)^y.

This distribution has mean \(\mu\) and variance \(\mu+\alpha\mu^2\). In the GLM, write \(\eta_{gj}=x_j^\mathsf{T}\beta_g+\log(s_j)\), so that \(\mu_{gj}=\exp(\eta_{gj})\). Taking the logarithm of the probability and summing over the independent samples gives the gene-specific log-likelihood:

\ell_g(\beta_g;\alpha_g)= \sum_{j=1}^{n} \left[ \log\Gamma(Y_{gj}+1/\alpha_g) -\log\Gamma(1/\alpha_g) -\log\Gamma(Y_{gj}+1) +Y_{gj}\log(\alpha_g\mu_{gj}) -(Y_{gj}+1/\alpha_g)\log(1+\alpha_g\mu_{gj}) \right].

The expression is a function of \(\beta_g\) when \(\alpha_g\), the counts, the design, and the offsets are held fixed. Substituting \(\mu_{gj}=\exp(x_j^\mathsf{T}\beta_g+\log(s_j))\) makes the coefficient dependence explicit. In particular, the terms involving \(\beta_g\) can be written as

\ell_g(\beta_g;\alpha_g) = \sum_{j=1}^{n} \left[ Y_{gj}\eta_{gj} -(Y_{gj}+1/\alpha_g)\log(1+\alpha_g\exp(\eta_{gj})) \right] +C_g,

where \(C_g=\sum_j[\log\Gamma(Y_{gj}+1/\alpha_g)-\log\Gamma(1/\alpha_g)-\log\Gamma(Y_{gj}+1)+Y_{gj}\log(\alpha_g)]\) does not depend on \(\beta_g\). Thus maximizing the full log-likelihood over \(\beta_g\) is equivalent to maximizing the displayed coefficient-dependent part. The constant can be omitted in coefficient optimization, but not if an absolute log-likelihood value is required for comparing models or reporting the full likelihood.

Worked Example: A Single Negative Binomial Count

Suppose a hypothetical gene has count \(y=2\) in one sample, with model mean \(\mu=3\) and fixed dispersion \(\alpha=0.5\). Check the probability and log-likelihood contribution using the stated parameterization.

The size parameter is \(1/\alpha=2\). The probability mass function gives

P(Y=2)= \frac{\Gamma(4)}{\Gamma(2)\Gamma(3)} \left(\frac{1}{1+0.5(3)}\right)^2 \left(\frac{0.5(3)}{1+0.5(3)}\right)^2 =3(0.4)^2(0.6)^2=0.1728.

The log-likelihood contribution is \(\log(0.1728)\approx-1.7556\). Equivalently, the full log expression gives \(\log(3)+2\log(1.5)-4\log(2.5)\approx-1.7556\). The two calculations agree. This verifies the contribution for one count; a gene’s likelihood adds the corresponding log contributions across its samples.

What the Coefficient Optimization Uses

For fixed dispersion, the coefficient optimizer uses the observed counts \(Y_{gj}\), the design rows \(x_j^\mathsf{T}\), and the known offsets \(\log(s_j)\) through \(\eta_{gj}\). It maximizes the sum of the coefficient-dependent terms above. The gamma-function terms and \(Y_{gj}\log(\alpha_g)\) are constant in \(\beta_g\), so they do not change the maximizing coefficient vector.

The offset must nevertheless remain in \(\eta_{gj}\). Although its value is not estimated, it changes the fitted mean and therefore changes the likelihood as a function of the coefficients. Dropping it would fit a different model. Likewise, the dispersion affects the coefficient-dependent logarithm, so fixing \(\alpha_g\) does not mean it is irrelevant to coefficient estimation.

Worked Example: Comparing Two Coefficient Choices

Consider hypothetical counts \(Y=(3,8)\), design indicators \(x=(0,1)\), size factors \(s=(1,2)\), and fixed dispersion \(\alpha=0.2\). Use the model \(\eta_j=\beta_0+x_j\beta_1+\log(s_j)\). Compare the coefficient-dependent likelihood at (A) \(\beta_0=\log(4)\), \(\beta_1=\log(1.5)\), and (B) \(\beta_0=\log(4)\), \(\beta_1=0\). These are candidate values for comparison, not claimed estimates.

For (A), the fitted means are \(\mu_1=4\) and \(\mu_2=2(4)(1.5)=12\). Since \(1/\alpha=5\), the coefficient-dependent contribution for sample 1 is \(3\log(4)-(3+5)\log(1+0.2(4))=3\log(4)-8\log(1.8)\approx-0.5434\). For sample 2 it is \(8\log(12)-(8+5)\log(1+0.2(12))=8\log(12)-13\log(3.4)\approx3.9702\). The total coefficient-dependent log-likelihood is approximately \(3.4268\).

For (B), the means are \(\mu_1=4\) and \(\mu_2=8\). The first contribution remains \(-0.5434\); the second is \(8\log(8)-13\log(1+0.2(8))=8\log(8)-13\log(2.6)\approx4.2139\). The total is approximately \(3.6705\), which is higher than the value at (A) by about \(0.2437\). Because the omitted constants are identical for these two coefficient choices, the full log-likelihood has the same difference. This comparison alone does not locate the maximum; optimization must search over the allowed coefficient values.

Score and Curvature: A Short Justification

The score is the derivative of the log-likelihood with respect to the coefficients. Its form explains how the fitted means are adjusted during optimization. For one sample, hold the dispersion fixed and differentiate the coefficient-dependent contribution with respect to its linear predictor \(\eta\):

\frac{\partial}{\partial\eta} \left[y\eta-(y+1/\alpha)\log(1+\alpha e^\eta)\right] = y-(y+1/\alpha)\frac{\alpha e^\eta}{1+\alpha e^\eta} = \frac{y-\mu}{1+\alpha\mu}, \qquad \mu=e^\eta.

Since \(\partial\eta_{gj}/\partial\beta_g=x_j\), the score vector and its derivative are

U_g(\beta_g)= \frac{\partial\ell_g}{\partial\beta_g} = \sum_{j=1}^{n}x_j\frac{Y_{gj}-\mu_{gj}}{1+\alpha_g\mu_{gj}}, \qquad \frac{\partial^2\ell_g}{\partial\beta_g\partial\beta_g^\mathsf{T}} = -\sum_{j=1}^{n}x_jx_j^\mathsf{T} \frac{\mu_{gj}(1+\alpha_gY_{gj})}{(1+\alpha_g\mu_{gj})^2}.

The score balances observed and fitted counts, with each difference weighted by \(1/(1+\alpha_g\mu_{gj})\). The Hessian is negative semidefinite because each scalar multiplying \(x_jx_j^\mathsf{T}\) is nonnegative. Thus the log-likelihood is concave in \(\beta_g\) for fixed positive dispersion. If the design has full column rank and the maximum is finite, this concavity supports a well-defined coefficient maximum; rank deficiency can still prevent unique coefficients.

Worked Example: Checking a Stationary Coefficient

For a hypothetical intercept-only model, let \(Y=(2,4)\), \(s=(1,2)\), \(\alpha=0.25\), and \(\beta_0=\log(2)\). The design value is \(x_j=1\) for both samples, so the fitted means are \(\mu_1=2\) and \(\mu_2=4\).

The intercept score is \(U(\beta_0)=(2-2)/(1+0.25(2))+(4-4)/(1+0.25(4))=0+0=0\). The second derivative is \[ -\frac{2(1+0.25(2))}{(1+0.25(2))^2} -\frac{4(1+0.25(4))}{(1+0.25(4))^2} =-\frac{3}{2.25}-\frac{8}{4}\approx-3.333. \] It is negative, and the log-likelihood is concave in this intercept. Therefore this stationary value is a maximum. The coefficient is \(\log(2)\) on the log-mean scale; exponentiating it gives a fitted baseline mean multiplier of \(2\), before applying each sample’s size factor.

Assumptions, Checks, and Alternatives

The likelihood calculation relies on specific assumptions. Counts are nonnegative integers, conditional sample counts are independent for the gene, the negative binomial mean-variance relation is appropriate, and the size factors and design are treated as known. The dispersion is positive and held fixed in this coefficient likelihood. The design must contain enough independent information to identify the coefficients; otherwise multiple coefficient vectors can give the same fitted linear predictors.

In practice, check that the model matrix has no exact redundant columns, the offsets are finite and correspond to positive size factors, and the optimizer reports convergence to a finite solution. Nonconvergence, extremely large coefficients, or strong sensitivity to one sample can signal weak identification, a problematic design, or a poor fit to the count distribution. These checks do not by themselves establish that the negative binomial model is adequate, but they can reveal failures that make a reported coefficient unreliable.

The Poisson likelihood is a main alternative count likelihood: it assumes variance equal to the mean, whereas the negative binomial model allows variance \(\mu+\alpha\mu^2\). If the extra variation is material, the Poisson likelihood can represent the count variability poorly. A quasi-likelihood approach uses a specified mean-variance relationship without treating the full negative binomial probability mass function as the objective. These are different inferential frameworks; the negative binomial log-likelihood derived here specifically depends on the stated distribution and its fixed dispersion.

Common Mistakes and Limitations

Optimizing the wrong terms. The coefficient-independent terms may be dropped when maximizing over \(\beta_g\), but terms involving \(\eta_{gj}\), including the offset, cannot be dropped.
Confusing dispersion with a coefficient. The derivation holds \(\alpha_g\) fixed. It affects the likelihood and score, but this coefficient optimization does not itself estimate dispersion.
Using an inconsistent parameterization. Here the variance is \(\mu+\alpha\mu^2\), and the negative binomial size parameter is \(1/\alpha\). A formula using a different dispersion convention must be translated before use.
Calling any stationary point a unique estimate. Concavity does not repair a rank-deficient design or guarantee a finite maximum in every data configuration. Check identification and convergence.
Reading a larger kernel value as an absolute likelihood. The coefficient-dependent part is enough to compare coefficient values for the same gene and fixed dispersion, but it omits constants needed for the full log-likelihood.

Key takeaway: For fixed gene-wise dispersion, the negative binomial GLM estimates coefficients by maximizing \(\sum_j[Y_{gj}\eta_{gj}-(Y_{gj}+1/\alpha_g)\log(1+\alpha_g e^{\eta_{gj}})]\), where \(\eta_{gj}=x_j^\mathsf{T}\beta_g+\log(s_j)\). The offset is fixed, but it remains part of the likelihood through the fitted mean.

Check Your Understanding

Answer the questions, then compare your reasoning with the brief solutions.

For the stated negative binomial model, what is the variance of a count with mean \(\mu\) and dispersion \(\alpha\)?
Why can the gamma-function terms be omitted when maximizing over \(\beta_g\), but not when an absolute full log-likelihood is required?
With \(\alpha=0.5\), \(y=1\), and \(\mu=2\), calculate the probability mass \(P(Y=1)\).
In a two-sample model, the second sample has size factor \(s_2=3\). Where does this factor enter the coefficient likelihood?
What design condition helps ensure that the coefficient vector is identifiable?

Solutions. (1) The variance is \(\mu+\alpha\mu^2\). (2) Those terms depend on observed counts and fixed dispersion, not on the coefficients, so they do not change the maximizing coefficient; they are nevertheless part of the full log-likelihood. (3) The size is \(1/\alpha=2\), so \(P(Y=1)=\frac{\Gamma(3)}{\Gamma(2)\Gamma(2)}(1/2)^2(1/2)=2(0.25)(0.5)=0.25\). (4) It enters as the known offset \(\log(3)\) in \(\eta_2=x_2^\mathsf{T}\beta_g+\log(3)\), and hence changes \(\mu_2\). (5) The design matrix should have full column rank; exact redundant columns can make distinct coefficient vectors observationally equivalent.