Profiling the Negative-Binomial Dispersion Parameter

Why Profile the Size Parameter?

In negative-binomial regression, the size parameter controls how much variation the model allows beyond the mean. A small size corresponds to greater overdispersion, while a large size corresponds to less. The regression coefficients and size are estimated together: changing the size changes the likelihood, and the best-fitting regression coefficients can change with it. A profile likelihood records the best likelihood available at each fixed size, after the coefficients have been re-estimated.

This distinction matters in practice. Evaluating the likelihood over sizes while holding the regression coefficients at one fit does not generally produce the profile. It instead traces a slice through the joint likelihood. The profile is the relevant one-dimensional summary when studying how the fitted model supports different finite size values.

Model and Profile Definition

Definition: Suppose $Y_1,\ldots,Y_n$ are independent negative-binomial counts with means $\mu_i>0$ and common size $k>0$, using the parameterization $$ P(Y_i=y_i)=\frac{\Gamma(y_i+k)}{\Gamma(k)y_i!} \left(\frac{k}{k+\mu_i}\right)^k \left(\frac{\mu_i}{k+\mu_i}\right)^{y_i}, \qquad y_i=0,1,2,\ldots. $$ In a log-link regression, $\log\mu_i=x_i^{\mathsf T}\beta+o_i$, where $x_i$ is the covariate vector, $\beta$ is the regression coefficient vector, and $o_i$ is a known offset. The log-likelihood is $\ell(\beta,k)$. The profile log-likelihood for $k$ is $$ \ell_p(k)=\sup_{\beta\in\mathcal B}\ell(\beta,k),\qquad k>0, $$ where $\mathcal B$ is the allowed coefficient space. When the supremum is attained, profiling at $k$ means finding a maximizing coefficient vector $\hat\beta(k)$ and evaluating $\ell(\hat\beta(k),k)$.

An equivalent overdispersion parameterization is $\alpha=1/k$, for which the NB2 variance is $\operatorname{Var}(Y_i)=\mu_i+\alpha\mu_i^2$. Profiling in $\alpha$ is the same calculation under the one-to-one transformation $k=1/\alpha$, with the horizontal axis relabelled. The maximized log-likelihood values do not change under that relabelling.

For one observation, omitting no terms that depend on $k$, the log-likelihood contribution is $$ \ell_i(\beta,k)=\log\Gamma(y_i+k)-\log\Gamma(k)-\log(y_i!) +k\log k-k\log(k+\mu_i)+y_i\log\mu_i-y_i\log(k+\mu_i). $$ For fixed $k$, the coefficient score under the log link is $$ \frac{\partial\ell}{\partial\beta} =\sum_{i=1}^n x_i\frac{k(y_i-\mu_i)}{k+\mu_i}. $$ Thus, at each candidate $k$, the coefficient fit solves a score equation that depends on $k$. The size score at fixed $\beta$, writing $\psi$ for the digamma function, is $$ \frac{\partial\ell}{\partial k} =\sum_{i=1}^n\left[ \psi(y_i+k)-\psi(k)+\log\left(\frac{k}{k+\mu_i}\right) +1-\frac{k+y_i}{k+\mu_i} \right]. $$ These formulas can be used in numerical optimization, but the profile score is evaluated at the coefficient fit corresponding to that same $k$.

Worked Profile Calculations

Worked Example: An Intercept-Only Fit at Three Sizes

Observe counts $0$ and $2$ from an intercept-only negative-binomial model, so both observations have mean $\mu$. At any fixed $k$, the maximizing mean is $\hat\mu=1$. To verify this, the derivative of the log-likelihood with respect to $\eta=\log\mu$ is $$ \frac{\partial\ell}{\partial\eta} =\sum_{i=1}^2\frac{k(y_i-\mu)}{k+\mu} =\frac{k(2-2\mu)}{k+\mu}, $$ which is zero at $\mu=1$, positive for $\mu<1$, and negative for $\mu>1$.

For $\mu=1$, the zero-count probability is $P(Y=0)=(k/(k+1))^k$, and $$ P(Y=2)=\frac{k(k+1)}{2}\left(\frac{k}{k+1}\right)^k\left(\frac{1}{k+1}\right)^2. $$ Multiplying these gives the profile likelihood for the two observations. At $k=1$, the two probabilities are $1/2$ and $1/8$, so $L_p(1)=1/16$. At $k=2$, they are $4/9$ and $4/27$, so $L_p(2)=16/243$. At $k=4$, they are $256/625$ and $512/3125$, so $L_p(4)=131072/1953125$.

The corresponding profile log-likelihood values are approximately $-2.773$, $-2.720$, and $-2.701$. They increase across these three candidate sizes, but this finite comparison alone does not establish where the profile has its maximum. It does show why every value must use the coefficient fit at that size; here the fitted mean happens to remain 1.

Worked Example: Re-estimating Regression Coefficients by Group

Consider four independent counts, with two observations in each of two covariate groups. The observed counts are $0,2$ in group A and $1,3$ in group B. Fit an intercept and a group indicator, so each group has its own mean. At every fixed $k$, the group means that maximize the likelihood are the within-group averages: $\hat\mu_A=1$ and $\hat\mu_B=2$. This follows by applying the log-mean score calculation from the previous example separately within each group.

For group A, the product of probabilities at $k=1,2,4$, respectively, is $1/16$, $16/243$, and $131072/1953125$, as calculated above. For group B, at mean 2, the probabilities for counts 1 and 3 give the following products:

At $k=1$, $P(Y=1)=2/9$ and $P(Y=3)=8/81$, so the product is $16/729$. At $k=2$, the probabilities are $1/4$ and $1/8$, so the product is $1/32$. At $k=4$, they are $64/243$ and $320/2187$, so the product is $20480/531441$.

Independence makes the full profile likelihood the product of the group-specific products. Consequently, $$ L_p(1)=\frac{1}{16}\frac{16}{729}=\frac{1}{729},\qquad L_p(2)=\frac{16}{243}\frac{1}{32}=\frac{1}{486}, $$ and $$ L_p(4)=\frac{131072}{1953125}\frac{20480}{531441}. $$ The coefficient estimates correspond to group means 1 and 2 at each candidate size, or to regression coefficients $\hat\beta_0(k)=0$ and $\hat\beta_1(k)=\log 2$. This example happens to have the same fitted means at all three sizes because the regression is saturated for the two groups. In a less saturated regression, the fitted means will generally vary with $k$, making repeated refitting essential.

Profile Score and Curvature

Theorem: Suppose $\ell(\beta,k)$ is twice continuously differentiable near $k_0>0$, and there is a unique interior maximizing coefficient vector $\hat\beta(k)$ near $k_0$. Suppose $\ell_{\beta\beta}(\hat\beta(k_0),k_0)$ is nonsingular and negative definite, and the implicit function theorem applies to the coefficient score equation. Then $$ \ell_p'(k_0)=\ell_k(\hat\beta(k_0),k_0), $$ and $$ \ell_p''(k_0)=\ell_{kk} -\ell_{k\beta}\ell_{\beta\beta}^{-1}\ell_{\beta k}, $$ with all derivatives on the right evaluated at $(\hat\beta(k_0),k_0)$.

Proof. The coefficient score vanishes at the profile fit: $\ell_\beta(\hat\beta(k),k)=0$. By the chain rule, $$ \ell_p'(k)=\ell_k+\ell_\beta^{\mathsf T}\hat\beta'(k)=\ell_k, $$ which proves the first identity. Differentiate the score equation with respect to $k$: $$ \ell_{\beta\beta}\hat\beta'(k)+\ell_{\beta k}=0, \qquad \hat\beta'(k)=-\ell_{\beta\beta}^{-1}\ell_{\beta k}. $$ Differentiating $\ell_p'(k)=\ell_k(\hat\beta(k),k)$ now gives $$ \ell_p''(k)=\ell_{kk}+\ell_{k\beta}\hat\beta'(k) =\ell_{kk}-\ell_{k\beta}\ell_{\beta\beta}^{-1}\ell_{\beta k}. $$ The stated conditions ensure the local maximizing branch is differentiable and these derivatives exist. $\square$

The first identity is the envelope property: the profile slope is the size score evaluated at the re-estimated coefficients. For curvature, define the observed-information blocks $I_{\beta\beta}=-\ell_{\beta\beta}$, $I_{\beta k}=-\ell_{\beta k}$, and $I_{kk}=-\ell_{kk}$. The negative profile curvature is $$ -\ell_p''(k)=I_{kk}-I_{k\beta}I_{\beta\beta}^{-1}I_{\beta k}. $$ When the joint observed information is positive definite, this is a Schur complement and is positive. It is the curvature remaining after allowing the regression coefficients to adjust, rather than the curvature obtained by freezing them.

Worked Example: How Coefficient Adjustment Changes Curvature

Suppose a fitted model has one regression coefficient and one size parameter, and the observed-information blocks at a particular fit are $I_{\beta\beta}=4$, $I_{\beta k}=I_{k\beta}=1$, and $I_{kk}=2$. The joint information matrix is $$ I=\begin{pmatrix}4&1\\1&2\end{pmatrix}. $$ It is positive definite because its first leading entry is 4 and its determinant is $4(2)-1(1)=7>0$. If the coefficient were held fixed, the negative curvature in $k$ would be $I_{kk}=2$. After profiling, it is $$ I_{kk}-I_{k\beta}I_{\beta\beta}^{-1}I_{\beta k} =2-(1)(1/4)(1)=\frac{7}{4}. $$ Allowing the coefficient to respond to changes in $k$ reduces the local curvature from 2 to 1.75. This calculation is local: it describes curvature at the fitted point, not the profile’s shape over its entire parameter range.

Finding and Interpreting the Profile Maximum

A practical profile calculation evaluates a grid of positive sizes, maximizes $\ell(\beta,k)$ over $\beta$ at each grid point, and then refines the search near the best values. A reliable implementation uses the previous nearby coefficient fit as a starting value when appropriate, but checks convergence and the coefficient score at each size. The size maximizing $\ell_p(k)$, if attained in the parameter space being considered, is the joint maximum-likelihood estimate of size; its associated coefficient fit is the regression estimate at that size.

Near an interior maximum $k_0$, the profile score is zero and the profile log-likelihood is locally concave when $\ell_p''(k_0)<0$. Greater negative curvature indicates a sharper local peak, so the likelihood changes more rapidly for nearby size values. A broad, flat peak indicates that the data distinguish nearby sizes less strongly. Curvature is not itself a standard error or a test result; converting a profile into inferential statements requires the relevant likelihood-inference framework and its conditions.

Common Pitfall: Holding the Coefficients Fixed

The incorrect argument is: “Fit the regression once, keep $\hat\beta$ fixed, and vary $k$; that curve is the profile likelihood.” That curve is only $\ell(\hat\beta,k)$. By definition, the profile is $\sup_\beta\ell(\beta,k)$, so at every candidate $k$ the coefficients must be optimized again. The fixed-coefficient curve may lie below the profile and can have a different slope and curvature.

A related numerical pitfall is to round log-likelihood values before calculating differences or slopes. For example, a small difference formed from rounded values can give a noticeably inaccurate secant slope. Keep full precision for optimization and intermediate calculations, then round only the reported results.

Key takeaway: A negative-binomial size profile is the maximized log-likelihood over regression coefficients at each fixed positive size. Its slope uses the size score at that refit, and its curvature accounts for coefficient adjustment through a Schur complement.

Check Your Understanding

Answer these questions using the definitions and calculations above.

In a negative-binomial regression, what is re-estimated at each fixed value of the size parameter when forming the profile log-likelihood?
For the counts $0,2$ in an intercept-only fit, verify that the fitted mean is 1 at any fixed positive size by examining the log-mean score.
Explain why $\ell(\hat\beta,k)$, with one fixed coefficient estimate, is generally not the profile log-likelihood.
State the regularity conditions used to derive the profile-curvature formula, and identify the adjustment term.
What does a broad profile maximum suggest about the local information for the size parameter, and why is that not by itself a test result?