为何statsmodels的NegativeBinomialP回归系数与R的计算结果不一致？

1. 未手动添加截距项 #

R的glm.nb拟合时默认自动为回归方程生成截距项，但statsmodels的所有回归模型默认不会自动添加截距，需要手动调用sm.add_constant()处理自变量矩阵。你当前的代码没有加截距，自然和R的输出系数对应不上。

2. 三分类分类变量编码错误 #

你把prog这个三分类无序变量直接编码为1、2、3作为连续变量输入模型，这是完全错误的。R在拟合时会自动将字符型分类变量转换为哑变量（虚拟变量），默认以第一个水平为参考组，原练习里prog的三个水平是Academic、Vocational、General，正确的做法是生成2个哑变量输入模型，而非当成连续数值处理。

3. 可选优化：数据读取方式不规范 #

手动逐行读取文件容易出现格式解析错误，建议直接用pandas的read_csv读取数据集，避免手动拆分字段出现的异常。

修正后可复现R结果的代码 #

import numpy as np
import pandas as pd
import statsmodels.api as sm
from statsmodels.discrete.discrete_model import NegativeBinomialP

# 读取数据
df = pd.read_csv("sampleNBdata2.dat")
# 处理分类变量prog，生成哑变量，默认以General为参考组和R逻辑保持一致
df = pd.get_dummies(df, columns=['prog'], drop_first=True)
# 构造自变量：math + prog的两个哑变量 + 手动添加截距项
exog = sm.add_constant(df[['math', 'prog_Academic', 'prog_Vocational']])
# 因变量
endog = df['daysabs']

# 拟合NB2模型，若bfgs收敛精度不够可换用newton方法
model_nb = NegativeBinomialP(endog, exog, p=2)
res_nb = model_nb.fit(method='bfgs', maxiter=5000)
print(res_nb.summary())

内容的提问来源于stack exchange，提问作者ankit agrawal