Why does my Poisson GLM with identity link converge, but a Negative Binomial GLM with identity link NOT converge?

Question

I am trying to fit a GLM to count data, with a single feature (also a count), and intercept term. Both the feature and response are nonnegative (but may be 0). I need to use identity link because I want to test the hypothesis that the coefficient in the linear model is 1.

I can fit a Poisson regression to the data no problem. However, as my data is definitely over-dispersed, I would like to fit a Negative Binomial regression instead. Yet when I do this, the optimization algorithm often fails to converge. I have noticed this behavior both in the glm.nb function from the MASS package in R, the vglm function from VGAM in R, as well as in the statsmodels package in Python, so this behavior seems unrelated to any particular software.

I have a number of datasets that I wish to fit these models on in sequence, and the same behavior is observed for most of them.

So my question is: Why would Poisson regression converge but Negative Binomial not? Given the nature of my data and the chosen link function, I understand the intercept and the coefficient both must be positive for the model to make any sense. However, I would think this would be a problem for both Poisson and Negative Binomial, but one works and the other does not.

My suspicion is that this is an issue related to whatever optimization algorithm these packages are using under the hood. But the fact that this happens across three packages in two different languages makes me doubt that this is the case. But if it is, is there anything I can do to rectify these issues, say, alter starting values or pass a different optimization method?

Editing to add a minimum reproducible example. The raw data is unavailable to share but I have taken a small sample. This converges with Poisson, but not with Negative Binomial, when using the statsmodels package in Python. Note the data point with an x value of 0 and a very large y.

import statsmodels.api as sm
import numpy as np

X = np.array([60, 16, 32, 0, 1148, 96, 16, 32, 208, 2, 23, 30, 60, 340, 16, 132, 51, 350, 4, 0])
X = np.vstack((X, [1]*X.shape[0])).T # add intercept
y = np.array([48, 14, 21, 1, 2779, 96, 16, 24, 79, 1, 12, 30, 1, 135, 16, 223, 37, 154, 4, 3279])

poiReg = sm.GLM(y, X, family = sm.families.Poisson(link = sm.families.links.identity)).fit()
poiReg.summary()

nbReg = sm.GLM(y, X, family = sm.families.NegativeBinomial(link = sm.families.links.identity)).fit()
nbReg.summary()

Also note that while this model fails to converge with statsmodels, it actually does with glm.nb in R. However, the full dataset does not, and I do not have the time to continue trying sample permutations until I find one that fails in all cases.

Answer 1

The linear link function does not impose non-negativity constraints on the prediction, i.e. conditional mean of the distribution.

It is possible that for some values of the explanatory variables, the predicted mean is negative even if the optimization has converged. This could be the case for example if we use OLS to estimate the model.

In this case some parameters of the distribution, loglikelihood, implied variance function or similar, will not be in a valid range and the evaluation of those functions will break or raise an exception. Then, the family-link combination cannot be estimated.

In other cases, it is possible that the nonnegativity constraints are not binding at the optimum, i.e. the maximum likelihood estimate. In those cases improving starting values can help in achieving convergence.

Specifially, in GLM families different from gaussian, the implied variance is a function of the mean, negative predicted mean implies therefore a negative variance which will in most cases cause numerical problems.