Negative-Binomial Regression Likelihood and Score Equations

From a Count Distribution to a Regression Likelihood

A negative-binomial regression model connects a count $Y_i$ to covariates through its mean. Its likelihood is not the Poisson likelihood with an extra variance correction: the negative-binomial probability mass function changes the score itself. Here the size parameter is held fixed, and the regression coefficients are estimated by maximizing the resulting likelihood.

Use the mean–size parameterization: conditional on the covariates, $Y_i$ has mean $\mu_i$ and fixed size $\kappa>0$. Its variance is $\mu_i+\mu_i^2/\kappa$. The observations are assumed independent. Let $x_i$ be a $p$-dimensional covariate vector, including a constant component if the model has an intercept. The log link specifies $\log\mu_i=x_i^T\beta$, where $\beta$ is the vector of regression coefficients.

Definition: In fixed-size negative-binomial regression, $Y_1,\ldots,Y_n$ are independent, with $$ P(Y_i=y_i)= \frac{\Gamma(y_i+\kappa)}{\Gamma(\kappa)y_i!} \left(\frac{\kappa}{\kappa+\mu_i}\right)^\kappa \left(\frac{\mu_i}{\kappa+\mu_i}\right)^{y_i}, \qquad y_i=0,1,2,\ldots, $$ where $\kappa>0$ is fixed and $\mu_i=\exp(x_i^T\beta)$. The likelihood for $\beta$ is the product of these probabilities over $i$.

The fixed-$\kappa$ condition matters: terms involving $\kappa$ may be treated as constants when differentiating with respect to $\beta$. Estimating $\kappa$ is a separate problem. The derivations below concern the coefficient likelihood conditional on a specified value of $\kappa$.

Log-Likelihood and Coefficient Score

Independence makes the likelihood a product and the log-likelihood a sum. Substituting the model mean into the probability mass function gives

L(\beta)=\prod_{i=1}^n \frac{\Gamma(y_i+\kappa)}{\Gamma(\kappa)y_i!} \left(\frac{\kappa}{\kappa+\mu_i}\right)^\kappa \left(\frac{\mu_i}{\kappa+\mu_i}\right)^{y_i}, \qquad \mu_i=\exp(x_i^T\beta).

Taking logarithms and collecting the terms that depend on $\mu_i$ yields

\ell(\beta)=\sum_{i=1}^n \left[ \log\Gamma(y_i+\kappa)-\log\Gamma(\kappa)-\log(y_i!) +\kappa\log\kappa +y_i\log\mu_i-(y_i+\kappa)\log(\kappa+\mu_i) \right].

The first four terms inside each bracket do not depend on $\beta$, but retaining them gives the full log-likelihood. To differentiate the remaining terms, first use the linear predictor $\eta_i=x_i^T\beta$, so that $\mu_i=\exp(\eta_i)$ and $d\mu_i/d\eta_i=\mu_i$. The derivative of the $i$th log-likelihood contribution with respect to $\eta_i$ is

\frac{\partial\ell_i}{\partial\eta_i} =y_i-(y_i+\kappa)\frac{\mu_i}{\kappa+\mu_i} =\frac{\kappa(y_i-\mu_i)}{\kappa+\mu_i}.

Since $\partial\eta_i/\partial\beta=x_i$, the full score vector is

U(\beta)=\frac{\partial\ell(\beta)}{\partial\beta} =\sum_{i=1}^n x_i\frac{\kappa(y_i-\mu_i)}{\kappa+\mu_i} =X^T r(\beta), \qquad r_i(\beta)=\frac{\kappa(y_i-\mu_i)}{\kappa+\mu_i}.

The maximum-likelihood score equations set each of the $p$ components of this vector to zero. The residual-like quantity is not simply $y_i-\mu_i$: it is scaled by $\kappa/(\kappa+\mu_i)$. Thus, in general, fitting the model does not make the usual covariate-weighted raw residual sums equal to zero.

Worked Example: One Observation’s Probability and Log-Likelihood

Suppose an observation has $y=2$, fitted mean $\mu=3$, and fixed size $\kappa=2$. Its probability is

P(Y=2)= \frac{\Gamma(4)}{\Gamma(2)2!} \left(\frac{2}{5}\right)^2 \left(\frac{3}{5}\right)^2 =3\cdot\frac{4}{25}\cdot\frac{9}{25} =\frac{108}{625}=0.1728.

The corresponding log-likelihood contribution, including all terms, is

\ell_i =\log 3+2\log(2/5)+2\log(3/5) \approx 1.098612-1.832581-1.021651 =-1.755620.

The contribution to the score with respect to the linear predictor is $\kappa(y-\mu)/(\kappa+\mu)=2(2-3)/5=-0.4$. This is the derivative for this observation before multiplying by its covariate vector; the full regression score adds the contributions from all observations.

Worked Example: Score Equation in an Intercept-Only Model

Consider three independent counts $y=(0,1,4)$, with $\kappa=3$ and an intercept-only model. Every fitted mean is $\mu=\exp(\beta_0)$. The intercept score is

U_0(\beta_0) =\sum_{i=1}^3\frac{3(y_i-\mu)}{3+\mu} =\frac{3}{3+\mu}\{(0+1+4)-3\mu\}.

The factor $3/(3+\mu)$ is positive, so the score is zero precisely when $5-3\mu=0$. Therefore $\hat\mu=5/3$, and the corresponding coefficient is

\hat\beta_0=\log(5/3)\approx 0.510826.

This closed-form result is special to the intercept-only model: all observations have the same denominator $3+\mu$, so it factors out of the score. With varying covariates, the means and denominators differ, and this simplification generally disappears.

Curvature and Why Numerical Solution Is Usually Needed

The score equations are nonlinear in $\beta$, because every $\mu_i=\exp(x_i^T\beta)$ appears inside both a numerator and a denominator. Their curvature can be made explicit by differentiating the score contribution once more with respect to $\eta_i$:

\frac{\partial}{\partial\eta_i} \left\{\frac{\kappa(y_i-\mu_i)}{\kappa+\mu_i}\right\} =-\frac{\kappa\mu_i(\kappa+y_i)}{(\kappa+\mu_i)^2}.

Applying the chain rule to the score vector gives the Hessian matrix

H(\beta)=\frac{\partial^2\ell(\beta)}{\partial\beta\,\partial\beta^T} =-\sum_{i=1}^n \frac{\kappa\mu_i(\kappa+y_i)}{(\kappa+\mu_i)^2} x_i x_i^T.

Each scalar weight in this sum is positive: $\kappa>0$, $\mu_i>0$, and $\kappa+y_i>0$. Consequently, for every vector $v$,

v^T H(\beta)v =-\sum_{i=1}^n \frac{\kappa\mu_i(\kappa+y_i)}{(\kappa+\mu_i)^2} (x_i^Tv)^2\leq 0.

Result: The log-likelihood is concave in $\beta$. If the design matrix $X$, whose rows are $x_i^T$, has full column rank, then the Hessian is negative definite at every finite $\beta$, and the log-likelihood is strictly concave.

To verify the strict statement, take any nonzero $v$. Full column rank means $Xv\neq 0$, so at least one $x_i^Tv\neq 0$. Its squared term has a strictly positive weight, making the displayed quadratic form strictly negative. For concavity itself, the quadratic form is nonpositive for every $v$, which is exactly the Hessian criterion for a twice-differentiable concave function.

If a finite solution to $U(\beta)=0$ exists, concavity makes it a global maximizer. Indeed, along the line from a stationary point $\beta^*$ to any $\beta$, the log-likelihood is a concave one-dimensional function whose derivative at $\beta^*$ is $U(\beta^*)^T(\beta-\beta^*)=0$. Concavity implies its value cannot increase along that line. Under full column rank, strict concavity also rules out a second finite maximizer. These facts do not by themselves guarantee that a finite maximizer exists; the data and design still matter.

Except in special cases such as the intercept-only example, the nonlinear score equations generally have no algebraic solution. A standard approach is Newton–Raphson: at a current value $\beta^{(t)}$, solve $H(\beta^{(t)})\Delta^{(t)}=-U(\beta^{(t)})$, then set $\beta^{(t+1)}=\beta^{(t)}+\Delta^{(t)}$. Implementations may use step control if a full update does not improve the likelihood. The score and Hessian above give the quantities required for this coefficient-fitting step.

Worked Example: One Newton–Raphson Update

Take three observations with covariates $x_i^T=(1,z_i)$, $z=(0,1,2)$, counts $y=(0,1,4)$, and fixed size $\kappa=2$. Start at $\beta^{(0)}=(0,0)^T$. Then every fitted mean is $\mu_i=1$, and the score weights are $2(y_i-1)/3$. Thus

U(\beta^{(0)}) =\frac{2}{3} \left[ \begin{pmatrix}1\\0\end{pmatrix}(-1) +\begin{pmatrix}1\\1\end{pmatrix}(0) +\begin{pmatrix}1\\2\end{pmatrix}(3) \right] =\begin{pmatrix}4/3\\4\end{pmatrix}.

For the Hessian, the positive weights before the minus sign are $2(2+y_i)/9$, which are $4/9,6/9,12/9$. Therefore

-H(\beta^{(0)}) =\frac{1}{9} \left[ 4\begin{pmatrix}1&0\\0&0\end{pmatrix} +6\begin{pmatrix}1&1\\1&1\end{pmatrix} +12\begin{pmatrix}1&2\\2&4\end{pmatrix} \right] =\frac{1}{9} \begin{pmatrix}22&30\\30&54\end{pmatrix}.

The Newton increment solves $(-H)\Delta=U$. Multiplying both sides by $9$ gives $\begin{pmatrix}22&30\\30&54\end{pmatrix}\Delta=(12,36)^T$. The determinant is $22\cdot54-30\cdot30=288$, so

\Delta =\frac{1}{288} \begin{pmatrix}54&-30\\-30&22\end{pmatrix} \begin{pmatrix}12\\36\end{pmatrix} =\begin{pmatrix}-3/2\\3/2\end{pmatrix}.

The first updated coefficient vector is therefore $\beta^{(1)}=(-1.5,1.5)^T$. This is one iteration, not a claim that the score is now zero. Further updates evaluate the score and Hessian at the new coefficients and continue until a specified convergence criterion is met.

Common Pitfall and Practical Meaning

A frequent mistake is to copy the Poisson score equation and set $\sum_i x_i(y_i-\mu_i)=0$. For fixed-size negative-binomial regression, differentiation of the actual mass function instead gives $\sum_i x_i\kappa(y_i-\mu_i)/(\kappa+\mu_i)=0$. The scaling depends on the fitted mean, so it cannot generally be removed from the sum. Even though the mean model uses the same log link, changing the response distribution changes the likelihood and its estimating equations.

Another error is to differentiate with respect to $\kappa$ while treating it as fixed, or to omit the $\kappa$-dependent normalization terms when writing the full likelihood. Those terms are constant for coefficient estimation at fixed $\kappa$, but they are part of the likelihood itself and cannot be discarded when the dispersion parameter is to be estimated. The coefficient score derived here supplies the regression part of the fitting problem; it does not solve the separate dispersion estimation problem.

Key takeaway: With a log-linked mean and fixed size $\kappa$, the negative-binomial coefficient score is $\sum_i x_i\kappa(y_i-\mu_i)/(\kappa+\mu_i)$. Its nonlinear equations are typically solved numerically; the Hessian is negative semidefinite, and is negative definite when the design matrix has full column rank.

Check Your Understanding

Use the fixed-size likelihood and score equations developed above.

Write the mean $\mu_i$ in terms of $x_i$ and $\beta$ under the log link.
Why can the terms involving $\kappa$ be retained in the log-likelihood but treated as constants when differentiating with respect to $\beta$?
For one observation, what is the derivative of its log-likelihood with respect to $\eta_i$?
What feature of an intercept-only model makes its score equation simplify to a sample-mean equation?
Explain why full column rank of $X$ makes the Hessian negative definite at finite coefficient values.

Negative-Binomial Regression Likelihood and Score Equations

What You'll Learn