多因素方差分析中预测因素的筛多_用回归来理解方差分析（二）：两因素方差分析...

1 两因素方差分析的形式

多因素方差分析针对的是多因素完全随机设计。包含两个及以上的自变量，为便于讲解，本文以两因素方差分析为例。

在一个两因素完全随机设计中，自变量

包含

，共

个水平。自变量

包含

，共

个水平。总共形成

个处理。各处理

之下有

个观测值。特别的，当各处理的观测值都为

时，称为等组设计。

等组设计的具体形式如下：

若各组观测值

数量不同，称为非等组设计，形式如下：

其中

各不相同

两因素方差分析的逻辑是对变异进行比较。通过对平方和(Sum of Squares)和自由度(Degree of Freedom)的分解，利用公式：

计算出各主效应方差，交互作用方差和组内方差（即误差方差）。各主效应和交互作用的方差与误差方差的差异（比值）体现了该效应的显著性。而且在零假设为真的情况下，该比值服从F分布，进而分别完成各效应的假设检验。

具体的分解如下图：

常规的流程大家都非常熟悉，不再赘述。接下来换个角度，从多元线性回归来看一看两因素方差分析。

2 两因素方差分析的线性模型

模型的基本形式

两因素方差分析的线性模型可以表示为：

其中

指A因素第

个水平的效应，

指

因素第

个水平的效应，

指A因素第

个水平与B因素第

个水平的交互作用效应，

表示

处理中，第

个观测值的随机误差。该式说明：某个观测值等于其所在处理的各主效应和交互作用效果叠加，并附加上一个随机误差组成。

通常，我们让:

其中

表示各平均效应。上述式子意味着：以各因素的平均效应为基准，用相对的效应值

来表示各水平的效应。可见，就是

是相对效应量，是各效应偏离平均效应的程度，即离均差。易知：

2.2.4和2.2.5式需要稍微推导一下，有兴趣的同学可以自己动手试试，没兴趣忽略就好

据此，可将2.1式写为：

其中:

是观测值
是总平均，相当于是截距
是随机误差，且独立同分布

模型表达的意思是：每一个观测值

实际上都是在某个总平均

上附加了一个

因素

水平的效果

，一个

因素

水平的效果

，

因素和

因素在各自水平上的交互作用

，并且还受到了随机误差

的影响。

例如，如果某被试接受的是自变量

的处理，那么他的观测值就是

是指该处理中观测值的编号，

引入虚拟变量

为了便于理解，我们可以引入虚拟变量

，令：

意味着：在

处理中，

,其余虚拟变量均为0。

例如：对第

处理中，

,其余均为0

可以发现，当把虚拟变量带入2.3式之后得到：

2.4式可看成是一个包含

个自变量的多元线性回归模型。

此时

的取值为：

此时

的取值为：

是不是眼花缭乱？哈哈。

如果从矩阵的角度来看，此时的

形成的矩阵就是一个

的单位矩阵，只有对角线元素为1。而且注意到：

的取值都可以由

组合得到。即

的每一列都可以是

列向量的线性组合。也就是自变量之间存在多重共线性，此时有效的自变量数量仅仅只有

的

个

。这就出问题了，因为我们需要估计的参数一共有

个

别担心，我们还有2.2.2-2.2.5这几个式子，利用这几个式子将上述取值改为：（又要眼花缭乱了）

此时利用2.2.1和2.2.2式将

和

中各去掉了一个自变量，表中去掉的是各自的最后一项，因为

，所以当

时，将

都取为

就可以了。实际上去掉任何一项都可以，不影响最终结果。

同理，可将交互作用的虚拟变量改为：

注意到对角线的子矩阵，以及最后一行的各矩阵。特别是最后一行

时的取值，用心的小伙伴简单推导就可以得到，这里就不再啰嗦了。

3 对单因素方差分析的线性模型进行多元回归

提示：该部分直接用公式推导实在太麻烦了，所以用一个数据例子来进行说明。

以一个3×2的两因素完全随机设计为例，其中A因素个水平，B因素两个水平，每个处理之下一共4个观测值。

将观测值

和虚拟变量都列出来为。为了偷懒，除了第一个处理的四个观测值都列出来之外，其余的每个处理只列出一个。

接下来，利用python的statemodels进行计算：

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.paramsIntercept    6.541667
XA1         -2.416667
XA2          0.458333
XB1         -0.708333
XAB11        0.083333
XAB21        0.708333
dtype: float64
#输出原始观测值的总均值
np.mean(df['Y'])6.541666666666667

可见，用常规的方差分析计算出来的结果与回归分析一致，决定系数

是相同的。此外，回归模型告诉我们：各处理的预测值就是该处理的组均值；模型的截距（总平均）就是所有观测值的平均值。实际上，决定系数就是回归平方和除以总平方和。而回归平方和等于观测值与预测值的平方和，对应的就是方差分析组间平方和，不过多因素方差分析中，组间平方和等于所有主效应与交互作用平方和之和。回归总平方和就是观测值与总平均的平方和，此处与方差分析的总平方和一致。

绕了这么半天，似乎说了好多废话。有同学可能会嫌弃，认为我把简单的事情搞复杂了。确实这样理解很复杂，但是只有这样理解了，我们才能搞清楚什么是Ⅰ型平方和以及Ⅲ型平方和。这两种平方和专门针对多因素方差分析，在等组设计中二者没有区别，但是非等组设计中二者区别就大了，反正我最终是通过回归才搞清楚。

关于Ⅰ型平方和以及Ⅲ型平方和，我会专门写一篇文章来介绍。

4 结论

多元回归的结果与方差分析的结果是一致的
将多因素方差分析转换为回归分析时，要注意虚拟变量的取值，特别是交互作用的虚拟变量。
通过回归来理解方差分析，有助于我们理解Ⅰ型平方和以及Ⅲ型平方和（另开文章介绍）

欢迎非商业转载，只需注明作者“AhaDad”和来源即可