Python statsmodel 双向方差分析给出的平方和不同于从头开始的手动计算

Question

我有一个包含一个连续变量和两个分类变量的数据集，通过使用它，我通过从头开始计算平方和和均方差来运行双向方差分析。然后我通过将结果与 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

这是使用我的代码的方差分析表： （很抱歉，您需要单击链接才能查看图片，我还没有足够的声誉将图片嵌入正文中。:( )

然后，我使用以下代码检查了 Python 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)

这是统计模型返回的方差分析表： （再次请点击链接查看图片）

如您所见，虽然 statsmodel 对第一个类别和残差平方和产生与我相同的结果，但第二个分类变量的平方和之间存在小而明显的差异，因此导致不同的平方和互动。

我认为与我的计算相比，statsmodel 仅针对两个因素中的一个给出不同的结果，这真的很奇怪，因为我使用基本相同的公式来计算两个因素的平方和。谁能帮我弄清楚为什么会出现这些差异？

谢谢！

我还检查了这个 Kaggle 代码，似乎 statsmodels 为残差项和交互项返回了不同的结果。我很困惑......请帮我解决这个问题！

Answer 1

想通了。

问题来自样本不平衡的事实。在这种情况下，ols 模型中因子变量的顺序很重要。也就是说，Y ~ a+b 与 Y ~ b+a 不同。

我做了额外的测试，调换了ols估计中的因子变量的顺序，结果表明只有最先出现在估计方程中的因子和残差在平方和方面与人工计算的结果相符。