从R迁移至Python实现含偏移量、案例权重与分类变量交互项的Gamma回归，寻求与R结果一致的解决方案

我在把复杂的GLM模型从R迁移到Python时遇到了点小麻烦，核心诉求是让两种语言的结果在本质上保持一致。我已经借助ChatGPT准备了一个可复现的例子，结果却发现Python的statsmodels包处理分类变量交互项的方式和R不一样（下面会贴出两边的输出对比）。

我现在差点就要手动硬编码模型矩阵，让它和R里的特征完全一致了，但还是不想太早放弃，相信肯定有更简单的方法能让Python输出和R（以及SAS）几乎一致的结果。

R实现 #

代码 #

# 创建包含3个变量、偏移量和案例权重的样本数据集
dt <- data.frame(target = c(1.5, 2.0, 3.0, 4.0, 2.5, 3.5, 4.5, 3.0, 4.5, 5.5, 4.0, 3.5, 2.5, 3.0, 2.0, 1.5, 2.5, 3.5, 4.5, 3.0),
var1 = factor(c("A", "A", "B", "B", "C", "C", "D", "D", "A", "A", "B", "B", "C", "C", "D", "D", "A", "A", "B", "B")),
var2 = factor(c("X", "Y", "X", "Y", "X", "Y", "X", "Y", "X", "Y", "X", "Y", "X", "Y", "X", "Y", "X", "Y", "X", "Y")),
var3 = factor(c("M", "N", "N", "M", "N", "M", "M", "N", "N", "M", "M", "N", "M", "N", "N", "M", "M", "N", "N", "M")),
offset = c(1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2),
case_weights = c(1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5))

# 拟合包含偏移量、案例权重和分类变量交互项的GLM模型
fit <- glm(formula = as.formula("target ~ var1 + var2:var3"),
offset = offset,
weights = case_weights,
family = Gamma(link="log"),
data = dt)

summary(fit)

输出 #

Call:
glm(formula = as.formula("target ~ var1 + var2:var3"), family = Gamma(link = "log"),
    data = dt, weights = case_weights, offset = offset)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.0494  -0.5148  -0.1940   0.4460   1.1210  

Coefficients: (1 not defined because of singularities)
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) -0.30812    0.34302  -0.898  0.38540   
var1B       -0.87547    0.38038  -2.302  0.03850 * 
var1C       -0.12718    0.44326  -0.287  0.77870   
var1D       -1.18807    0.40693  -2.920  0.01190 * 
var2X:var3M -0.04003    0.40421  -0.099  0.92260   
var2Y:var3M  0.13896    0.40131   0.346  0.73470   
var2X:var3N -0.13315    0.40397  -0.330  0.74700   
var2Y:var3N       NA         NA      NA       NA   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.5822375)

    Null deviance: 13.6545  on 19  degrees of freedom
Residual deviance:  6.8235  on 13  degrees of freedom
AIC: 125.64

Number of Fisher Scoring iterations: 14

Python实现 #

代码 #

import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 创建包含3个变量、偏移量和案例权重的样本数据集
dt = pd.DataFrame({'target': [1.5, 2.0, 3.0, 4.0, 2.5, 3.5, 4.5, 3.0, 4.5, 5.5, 4.0, 3.5, 2.5, 3.0, 2.0, 1.5, 2.5, 3.5, 4.5, 3.0],
'var1': ['A', 'A', 'B', 'B', 'C', 'C', 'D', 'D', 'A', 'A', 'B', 'B', 'C', 'C', 'D', 'D', 'A', 'A', 'B', 'B'],
'var2': ['X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y', 'X', 'Y'],
'var3': ['M', 'N', 'N', 'M', 'N', 'M', 'M', 'N', 'N', 'M', 'M', 'N', 'M', 'N', 'N', 'M', 'M', 'N', 'N', 'M'],
'offset': [1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2, 1, 2, 3, 2],
'case_weights': [1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5, 1, 1.5, 2, 1.5]})

# 拟合包含偏移量、案例权重和分类变量交互项的GLM模型
fit = smf.glm(formula='target ~ var1 + var2:var3',
offset=dt['offset'],
data=dt,
family=sm.families.Gamma(link=sm.families.links.log()),
freq_weights=dt['case_weights']).fit()

print(fit.summary())

输出 #

Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:                 target   No. Observations:                   20
Model:                            GLM   Df Residuals:                       23
Model Family:                   Gamma   Df Model:                            6
Link Function:                    log   Scale:                         0.32910
Method:                          IRLS   Log-Likelihood:                -55.981
Date:                Thu, 20 Apr 2023   Deviance:                       6.8235
Time:                        15:48:56   Pearson chi2:                     7.57
No. Iterations:                    36   Pseudo R-squ. (CS):             0.6458
Covariance Type:            nonrobust
=====================================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
-------------------------------------------------------------------------------------
Intercept            -0.3482      0.281     -1.240      0.215      -0.899       0.202
var1[T.B]            -0.8754      0.286     -3.061      0.002      -1.436      -0.315
var1[T.C]            -0.1272      0.333     -0.382      0.703      -0.780       0.526
var1[T.D]            -1.1880      0.306     -3.883      0.000      -1.788      -0.588
var3[T.N]            -0.0931      0.302     -0.309      0.758      -0.684       0.498
var2[T.Y]:var3[M]     0.1790      0.304      0.589      0.556      -0.416       0.774
var2[T.Y]:var3[N]     0.1332      0.304      0.439      0.661      -0.462       0.728
=====================================================================================

补充：SAS实现（用于复现R结果） #

代码 #

/* 创建包含3个变量、偏移量和案例权重的样本数据集 */
data dt;
input target var1 $ var2 $ var3 $ offset case_weights;
datalines;
1.5 A X M 1 1
2.0 A Y N 2 1.5
3.0 B X N 3 2
4.0 B Y M 2 1.5
2.5 C X N 1 1
3.5 C Y M 2 1.5
4.5 D X M 3 2
3.0 D Y N 2 1.5
4.5 A X N 1 1
5.5 A Y M 2 1.5
4.0 B X M 3 2
3.5 B Y N 2 1.5
2.5 C X M 1 1
3.0 C Y N 2 1.5
2.0 D X N 3 2
1.5 D Y M 2 1.5
2.5 A X M 1 1
3.5 A Y N 2 1.5
4.5 B X N 3 2
3.0 B Y M 2 1.5
;
run;

/* 拟合包含偏移量、案例权重和分类变量交互项的GLM模型 */
proc genmod data=dt;
class var1(REF="A") var2 var3;
model target = var1 var2*var3 / dist=gamma link=log offset=offset;
weight case_weights;
run;

/* 拟合包含偏移量、案例权重和分类变量交互项的GLM模型 */
proc hpgenselect data=dt;
class var1(REF="A") var2 var3 /param=GLM;
model target = var1 var2*var3 / Distribution=Gamma offset=offset link=log;
selection method=NONE details=all;
run;

说明 #

SAS的proc genmod和proc hpgenselect输出结果略有差异，但这是正常情况，我明白不能要求数值完全精确。

我非常希望能得到一个Python解决方案，让它复现R（以及SAS）的结果，感谢各位的帮助！

备注：内容来源于stack exchange，提问作者James Black