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.