sklearn构建类statsmodels的ANOVA表及StandardScaler使用问题

你代码中调用sm.stats.anova_lm(cancer, typ=1)输出的是Type I（序贯）平方和，iv1到iv5对应的平方和完全按照你在回归公式中书写的变量顺序，逐次计算新增变量带来的残差平方和减量，具体计算规则如下：

• 第一步：拟合仅包含截距项的空模型，此时模型预测值全为因变量dv的均值，计算空模型残差平方和SSE0 = np.sum((y_train - np.mean(y_train))**2)，也就是总平方和SST
• 第二步：仅加入第一个自变量iv1，拟合dv ~ iv1的一元线性回归模型，得到该模型的残差平方和SSE1，iv1对应的平方和SS_iv1 = SSE0 - SSE1，即加入iv1后模型残差减少的量
• 第三步：在已有截距、iv1的基础上加入第二个自变量iv2，拟合dv ~ iv1 + iv2的二元回归模型，得到残差平方和SSE2，iv2对应的平方和SS_iv2 = SSE1 - SSE2
• 第四步：按照iv3→iv4→iv5的顺序重复上述操作，每加入一个新变量，该变量对应的平方和 = 加入前模型的残差平方和 - 加入该变量后新模型的残差平方和
• 所有变量都加入后得到的全模型残差平方和，就是你已经掌握计算方法的残差项平方和，和你输出结果中的0.175166完全对应。按照这个规则逐次拟合嵌套模型计算RSS差值，就能完全复现statsmodels输出的Type I ANOVA表结果。

注意：Type I平方和的计算结果高度依赖变量进入模型的顺序，调整公式中自变量的书写顺序，各变量对应的平方和数值会发生变化。

• 线性回归场景下使用StandardScaler()做变量标准化完全合理，不存在StandardScaler仅适用于分类任务的说法。
• 线性回归中使用标准化的常见场景：
  • 需要对比不同自变量对因变量的相对影响大小时，标准化消除了变量量纲差异，回归系数的绝对值可以直接比较
  • 使用带正则项的线性模型（Ridge、Lasso、ElasticNet）时必须做标准化，否则正则项对不同量纲变量的惩罚强度不一致，结果完全失效
  • 使用梯度下降类算法求解线性回归时，标准化能大幅提升算法收敛速度、避免收敛不稳定
• 对于无正则的普通最小二乘线性回归，标准化不会改变模型的预测结果、ANOVA表的F值、p值、R²等核心统计结论，仅会改变回归系数的数值尺度，模型整体解释力完全不变。如果需要输出原始量纲下的系数做业务解释，可以选择不做标准化，或是在标准化训练完成后将系数转换回原始尺度即可。

import numpy as np
import pandas as pd
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.pipeline import Pipeline
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats

dv = [0.75, 1.00, 1.00, 0.75, 0.50, 0.75, 1.00, 1.00, 0.75, 0.50]
iv1 = [0.75, 1.00, 1.00, 0.75, 0.75, 1.00, 0.50, 0.50, 0.75, 0.25]
iv2 = [0.882, 0.867, 0.900, 0.333, 0.875, 0.500, 0.882, 0.875, 0.778, 0.867]
iv3 = [1.000, 0.067, 1.000, 0.933, 0.875, 0.500, 0.588, 0.875, 1.000, 0.467]
iv4 = [0.889, 1.000, 0.905, 0.938, 0.833, 0.882, 0.444, 0.588, 0.895, 0.812]
iv5 = [ 18, 16, 21, 16, 18, 17, 18, 17, 19, 16]
d = {'iv1': iv1, 'iv2': iv2,'iv3': iv3, 'iv4': iv4,'iv5': iv5}
d2 = {'dv':dv,'iv1': iv1, 'iv2': iv2,'iv3': iv3, 'iv4': iv4,'iv5': iv5}
df = pd.DataFrame(data=d)

x_train=df.to_numpy()
y_train=np.array(dv)

linear_model = LinearRegression(fit_intercept=True)
#pipe=Pipeline([  ('StandardScaler',StandardScaler()), ('regression',linear_model)])
pipe=Pipeline([  ('regression',linear_model)])
pipe.fit(x_train, y_train)

LR= pipe.named_steps['regression']
beta=LR.coef_
intercept=LR.intercept_
print('beta',beta)
print('intercept',intercept)

yhat=pipe.predict(x_train)
residual_ss=np.sum((y_train-yhat)**2) 

cancer=smf.ols('dv ~ iv1 + iv2 + iv3 + iv4 + iv5', data=d2).fit()
anova_table = sm.stats.anova_lm(cancer, typ=1)
print('ANOVA',anova_table) 

代码运行输出的ANOVA表结果如下：

ANOVA            df    sum_sq   mean_sq         F    PR(>F)
iv1       1.0  0.033989  0.033989  0.776151  0.428075
iv2       1.0  0.022435  0.022435  0.512307  0.513703
iv3       1.0  0.003048  0.003048  0.069596  0.804959
iv4       1.0  0.115143  0.115143  2.629357  0.180222
iv5       1.0  0.000220  0.000220  0.005023  0.946900
Residual  4.0  0.175166  0.043791       NaN       NaN

内容的提问来源于stack exchange，提问作者ankit agrawal