R:aov和lm方差分析的区别
在R中经常会用aov()和lm()两个函数进行方差分析，aov 函数的内核使用了lm算法，但二者有一定的区别。
aov() 默认（summary） 结果是基于Type I 平方和，而 lm() 默认（summary）的结果是Type III平方和。aov()分析的结果受自变量输入顺序的影响，而lm()与自变量输入顺序无关。当然这种差异是针对非平衡数据而言。对于平衡全处理的数据结构，二者分析的结果是一致的。
Ⅲ型平方和与Ⅰ型平方和
如果是等组设计，这几种平方和没有任何区别。不同类型的平方和是针对非等组设计提出的。需要通过线性模型来理解。
I型平方和也叫顺序平方和，需要考虑各主效应的进入顺序。但I型平方和不正交，顺序不同计算出的结果也不同，不推荐使用。
II型平方和叫偏序平方和，不考虑主效应顺序，但在计算主效应时不考虑交互作用。
如果各效应之间不存在交互作用，II型平方和效力较高III型平方和也叫正交平方和，不考虑顺序，但不能将总平方和完全分解。目前非等组设计中推荐使用的是III型平方和
from statsmodels.formula.api import ols #拟合线性模型的包
from statsmodels.stats.anova import anova_lm#进行方差分析的包
import pandas as pd
import numpy as np
data={
'XA1': [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,-1,-1,-1,-1,-1,-1,-1],
'XA2': [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1],
'XB1': [1,1,1,-1,-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,-1,1,1,1,-1,-1,-1,-1],
'XAB11':[1,1,1,-1,-1,-1,-1,-1, 0,0,0,0, 0,0,0,0,0,-1,-1,-1,1,1,1,1],
'XAB21':[0,0,0, 0,0,0,0,0,1,1,1,1, -1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1],
'Y':[6,3,3,6,3,5,5,2,7,6,7,5,8,7,8,9,8,9,7,5,10,12,11,7],
'A':[1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3],
'B':[1,1,1,2,2,2,2,2,1,1,1,1,2,2,2,2,2,1,1,1,2,2,2,2]
}#其中前5行是虚拟变量，Y是观测值，A,B是分类变量
df=pd.DataFrame(data)
model_none=ols('Y~1',data=df).fit()#零模型，即没有任何效应
#-----------------首先进场的三种情况-----------------
model_A=ols('Y~XA1+XA2',data=df).fit()#只有A因素，即A先入场
model_B=ols('Y~XB1',data=df).fit()#只有B因素，即B先入场
#------------------其次进场的情况----------------------
model_AB=ols('Y~XA1+XA2+XB1',data=df).fit()#A，B主效应同时存在
#------------------最后进场，全模型-------------------------
model_full=ols('Y~XA1+XA2+XB1+XAB11+XAB21',data=df).fit()#全模型
#------------------以下为方差分析结果---------------------
model_anova_ABR=ols('Y~C(A,Sum)+C(B,Sum)+C(A,Sum):C(B,Sum)',data=df).fit()
#用分类变量进行方差分析，进场顺序为ABR
anova_Result_ABR=anova_lm(model_anova_ABR,typ=1)#此时选择Ⅰ型平方和
model_anova_BAR=ols('Y~C(B,Sum)+C(A,Sum)+C(A,Sum):C(B,Sum)',data=df).fit()
#用分类变量进行方差分析，进场顺序为BAR
anova_Result_BAR=anova_lm(model_anova_BAR,typ=1)#此时选择Ⅰ型平方和先来计算进场顺序依次为A主效应，B主效应，AB交互作用的Ⅰ型平方和。看看结果如何
#以下计算进场顺序A-B-R的Ⅰ型平方和
SSEN=model_none.mse_resid*model_none.df_resid#零模型残差
SSEA=model_A.mse_resid*model_A.df_resid#A先进场模型残差
SSEAB=model_AB.mse_resid*model_AB.df_resid#AB模型残差,此时B第二个进场
SSEF=model_full.mse_resid*model_full.df_resid#全模型残差
SSA=SSEN-SSEA#A主效应
SSB=SSEA-SSEAB#B主效应
SSR=SSEAB-SSEF#交互作用效应
print('用回归计算的Ⅰ型平方和输出结果为：'
'\nSSA:%.3f,\nSSB:%.3f,\nSSR:%.3f'%(SSA,SSB,SSR))
#依次打印A主效应平方和，B主效应平方和，交互作用平方和
用回归计算的Ⅰ型平方和输出结果为：
SSA:83.766,
SSB:15.226,
SSR:7.083
print('Ⅰ型平方和方差分析输出结果为：\n',anova_Result_ABR)#打印方差分析结果
Ⅰ型平方和方差分析输出结果为：
df sum_sq mean_sq F PR(>F)
C(A, Sum) 2.0 83.765873 41.882937 17.310973 0.000064
C(B, Sum) 1.0 15.226146 15.226146 6.293241 0.021911
C(A, Sum):C(B, Sum) 2.0 7.082981 3.541490 1.463762 0.257627
Residual 18.0 43.550000 2.419444 NaN NaN
print('各主效应与交互作用平方和之和为：%.3f'%(np.sum(anova_Result_ABR.sum_sq[0:3])))
各主效应与交互作用平方和之和为：106.075可见，此时用回归模型计算的结果与用方差分析计算的结果相同，注意此时的方差分析使用的是Ⅰ型平方和。而且Ⅰ型平方和的主效应与交互作用之和刚好等于组间平方和106.075
接下来再计算进场顺序为B主效应，A主效应，AB交互作用的Ⅰ型平方和看看
#以下计算进场顺序B-A-R的Ⅰ型平方和
SSEN=model_none.mse_resid*model_none.df_resid#零模型残差
SSEB=model_B.mse_resid*model_B.df_resid#B先进场模型残差
SSEAB=model_AB.mse_resid*model_AB.df_resid#AB模型残差,此时A第二个进场
SSEF=model_full.mse_resid*model_full.df_resid#全模型残差
SSA=SSEN-SSEB#A主效应
SSB=SSEB-SSEAB#B主效应
SSR=SSEAB-SSEF#交互作用效应
print('用回归计算的Ⅰ型平方和输出结果为：'
'\nSSA:%.3f,\nSSB:%.3f,\nSSR:%.3f'%(SSA,SSB,SSR))
#依次打印A主效应平方和，B主效应平方和，交互作用平方和
用回归计算的Ⅰ型平方和输出结果为：
SSA:11.668,
SSB:87.324,
SSR:7.083
print('Ⅰ型平方和方差分析输出结果为：\n',anova_Result_BAR)#打印方差分析结果
Ⅰ型平方和方差分析输出结果为：
df sum_sq mean_sq F PR(>F)
C(B, Sum) 1.0 11.667857 11.667857 4.822536 0.041435
C(A, Sum) 2.0 87.324162 43.662081 18.046325 0.000050
C(A, Sum):C(B, Sum) 2.0 7.082981 3.541490 1.463762 0.257627
Residual 18.0 43.550000 2.419444 NaN NaN
print('各主效应与交互作用平方和之和为：%.3f'%(np.sum(anova_Result_BAR.sum_sq[0:3])))
各主效应与交互作用平方和之和为：106.075
详细阅读：用回归理解方差分析
为了更好理解这篇文章，你可能需要了解：
- 两因素方差分析
- 平方和的分解
- 方差分析模型
- 虚拟变量
双因素方差分析python
from statsmodels.formula.api import ols #拟合线性模型的包
from statsmodels.stats.anova import anova_lm#进行方差分析的包
import pandas as pd
import numpy as np
data={
'XA1': [1,1,1,1, 1, 1, 1, 1,0,0,0,0,0,0,0,0,-1,-1,-1,-1,-1,-1,-1,-1],
'XA2': [0,0,0,0, 0, 0, 0, 0,1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1,-1],
'XB1': [1,1,1,1,-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1],
'XAB11':[1,1,1,1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,-1,-1,-1,-1,1,1,1,1],
'XAB21':[0,0,0,0,0,0,0,0,1,1,1,1,-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1],
'Y':[5,3,4,2,6,3,5,5,7,6,7,8,8,7,8,5,9,7,5,7,10,12,11,7],
'A':[1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3],
'B':[1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2,1,1,1,1,2,2,2,2]
}#其中前5行是虚拟变量，Y是观测值，A,B是分类变量
df=pd.DataFrame(data)
model_lin=ols('Y~XA1+XA2+XB1+XAB11+XAB21',data=df).fit()#用虚拟变量进行回归
model_anova=ols('Y~C(A)*C(B)',data=df).fit()#用分类变量进行方差分析
anova_Result=anova_lm(model_anova)
print(anova_Result)
df sum_sq mean_sq F PR(>F)
C(A) 2.0 79.083333 39.541667 17.905660 0.000052
C(B) 1.0 12.041667 12.041667 5.452830 0.031315
C(A):C(B) 2.0 9.083333 4.541667 2.056604 0.156887
Residual 18.0 39.750000 2.208333 NaN NaN
#利用方差分析结果计算决定系数，组间平方和（主效应与交互作用之和）/总平方和
R_square=np.sum(anova_Result.sum_sq[0:3])/np.sum(anova_Result.sum_sq)
print(R_square)
0.7159869008633524
#打印回归分析结果
print(model_lin.summary())
OLS Regression Results
==============================================================================
Dep. Variable: Y R-squared: 0.716
Model: OLS Adj. R-squared: 0.637
Method: Least Squares F-statistic: 9.075
Date: Sat, 11 Apr 2020 Prob (F-statistic): 0.000191
Time: 18:02:59 Log-Likelihood: -40.109
No. Observations: 24 AIC: 92.22
Df Residuals: 18 BIC: 99.29
Df Model: 5
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
Intercept 6.5417 0.303 21.566 0.000 5.904 7.179
XA1 -2.4167 0.429 -5.633 0.000 -3.318 -1.515
XA2 0.4583 0.429 1.068 0.299 -0.443 1.360
XB1 -0.7083 0.303 -2.335 0.031 -1.346 -0.071
XAB11 0.0833 0.429 0.194 0.848 -0.818 0.985
XAB21 0.7083 0.429 1.651 0.116 -0.193 1.610
==============================================================================
Omnibus: 1.556 Durbin-Watson: 2.330
Prob(Omnibus): 0.459 Jarque-Bera (JB): 1.295
Skew: -0.535 Prob(JB): 0.523
Kurtosis: 2.614 Cond. No. 1.73
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
#利用回归模型输出各预测值
model_lin.predict()
array([ 3.5 , 3.5 , 3.5 , 3.5 , 4.75, 4.75, 4.75, 4.75, 7. ,
7. , 7. , 7. , 7. , 7. , 7. , 7. , 7. , 7. ,
7. , 7. , 10. , 10. , 10. , 10. ])
#利用回归模型输出参数估计值
model_lin.params
Intercept 6.541667
XA1 -2.416667
XA2 0.458333
XB1 -0.708333
XAB11 0.083333
XAB21 0.708333
dtype: float64
#输出原始观测值的总均值
np.mean(df['Y'])
6.5416666666666672 不等组设计的平方和Ⅲ型
from statsmodels.formula.api import ols #拟合线性模型的包
from statsmodels.stats.anova import anova_lm#进行方差分析的包
import pandas as pd
import numpy as np
data={
'XA1': [1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,-1,-1,-1,-1,-1,-1,-1],
'XA2': [0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,-1,-1,-1,-1,-1,-1,-1],
'XB1': [1,1,1,-1,-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,-1,1,1,1,-1,-1,-1,-1],
'XAB11':[1,1,1,-1,-1,-1,-1,-1, 0,0,0,0, 0,0,0,0,0,-1,-1,-1,1,1,1,1],
'XAB21':[0,0,0, 0,0,0,0,0,1,1,1,1, -1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1],
'Y':[6,3,3,6,3,5,5,2,7,6,7,5, 8,7,8,9,8, 9,7,5,10,12,11,7],
'A':[1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3],
'B':[1,1,1,2,2,2,2,2,1,1,1,1,2,2,2,2,2,1,1,1,2,2,2,2]
}#其中前5行是虚拟变量，Y是观测值，A,B是分类变量
df=pd.DataFrame(data)
model_full=ols('Y~XA1+XA2+XB1+XAB11+XAB21',data=df).fit()#全模型回归
model_A=ols('Y~XB1+XAB11+XAB21',data=df).fit()#去掉A因素主效应
model_B=ols('Y~XA1+XA2+XAB11+XAB21',data=df).fit()#去掉B因素主效应
model_R=ols('Y~XA1+XA2+XB1',data=df).fit()#去掉交互作用
model_anova=ols('Y~C(A,Sum)*C(B,Sum)',data=df).fit()#用分类变量进行方差分析
anova_Result=anova_lm(model_anova,typ=3)
SSE=model_full.mse_resid*model_full.df_resid#全模型残差
SSEA=model_A.mse_resid*model_A.df_resid#去A模型残差
SSEB=model_B.mse_resid*model_B.df_resid#去B模型残差
SSER=model_R.mse_resid*model_R.df_resid#去交互模型残差
print(SSEA-SSE,SSEB-SSE,SSER-SSE)
#依次打印A主效应平方和，B主效应平方和，交互作用平方和
74.03142710822807 15.639893617021286 7.082980539433265
print(anova_Result)#打印方差分析结果
sum_sq df F PR(>F)
Intercept 993.384574 1.0 410.583751 7.687336e-14
C(A, Sum) 74.031427 2.0 15.299262 1.311731e-04
C(B, Sum) 15.639894 1.0 6.464250 2.041777e-02
C(A, Sum):C(B, Sum) 7.082981 2.0 1.463762 2.576274e-01
Residual 43.550000 18.0 NaN NaN
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
