为啥 NegativeBinomialP 与 R 相比给出不同的系数？

Posted 2023-02-14

【中文标题】为啥 NegativeBinomialP 与 R 相比给出不同的系数？

【英文标题】：Why NegativeBinomialP gives different coefficients compared to R?为什么 NegativeBinomialP 与 R 相比给出不同的系数？

【发布时间】：2022-01-08 12:48:06

【问题描述】：

在 python 中重复以下 R 练习以达到相同的结果，我几乎没有困难。我错过了什么？

R 练习 https://stats.idre.ucla.edu/r/dae/negative-binomial-regression/

数据链接 https://www.dropbox.com/s/mz4stp72eco3rfq/sampleNBdata2.dat?dl=0

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import statsmodels.distributions.discrete as distr
from statsmodels.discrete.discrete_model import NegativeBinomialP, NegativeBinomial, Poisson, GeneralizedPoisson
from statsmodels.discrete.count_model import (ZeroInflatedNegativeBinomialP, ZeroInflatedPoisson,
                                              ZeroInflatedGeneralizedPoisson)
import statsmodels.discrete._diagnostics_count as dia

import statsmodels.api as sm


f=open('sampleNBdata2.dat')
id=[]
gender=[]
math=[]
daysabs=[]
prog=[]
x=[]
f.readline()
d=
d['Academic']=1
d['Vocational']=2
d['General']=3

for line in f:
    l=line.split(',')
    id.append(l[1])
    gender.append(l[2])
    math.append(l[3]) #independent
    daysabs.append(int(l[4])) #dependent y
    prog.append(l[5]) #independent
    #x.append([int(l[3]),d[l[5]], ] )
    x.append([int(l[3]),int(l[5]), ] )

print(x,daysabs)
endog=np.array(daysabs)
exog=np.array(x)


print("endog",endog.shape)
print("exog",exog.shape)



#model_nb = NegativeBinomial(endog, exog, loglike_method='nb2')
model_nb = NegativeBinomialP(endog, exog, p=2)

res_nb = model_nb.fit(method='bfgs', maxiter=5000, maxfun=5000)

print(endog)
print(exog)
print(res_nb.summary())

Python 输出 R输出

【问题讨论】：

看起来您没有向 statsmodels 中的 exog 添加拦截。它不会自动添加，除非使用公式。 prog 在 R 中是分类的，而您将其转换为数字 AFAICS。使用 pandas 和公式，则默认行为类似于 R。

【参考方案1】：

以下代码以几乎相似的系数再现 R 的结果。

df=pd.read_csv('sampleNBdata.dat')
data=pd.concat((df,pd.get_dummies(df['prog'],drop_first=False)),axis=1)
endog=data['daysabs']
data['intercept'] = 1
exog=data.drop(['prog','daysabs','id','gender','Unnamed: 0','General'],axis=1)
model_nb = NegativeBinomialP(endog, exog, p=2)
res_nb = model_nb.fit(method='bfgs', maxiter=5000, maxfun=5000)
print(res_nb.summary())