statsmodels为何包含参考组交互项？如何实现SAS式无参考组交互回归？

问题背景 #

使用Statsmodels构建包含三个4水平分类变量的Tweedie回归模型时，发现Statsmodels会输出包含参考组的交互项（如V2=W与V3=K的交互项），但此前用SAS做同类回归时，SAS会自动排除这类参考组交互项。

Statsmodels示例代码 #

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

data = pd.DataFrame({
    'V1': pd.Categorical(['A', 'B', 'C', 'D', 'A', 'B', 'C', 'D']),
    'V2': pd.Categorical(['W', 'X', 'Y', 'Z', 'W', 'X', 'Y', 'Z']),
    'V3': pd.Categorical(['K', 'L', 'M', 'N', 'K', 'L', 'M', 'N']),
    'y': [5.1, 7.3, 6.9, 8.0, 5.4, 7.1, 6.8, 8.2]
})

formula = 'y ~ C(V1, Treatment(''A'')) + C(V2, Treatment(''W'')):C(V3, Treatment(''K''))'
model = smf.GLM.from_formula(formula, data, family=sm.families.Tweedie())
result = model.fit()

print(result.summary())

SAS示例代码 #

data example;
    input V1 $ V2 $ V3 $ y;
    datalines;
A W K 5.1
B X L 7.3
C Y M 6.9
D Z N 8.0
A W K 5.4
B X L 7.1
C Y M 6.8
D Z N 8.2
;
run;

proc hpgenselect data=example;
    class V1 (ref='A') V2 (ref='W') V3 (ref='K'); 
    model y = V1 V2*V3 / dist=tweedie link=log; 
run;

差异原因 #

• SAS的默认行为：当通过class语句指定分类变量及参考组后，SAS在生成交互项时会自动排除两个变量参考组的交叉组合。这类项属于冗余项，其效应已被模型截距或主效应的参考组基准覆盖，保留会引发共线性，因此SAS默认移除。
• Statsmodels的默认行为：Statsmodels基于patsy库处理公式语法，当手动用C(var, Treatment(ref)):C(var2, Treatment(ref2))指定交互项时，patsy会生成所有可能的交叉组合，包括参考组的交叉项，不会自动过滤冗余项，需要用户显式控制。

在Statsmodels中实现SAS式的无参考组交互项回归 #

可以通过在公式中显式排除参考组的交互项来实现，示例代码如下：

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

data = pd.DataFrame({
    'V1': pd.Categorical(['A', 'B', 'C', 'D', 'A', 'B', 'C', 'D']),
    # 提前指定分类变量的顺序，确保参考组排在首位
    'V2': pd.Categorical(['W', 'X', 'Y', 'Z', 'W', 'X', 'Y', 'Z'], categories=['W', 'X', 'Y', 'Z']),
    'V3': pd.Categorical(['K', 'L', 'M', 'N', 'K', 'L', 'M', 'N'], categories=['K', 'L', 'M', 'N']),
    'y': [5.1, 7.3, 6.9, 8.0, 5.4, 7.1, 6.8, 8.2]
})

# 生成V2与V3的交互项，并排除参考组(W和K)的交叉项
formula = 'y ~ C(V1, Treatment("A")) + V2:V3 - V2[W]:V3[K]'
# 指定log链接，与SAS示例保持一致
model = smf.GLM.from_formula(formula, data, family=sm.families.Tweedie(link=sm.families.links.log()))
result = model.fit()

print(result.summary())

说明 #

• 提前指定分类变量的categories参数，确保参考组处于首位，方便后续精准排除目标交互项。
• 公式中V2:V3生成所有交互项，- V2[W]:V3[K]显式移除参考组的交叉组合，最终输出的交互项数量与SAS完全一致。

内容的提问来源于stack exchange，提问作者Felippe Trigueiro