手动计算双因素ANOVA平方和与statsmodels结果不符的问题排查

我有一个包含1个连续变量和2个分类变量的数据集，手动从零计算双因素ANOVA的平方和与均方后，用Python的statsmodels包验证时，发现结果存在无法解释的差异。

手动计算平方和的代码 #

import pandas as pd

p, q = len(pd.unique(df.category1)), len(pd.unique(df.category2))
n1, n2 = df.groupby('category1').size(), df.groupby('category2').size()

y_groupmean1, y_groupmean2 = df.groupby('category1')['feature'].mean(), df.groupby('category2')['feature'].mean()
y_totalmean = df['feature'].mean()
y_group12 = df.groupby(['category1','category2'])['feature']

ss_b1, ss_b2 = sum(n1*(y_groupmean1 - y_totalmean)**2), sum(n2*(y_groupmean2 - y_totalmean)**2)
ss_e = sum([sum((df['feature'][(df['category1']==f'G{i}') & (df['category2']==f'H{j}')] - y_group12.get_group((f'G{i}',f'H{j}')).mean())**2) for i in range(p) for j in range(q)])

ss_t = sum((df['feature']-y_totalmean)**2)

ss_i = ss_t-ss_b1-ss_b2-ss_e

手动计算得到的ANOVA表中，各平方和有对应结果。

statsmodels验证代码 #

import statsmodels.api as sm
from statsmodels.formula.api import ols

model = ols('feature ~ C(category1) + C(category2) + C(category1):C(category2)', data=df).fit()
sm.stats.anova_lm(model, type=2)

差异现象 #

对比两者结果：

• 第一个分类变量的平方和、残差平方和完全一致
• 第二个分类变量的平方和存在明显差异，进而导致交互项平方和不同

我计算两个因素平方和的公式逻辑基本一致，却找不到差异原因。甚至参考其他双因素ANOVA的Python代码时，也发现statsmodels返回的残差和交互项结果有差异，对此十分困惑，希望得到解答。

内容的提问来源于stack exchange，提问作者Hye Jin Lee