如何在Python statsmodels中实现类似R的二项式回归？

你遇到的PatsyError是因为statsmodels公式语法不支持(succ, fail)这种写法——Patsy会把这两个变量解析成一个长度为2的序列，和你的6行数据行数不匹配，导致报错。以下是两种可行的实现方式：

方法1：使用公式接口（适配statsmodels语法） #

statsmodels的Binomial家族支持两种响应变量形式：单个0/1观测，或者成功次数占总试验次数的比例，此时需要通过weights参数传入总试验次数（成功数+失败数）。代码如下：

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
binom = sm.families.Binomial()

tbTest = pd.DataFrame({
  'x'    : [1, 2, 3, 4, 5, 6],
  'succ' : [5, 4, 3, 2, 1, 0],
  'fail' : [1, 0, 1, 2, 2, 3],
})

# 计算总试验次数
tbTest['total'] = tbTest['succ'] + tbTest['fail']
# 用比例作为响应变量，指定weights
lrTest = smf.glm('succ / total ~ x', tbTest, family=binom, weights=tbTest['total']).fit()
print(lrTest.summary())

也可以用更简洁的写法，直接传入成功次数作为响应变量，statsmodels会结合weights自动计算比例：

lrTest = smf.glm('succ ~ x', tbTest, family=binom, weights=tbTest['total']).fit()

方法2：直接使用GLM类（对应R的矩阵输入形式） #

如果想完全复刻R中用矩阵传入成功数和失败数的写法，可以绕过公式接口，直接使用sm.GLM，传入二维数组作为响应变量：

import pandas as pd
import statsmodels.api as sm
binom = sm.families.Binomial()

tbTest = pd.DataFrame({
  'x'    : [1, 2, 3, 4, 5, 6],
  'succ' : [5, 4, 3, 2, 1, 0],
  'fail' : [1, 0, 1, 2, 2, 3],
})

# 构造带截距项的设计矩阵
X = sm.add_constant(tbTest['x'])
# 构造响应变量矩阵：每行对应[成功数, 失败数]
y = tbTest[['succ', 'fail']].values

model = sm.GLM(y, X, family=binom)
lrTest = model.fit()
print(lrTest.summary())

两种方法的计算结果和你在R中执行的glm(matrix(c(succ, fail), ncol=2) ~ x, data=tbTest, family="binomial")完全一致。

内容的提问来源于stack exchange，提问作者Igor F.