I need to estimate the log-likelihood of the negative binomial regression. I mean full log-likelihood, including the dispersion parameter.

The problem is: when I start LBFGS, BFGS, or gradient descent procedures, they estimate the dispersion parameter all differently depending on the starting point. I use this for the derivatives: https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Negative_Binomial_Regression.pdf

Is the negative binomial regression log-likelihood actually concave?

PS. I know that glm.nb in R estimates the dispersion in steps - first it estimates the linear coefficients with the Newton method, then it estimates the dispersion, then it comes back to linear coefficients, and so forth. This seems to give the same result not depending on the starting point. But I do not understand why it is done this way.

PPS. Links are welcome.


2 Answers 2


It is relatively simple to write the log-density of the negative binomial distribution in terms of its mean, and then use this to get a log-likelihood expression for the negative binomial GLM. For all values $y=0,1,2,...$ the log-density of the negative binomial distribution is:

$$\log f(y|r,\theta) = \log {y+r-1 \choose y} + r \log (1-\theta) + y \log (\theta).$$

If we let $\mu \equiv \mathbb{E}(Y)$ denote the mean of the distribution then we have the relationships:

$$\mu = r \cdot \frac{\theta}{1-\theta} \quad \quad \quad \theta = \frac{\mu}{r+\mu}.$$

Thus, we can re-parameterise the log-density in terms of its mean parameter to obtain the log-likelihood for a single data point:

$$\ell_y(r,\mu) \equiv \log f(y|r,\mu) = \log {y+r-1 \choose y} + r \log (r) + y \log (\mu) - (r+y) \log (r+\mu).$$

The derivatives with respect to $\mu$ are:

$$\begin{aligned} \frac{d\ell_y}{d \mu} (r,\mu) &= \frac{y}{\mu} - \frac{r+y}{r+\mu}, \\[6pt] \frac{d^2 \ell_y}{d \mu^2} (r,\mu) &= -\frac{y}{\mu^2} + \frac{r+y}{(r+\mu)^2}. \\[6pt] \end{aligned}$$

Letting $\hat{\mu}$ be a critical point of the log-density, we obtain:

$$\frac{d^2 \ell_y}{d \mu^2} (r,\hat{\mu}) = -\frac{y}{\hat{\mu}} \Bigg[ \frac{1}{\hat{\mu}} + \frac{1}{r+\hat{\mu}} \Bigg] <0.$$

This shows that all the log-likelihood for a single data point has a single maximising value for its mean parameter, and so it is quasi-concave with respect to the mean parameter. Unfortunately, this simple result gets more muddied when we transition to the log-likelihood for multiple data points in the negative binomial GLM.

Negative binomial GLM: In the negative binomial GLM we have data $Y_i \sim \text{NegBin}(r, \mu_i)$ where the value $\mu_i$ is determined from a corresponding explanatory variable $\mathbf{x}_i$ using a link function. Written in general terms, using the mean vector $\boldsymbol{\mu}$ the log-likelihood of the data $\mathbf{y}=(y_1,...,y_n)$ is:

$$\ell_\mathbf{y} (r, \boldsymbol{\mu}) = \sum_i \log {y_i+r-1 \choose y_i} + n r \log (r) + \sum_i y_i \log (\mu_i) - \sum_i (r+y_i) \log (r+\mu_i).$$

In the standard negative binomial GLM we use the link function $\log \mu_i(\boldsymbol{\beta}) = \sum_k \beta_k \cdot x_{i,k}$, which gives the likelihood function:

$$\begin{aligned} \ell_\mathbf{y} (r, \boldsymbol{\beta}) &= \sum_i \log {y_i+r-1 \choose y_i} + n r \log (r) \\ &\quad + \sum_i y_i \sum_k \beta_k \cdot x_{i,k} - \sum_i (r+y_i) \log (r+\mu_i(\boldsymbol{\beta})). \end{aligned}$$

The derivatives of this function with respect to the model coefficients are:

$$\begin{aligned} \frac{\partial \ell_\mathbf{y}}{\partial \beta_k} (r, \boldsymbol{\beta}) &= \sum_i y_i \cdot x_{i,k} - \sum_i (r+y_i) \cdot x_{i,k} \cdot \frac{\mu_i(\boldsymbol{\beta})}{r+\mu_i(\boldsymbol{\beta})}, \\[6pt] \frac{\partial^2 \ell_\mathbf{y}}{\partial \beta_k \partial \beta_l} (r, \boldsymbol{\beta}) &= - \sum_i (r+y_i) \cdot \frac{\mu_i(\boldsymbol{\beta})}{r+\mu_i(\boldsymbol{\beta})} \Bigg( 1- \frac{1}{r+\mu_i(\boldsymbol{\beta})} \Bigg) \cdot x_{i,k} \cdot x_{i,l}. \\[6pt] \end{aligned}$$

As you can see, these derivatives are more complicated than in the case of a single data point. This is owing to the fact that we now have $n$ data points, and we are looking at derivatives with respect to the model coefficients, rather than the mean value.

  • $\begingroup$ True, but even for many points, the log-likelihood is still concave. Which can be seen from writing neg.bin regression in the exponential form if the parameter $r$ is preset (I use the NB2 parameterization with variance $\mu + \alpha \mu^2$). However, if the $r$ is also estimated, then the function still has one maximum, but it is not concave anymore. So it breaks all 2nd order methods if they accidentally fall into the wrong region. It is a shame such a basic fact is not explained anywhere, and all the books talk about the optimization with Newton method, overlooking this basic fact. $\endgroup$ Jan 20, 2020 at 22:47

The answer is: NO, negative binomial regression is NOT CONCAVE when estimated with dispersion.

I do not know how to show this mathematically (since negative binomial WITH dispersion parameter does not belong to the exponential family), but I simply found a line in my parameter space, which is not concave. Here it is:

enter image description here

Target is the log-likelihood of the negative binomial model with 28 parameters, and $x$ ranges from 0 to 1 along the x-axis. I do not provide the numbers, but you can easily see from the picture, that the log-likelihood is not concave along the x-axis (concavity changes into convexity around 0.01).

In practice it means that during Newton method optimization, some of the Hessians become indefinite, which breaks the method (it tries to find the dispersion at infinity as seen from the picture). This is why glm.nb uses coordinate descent.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.