I'm trying to implement a negative binomial regression model in R which was originally implemented in Python with statsmodels.
In some cases, the estimated coefficients are very similar:
Python code:
import statsmodels.api as sm
import numpy as np
response = [17, 18, 10, 9, 8, 5, 6, 5, 15351, 9637, 9981, 9306, 16752, 11993, 13622, 9800]
design = np.array(
[[1, 1, 1, 0],
[1, 1, 1, 0],
[1, 1, 1, 0],
[1, 1, 1, 0],
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 0, 0, 0],
[1, 0, 1, 1],
[1, 0, 1, 1],
[1, 0, 1, 1],
[1, 0, 1, 1],
[1, 0, 0, 1],
[1, 0, 0, 1],
[1, 0, 0, 1],
[1, 0, 0, 1]])
theta = 0.00792199771866427
sm.GLM(response, design, family=sm.families.NegativeBinomial(alpha=theta)).fit().params
Which gives:
array([ 1.79175947, 0.97496014, -0.16402993, 7.68415156])
And the equivalent model in R:
library(MASS)
response = c(17, 18, 10, 9, 8, 5, 6, 5, 15351, 9637, 9981, 9306, 16752, 11993, 13622, 9800)
design = as.data.frame(matrix(c(1, 1, 1, 1, 1, 1, 1, 1,1 , 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 1, 1, 1, 1, 1, 1, 1), ncol = 4, byrow = F))
theta = 0.00792199771866427
coef(glm(response ~. +0, design, family = negative.binomial(theta)))
Which gives:
V1 V2 V3 V4
1.7917595 0.9749601 -0.1640299 7.6841516
So for these two models, the estimated coefficients are very similar, down to the second decimal place. I am also fitting a reduced model dropping the second column of the model matrix, however, here the estimated coefficients are quite different between R and Python.
Python
sm.GLM(response, design[:, [0,2,3]], family=sm.families.NegativeBinomial(alpha=theta)).fit().params
array([ 2.32804838, -0.10095997, 7.11684136])
R
coef(glm(response ~. +0, design[, c(1,3,4)] , family = negative.binomial(theta)))
V1 V3 V4
2.0650897 0.3232355 7.1965690
Why does this occur? I have noticed this characteristic for a number of different feature sets.