2020.6.29 概率统计-Task04-方差分析

方差分析

一、概要
- 1 从独立样本t检验到方差分析
- 2 方差分析定义与分类
二、单因素方差分析
- 1 推导过程
- 2 代码实现
三、双（多因素）因素方差分析
- 1 推导过程
- 2 代码实现
最后

概率统计的最后一节内容，有关方差分析的介绍。Datawhale的课件讲F检验的推导过程讲解的十分详细。在模型实现方面，只给出了R语言的实现途径。笔者将结合另外一个例子，介绍SPSS对于方差分析的实现过程，并对方差分析种类的介绍进行扩展，后期会补充python的实现。

一、概要

1 从独立样本t检验到方差分析

首先举个例子：

如果你的数学成绩和班里的小王差了1分，那么你和小王的分数有没有很大的差异呢？
你们班和隔壁班的数学平均成绩差了1分，那么两个班的水平有没有很大的差异呢？
你们省和隔壁省，用的一套卷子，数学平均分差了1分

因此，当总体很大的情况下，我们就需要一个固定且客观的标准用以判断两个群体的数值是否有差异，而不是主观臆断。在这一背景下，t检验被发明了出来。

1 + 2：单样本t检验——一批样本和一个固定的数值，其他数据也是和1比较
2 + 3：独立样本t检验——独立的两批人，4 + 5也是独立样本t检验
2 + 4：配对样本t检验——配对的数据，3 + 5也是独立样本t检验

但是，独立样本t检验有其局限性，每次只能确定两组是否存在差异，因此，就需要引入方差分析

2 方差分析定义与分类

方差分析(Analysis of variance, ANOVA) 主要研究分类变量作为自变量时，对因变量的影响是否是显著的。

方差分析的方法是由20世纪的统计学家Ronald Aylmer Fisher在1918年到1925年之间提出并陆续完善起来的，该方法刚开始是用于解决田间实验的数据分析问题，因此，方差分析的学习是和实验设计、实验数据的分析密不可分的。

实验设计和方差分析都有自己相应的语言。因此，在这里我们通过一个焦虑症治疗的实例，先了解一些术语，并且思考一下，方差分析主要用于解决什么样的问题。

以焦虑症治疗为例，现有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10位焦虑症患者作为志愿者，随机分配一半的人接受为期五周的CBT，另外一半接受为期五周的EMDR，设计方案如表1-1所示。在治疗结束时，要求每位患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。

表1-1 单因素组件方差分析

CBT	EMDR
s1	s4
s2	s5
s3	s6

在这个实验设计中，治疗方案是两水平（CBT、EMDR）的组间因子。之所以称其为组间因子，是因为每位患者都仅被分配到一个组别中，没有患者同时接受CBT和EMDR。表中字母s代表受试者（患者）。STAI是因变量，治疗方案是自变量。由于在每种治疗方案下观测数相等，因此这种设计也称为均衡设计（balanced design）；若观测数不同，则称作非均衡设计（unbalanced design）。

因为仅有一个类别型变量，表1的统计设计又称为单因素方差分析（one-way ANOVA），或进一步称为单因素组间方差分析。方差分析主要通过F检验来进行效果评测，若治疗方案的F检验显著，则说明五周后两种疗法的STAI得分均值不同。

假设你只对CBT的效果感兴趣，则需将10个患者都放在CBT组中，然后在治疗五周和六个月后分别评价疗效，设计方案如表1-2所示。

表1-2 单因素组内方差分析

	时间
患者	5周	6个月
s1
s2
s3
s4
s5
s6

此时，时间（time）是两水平（五周、六个月）的组内因子。因为每位患者在所有水平下都进行了测量，所以这种统计设计称单因素组内方差分析；又由于每个受试者都不止一次被测量，也称作重复测量方差分析。当时间的F检验显著时，说明患者的STAI得分均值在五周和六个月间发生了改变。

现假设你对治疗方案差异和它随时间的改变都感兴趣，则将两个设计结合起来即可：随机分配五位患者到CBT，另外五位到EMDR，在五周和六个月后分别评价他们的STAI结果（见表1-3）。

表1-3 含组间和组内因子的双因素方差分析

		时间
疗法	患者	5周	6个月
CBT	s1
	s2
	s3
EMDR	s4
	s5
	s6

疗法（therapy）和时间（time）都作为因子时，我们既可分析疗法的影响（时间跨度上的平均）和时间的影响（疗法类型跨度上的平均），又可分析疗法和时间的交互影响。前两个称作主效应，交互部分称作交互效应。

当设计包含两个甚至更多的因子时，便是因素方差分析设计，比如两因子时称作双因素方差分析，三因子时称作三因素方差分析，以此类推。若因子设计包括组内和组间因子，又称作混合模型方差分析，当前的例子就是典型的双因素混合模型方差分析。

本例中，你将做三次F检验：疗法因素一次，时间因素一次，两者交互因素一次。若疗法结果显著，说明CBT和EMDR对焦虑症的治疗效果不同；若时间结果显著，说明焦虑度从五周到六个月发生了变化；若两者交互效应显著，说明两种疗法随着时间变化对焦虑症治疗影响不同（也就是说，焦虑度从五周到六个月的改变程度在两种疗法间是不同的）。

现在，我们对上面的实验设计稍微做些扩展。众所周知，抑郁症对病症治疗有影响，而且抑郁症和焦虑症常常同时出现。即使受试者被随机分配到不同的治疗方案中，在研究开始时，两组疗法中的患者抑郁水平就可能不同，任何治疗后的差异都有可能是最初的抑郁水平不同导致的，而不是由于实验的操作问题。抑郁症也可以解释因变量的组间差异，因此它常称为混淆因素（confounding factor）。由于你对抑郁症不感兴趣，它也被称作干扰变数（nuisance variable）。

假设招募患者时使用抑郁症的自我评测报告，比如白氏抑郁症量表（BDI），记录了他们的抑郁水平，那么你可以在评测疗法类型的影响前，对任何抑郁水平的组间差异进行统计性调整。本案例中，BDI为协变量，该设计为协方差分析（ANCOVA）。

以上设计只记录了单个因变量情况（STAI），为增强研究的有效性，可以对焦虑症进行其他的测量（比如家庭评分、医师评分，以及焦虑症对日常行为的影响评价）。当因变量不止一个时，设计被称作多元方差分析（MANOVA），若协变量也存在，那么就叫多元协方差分析（MANCOVA）。

在此基础上，我们通过另外一个更为直观的例子来了解这些不同的种类

①+⑤：单因素方差分析（One-Way Analysis of Variance，One-Way ANOVA）
①+⑤+⑦+⑨：单因素重复测量方差分析（One-Way repeated measures ANOVA）
①+②+⑤：两因素方差分析（Two-Way ANOVA）
还可以四因素、五因素、六因素，统称为多因素方差分析（Multi-way ANOVA）
①+②+⑤+⑦+⑨：两因素重复测量方差分析（Two-Way Repeated Measures ANOVA）
①+④+⑤：协方差分析（Analysis of Covariance，ANCOVA）
为啥①+②+⑤叫做两因素方差分析，把②换成④，就成了协方差分析呢？②是影响因素，④也是影响因素啊！解释：因为协变量往往是和自变量有着很高相关性的变量，它并不是我们试验研究的重点，我们只是把它当作控制变量来处理
协变量（covariate）也是自变量的一种
①+⑤+⑥：多元方差分析（Multivariate ANOVA，MANOVA）
多因素方差分析：Multi-way
多元方差分析：Multivariate

二、单因素方差分析

1 推导过程

接下来我们使用种小麦的例子，去帮助理解方差分析里涉及的一些变量。

假设我们现在有若干品种的小麦，要在某一地区播种，我们想知道这些品种的产量有没有显著区别，为此我们先设计了一个田间实验，取一大块地将其分成形状大小都相同的nnn小块．设供选择的品种有kkk个，我们打算其中的n1n_1n1小块种植品种1, n2n_2n2小块种植品种2，等等，n1+n2+...nk=nn_1 + n_2 + ... n_k = nn1+n2+...nk=n.

接下来，我们使用方差分析的方法去看不同小麦品种的产量是否有显著差异。

设问题中涉及一个因素AAA，有kkk个水平，如上例的kkk个种子品种，以YijY_{ij}Yij记第iii个水平的第jjj个观察值，如上例YijY_{ij}Yij是种植品种iii的第jjj小块地上的亩产量。模型为
Yij=ai+eij,j=1,...,ni,i=1,...,k(2.1)Y_{ij} = a_i + e_{ij}, j = 1,...,n_i, i = 1,...,k\qquad(2.1) Yij=ai+eij,j=1,...,ni,i=1,...,k(2.1)
aia_iai表示水平iii的理论平均值，称为水平iii的效应。在小麦例子中，aia_iai就是品种iii的平均亩产量，eije_{ij}eij就是随机误差。并且我们假定：
E(eij)=0,0<Var(eij)=σ2<∞,一切eij独立同分布(2.2)E(e_{ij})=0, 0<Var(e_{ij})={\sigma}^2<\infty,一切e_{ij}独立同分布\qquad(2.2) E(eij)=0,0<Var(eij)=σ2<∞,一切eij独立同分布(2.2)
因素AAA的各水平的高低优劣，取决于其理论平均aia_{i}ai的大小。故对模型(2.1)，我们头一个关心的事情，就是诸aia_{i}ai是否全相同。如果是，则表示因素AAA对所考察的指标YYY其实无影响．这时我们就说因素A的效应不显著，否则就说它显著。当然，在实际应用中，所谓“显著”，是指诸aia_{i}ai之间的差异要大到一定的程度．这个 “一定的程度”，是从其实用上的意义着眼，而“统计显著性”，则是与随机误差相比而言．这点在下文的讨论中会有所体现．我们把所要检验的假设写为：
H0:a1=a2=⋯=ak(2.3)H_0:a_1=a_2=\cdots=a_k \qquad (2.3) H0:a1=a2=⋯=ak(2.3)
为检验该假设，我们需要分析，为什么各个YijY_{ij}Yij会有差异？从模型(2.1)来看，无非两个原因：一是各aia_{i}ai可能有差异．例如，若 a1>a2a_1>a_2a1>a2, 这就使Y1jY_{1j}Y1j倾向于大于Y2jY_{2j}Y2j；二是随机误差的存在。这一分析启发了如下的想法：找一个衡量全部yijy_{ij}yij的变异的量：
SS=∑i=1k∑j=1ni(Yij−Yˉ)2,Yˉ=∑i=1k∑j=1niYij/n(2.4)SS= \sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y} \right )^2, \qquad \bar{Y}=\sum_{i=1}^{k}\sum_{j=1}^{n_i}Y_{ij}/n \qquad (2.4) SS=i=1∑kj=1∑ni(Yij−Yˉ)2,Yˉ=i=1∑kj=1∑niYij/n(2.4)
SSSSSS愈大，表示YijY_{ij}Yij之间的差异越大。

接下来，把SSSSSS分为两部分，一部分表示随机误差的影响，记为SSeSS_eSSe；另一部分表示因素AAA的各水平理论平均值aia_iai不同带来的影响，记为SSASS_ASSA。

关于SSeSS_eSSe，先固定一个iii，此时对应的所有观测值Yi1,Yi2,⋯,YinY_{i1},Y_{i2},\cdots,Y_{in}Yi1,Yi2,⋯,Yin，他们之间的差异与每个水平的理论平均值不等无关，而是取决于随机误差，反映这些观察值差异程度的量是∑j=1ni(Yij−Yiˉ)2\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2∑j=1ni(Yij−Yiˉ)2，其中
Yiˉ=(Yi1+Yi2+⋯+Yin)/ni,i=1,2,⋯,n(2.5)\bar{Y_i}=(Y_{i1}+Y_{i2}+\cdots+Y_{in})/n_i,\quad i=1, 2,\cdots,n \qquad (2.5) Yiˉ=(Yi1+Yi2+⋯+Yin)/ni,i=1,2,⋯,n(2.5)
Yiˉ\bar{Y_i}Yiˉ可以视为对aia_iai的估计。把上述平方和做累加得：
SSe=∑i=1k∑j=1ni(Yij−Yiˉ)2(2.6)SS_e=\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2 \qquad (2.6) SSe=i=1∑kj=1∑ni(Yij−Yiˉ)2(2.6)
可求得SSASS_ASSA:
SSA=SS−SSe=∑i=1k∑j=1ni(Yij−Yˉ)2−∑i=1k∑j=1ni(Yij−Yiˉ)2=∑i=1k∑j=1ni((Yij−Yiˉ)−(Yiˉ−Yˉ))2−∑i=1k∑j=1ni(Yij−Yiˉ)2=∑i=1k∑j=1ni((Yij−Yiˉ)2−2(Yij−Yiˉ)(Yiˉ−Yˉ)+(Yiˉ−Yˉ)2)−∑i=1k∑j=1ni(Yij−Yiˉ)2=∑i=1k∑j=1ni(Yij−Yiˉ)2−2∑i=1k∑j=1ni((Yij−Yiˉ)(Yiˉ−Yˉ))+∑i=1k∑j=1ni(Yiˉ−Yˉ)2−∑i=1k∑j=1ni(Yij−Yiˉ)2=∑i=1k∑j=1ni(Yiˉ−Yˉ)2−2∑i=1k((Yiˉ−Yˉ)∑j=1ni(Yij−Yiˉ))(ps:∑j=1ni(Yij−Yiˉ)=0)=∑i=1kni(Yiˉ−Yˉ)2(2.7)\begin{aligned} SS_A &= SS-SS_e \\ &=\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y} \right )^2-\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2 \\ &=\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( (Y_{ij}-\bar{Y_i})-(\bar{Y_i}-\bar{Y}) \right )^2-\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2 \\ &=\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( (Y_{ij}-\bar{Y_i})^2-2(Y_{ij}-\bar{Y_i})(\bar{Y_i}-\bar{Y})+(\bar{Y_i}-\bar{Y})^2 \right )-\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2 \\ &=\sum_{i=1}^{k}\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y_i})^2 - 2\sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( (Y_{ij}-\bar{Y_i})(\bar{Y_i}-\bar{Y}) \right ) +\sum_{i=1}^{k}\sum_{j=1}^{n_i}(\bar{Y_i}-\bar{Y})^2 - \sum_{i=1}^{k}\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2 \\ &= \sum_{i=1}^{k}\sum_{j=1}^{n_i}(\bar{Y_i}-\bar{Y})^2 - 2\sum_{i=1}^{k}\left ( (\bar{Y_i}-\bar{Y})\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y_i}) \right ) \quad (ps:\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y_i})=0) \\ &= \sum_{i=1}^{k}n_i(\bar{Y_i}-\bar{Y})^2 \qquad (2.7) \end{aligned} SSA=SS−SSe=i=1∑kj=1∑ni(Yij−Yˉ)2−i=1∑kj=1∑ni(Yij−Yiˉ)2=i=1∑kj=1∑ni((Yij−Yiˉ)−(Yiˉ−Yˉ))2−i=1∑kj=1∑ni(Yij−Yiˉ)2=i=1∑kj=1∑ni((Yij−Yiˉ)2−2(Yij−Yiˉ)(Yiˉ−Yˉ)+(Yiˉ−Yˉ)2)−i=1∑kj=1∑ni(Yij−Yiˉ)2=i=1∑kj=1∑ni(Yij−Yiˉ)2−2i=1∑kj=1∑ni((Yij−Yiˉ)(Yiˉ−Yˉ))+i=1∑kj=1∑ni(Yiˉ−Yˉ)2−i=1∑kj=1∑ni(Yij−Yiˉ)2=i=1∑kj=1∑ni(Yiˉ−Yˉ)2−2i=1∑k((Yiˉ−Yˉ)j=1∑ni(Yij−Yiˉ))(ps:j=1∑ni(Yij−Yiˉ)=0)=i=1∑kni(Yiˉ−Yˉ)2(2.7)
因为Yiˉ\bar{Y_i}Yiˉ可以视为对aia_iai的估计，aia_iai的差异越大，Yiˉ\bar{Y_i}Yiˉ之间的差异也越大，所以SSASS_ASSA可以用来衡量不同水平之间的差异程度。

在统计学上，通常称SSSSSS为总平方和，SSASS_ASSA为因素AAA的平方和，SSeSS_eSSe为误差平方和，分解式SS=SSA+SSeSS=SS_A+SS_eSS=SSA+SSe为该模型的方差分析。

基于上面的分析，我们可以得到假设（5.3）的一个检验方法：当比值SSA/SSeSS_A/SS_eSSA/SSe大于某一给定界限时，否定H0H_0H0，不然就接受H0H_0H0。为了构造FFF分布的检验统计量，我们假定随机误差eije_{ij}eij满足正态分布N(0,σ2)N(0, \sigma^2)N(0,σ2)，同时我们也假定观察值YijY_{ij}Yij符合正态分布，此时，记
MSA=SSA/(k−1),MSe=SSe/(n−k)(2.8)MS_A = SS_A/(k-1), \quad MS_e = SS_e/(n-k) \qquad (2.8) MSA=SSA/(k−1),MSe=SSe/(n−k)(2.8)
当H0H_0H0成立时，有：
MSA/MSe∼Fk−1,n−k(2.9)MS_A / MS_e \sim F_{k-1, n-k} \qquad (2.9) MSA/MSe∼Fk−1,n−k(2.9)
据（5.9），在给定显著性水平α\alphaα时，即得（5.3）的假设H0H_0H0的检验如下：
当MSA/MSe⩽Fk−1,n−k(α)时，接受H0，不然就拒绝H0(2.10)当MS_A / MS_e \leqslant F_{k-1, n-k}(\alpha)时，接受H_0，不然就拒绝H_0 \qquad (2.10) 当MSA/MSe⩽Fk−1,n−k(α)时，接受H0，不然就拒绝H0(2.10)
MSAMS_AMSA和MSeMS_eMSe分别被称为因素AAA和随机误差的平均平方和。被除数k−1k-1k−1和n−kn-kn−k，分别称为这两个平方和的自由度。MSeMS_eMSe的自由度为什么是n−kn-kn−k呢？因为平方和∑j=1ni(Yij−Yiˉ)2\sum_{j=1}^{n_i}\left ( Y_{ij}-\bar{Y_i} \right )^2∑j=1ni(Yij−Yiˉ)2的自由度为ni−1n_i-1ni−1，故对iii求和，SSeSS_eSSe的自由度就是n−kn-kn−k。那么，MSAMS_AMSA的自由度为什么是k−1k-1k−1呢？因为一共有kkk个平均值a1,⋯,aka_1,\cdots,a_ka1,⋯,ak等k−1k-1k−1个，故自由度为k−1k-1k−1，两者自由度之和为n−1n-1n−1，恰好是总平方和的自由度。

到这里，我们可以做出方差分析表如表2-1

2-1 单因素方差分析的方差分析表

项目	SSSSSS	自由度	MSMSMS	FFF比	显著性
AAA	SSASS_ASSA	k−1k-1k−1	MSAMS_AMSA	MSA/MSeMS_A / MS_eMSA/MSe	, *, 或无
误差	SSeSS_eSSe	n−kn-kn−k	MSeMS_eMSe
总和	SSSSSS	n−1n-1n−1

在上表中，对于显著性一栏，一般来说，我们把算出的FFF比，即MSA/MSeMS_A / MS_eMSA/MSe，与Fk−1,n−k(0.05)=c1F_{k-1, n-k}(0.05)=c_1Fk−1,n−k(0.05)=c1和Fk−1,n−k(0.01)=c2F_{k-1, n-k}(0.01)=c_2Fk−1,n−k(0.01)=c2比较。若MSA/MSe<c2MS_A / MS_e<c_2MSA/MSe<c2，用**表示，表明A因素的效应是高度显著的，即在α=0.01\alpha=0.01α=0.01的显著性水平下，拒绝原假设（5.3）。同理，c2<MSA/MSe<c1c_2<MS_A / MS_e<c_1c2<MSA/MSe<c1用∗\ast∗表示，MSA/MSe>c1MS_A / MS_e>c_1MSA/MSe>c1时不显著。

2 代码实现

由于目前还没有接触过R，因此这部分的坑之后再来填，给出spss的实现过程和结果

三、双（多因素）因素方差分析

1 推导过程

在很多种情况下，只考虑一个指标对观察值的影响，显然是不够的，这时就会用到多因素方差分析。双因素方差分析和多因素方差分析在原理上是相似的，这里为了书写简便，我们只以双因素方差分析为例进行推导。

还是以田间实验的例子帮助理解推导过程，我们设有两个因素A,BA, BA,B，分别有k,lk, lk,l个水平（例如AAA为品种，有kkk个；BBB为播种量，考虑lll种不同的数值，如20斤／亩，25斤／亩，……）．AAA的水平iii与BBB的水平jjj的组合记为(i,j)(i,j)(i,j)，其试验结果记为 Yij,i=1,⋅⋅⋅,k,j=1,…,lY_{ij}, i = 1, · · ·, k,j = 1,…, lYij,i=1,⋅⋅⋅,k,j=1,…,l．统计模型定为
Yij=μ+ai+bj+eij，i=1,⋅⋅⋅,k,j=1,⋅⋅⋅,l(3.1)Y_{ij} = \mu + a_i + b_j + e_{ij}，i= 1, · · ·, k,j = 1,· · ·, l\qquad (3.1) Yij=μ+ai+bj+eij，i=1,⋅⋅⋅,k,j=1,⋅⋅⋅,l(3.1)
为解释这模型，首先把右边分成两部分：eije_{ij}eij为随机误差，它包含了未加控制的因素(A,BA,BA,B以外的因素）及大量随机因素的影响．假定
E(eij)=0,0<Var(eij)=σ2<∞,一切eij独立同分布(3.2)E(e_{ij})=0, 0<Var(e_{ij})={\sigma}^2<\infty,一切e_{ij}独立同分布\qquad(3.2) E(eij)=0,0<Var(eij)=σ2<∞,一切eij独立同分布(3.2)
另一部分μ+ai+bj\mu + a_i + b_jμ+ai+bj，它显示水平组合(i,j)(i,j)(i,j)的平均效应．它可以又分解为三部分：μ\muμ是总平均（一切水平组合效应的平均），是一个基准．aia_iai表示由AAA的水平iii带来的增加部分，称为因素AAA的水平iii的效应．bjb_jbj有类似的解释．调整μ\muμ的值，我们可以补充要求：
a1+⋅⋅⋅+ak=0,b1+⋅⋅⋅+bl=0(3.3)a_1+···+a_k=0,b_1+···+b_l=0 \qquad (3.3) a1+⋅⋅⋅+ak=0,b1+⋅⋅⋅+bl=0(3.3)
如果(3.3)(3.3)(3.3)式不成立，则分别把μ\muμ换为 μ+aˉ+bˉ\mu + \bar{a}+\bar{b}μ+aˉ+bˉ，aia_iai换为ai−aˉa_i-\bar{a}ai−aˉ，bjb_jbj换为bj−bˉb_j-\bar{b}bj−bˉ，则(3.1)(3.1)(3.1)式不变，而(3.3)(3.3)(3.3)式成立。

约束条件(3.3)(3.3)(3.3)给了ai，bja_i，b_jai，bj的意义一种更清晰的解释：ai>0a_i>0ai>0 表示A的水平iii的效应在AAA的全部水平的平均效应之上，ai<0a_i<0ai<0 则相反。另外，这个约束条件也给了μ，ai,bj\mu，a_i,b_jμ，ai,bj的一个适当的估计法：把YijY_{ij}Yij对一切i,ji,ji,j相加．注意到(3.3)(3.3)(3.3)，有
∑i=1k∑j=1lYij=klμ+∑i=1k∑j=1leij(3.4)\sum_{i=1}^{k}\sum_{j=1}^{l}Y_{ij}= kl\mu+\sum_{i=1}^{k}\sum_{j=1}^{l}e_{ij} \qquad (3.4) i=1∑kj=1∑lYij=klμ+i=1∑kj=1∑leij(3.4)
由(3.2)(3.2)(3.2)得，
Yˉ=∑i=1k∑j=1lYij/kl(3.5)\bar{Y}=\sum_{i=1}^{k}\sum_{j=1}^{l}Y_{ij}/kl \qquad (3.5) Yˉ=i=1∑kj=1∑lYij/kl(3.5)
是μ\muμ的一个无偏估计。其次，有
∑j=1lYij=lμ+la+∑j=1leij(3.6)\sum_{j=1}^{l}Y_{ij}=l\mu+la+\sum_{j=1}^{l}e_{ij} \qquad (3.6) j=1∑lYij=lμ+la+j=1∑leij(3.6)
于是，记
Yiˉ=∑j=1lYij/l,Yjˉ=∑i=1kYij/k(3.7)\bar{Y_i}=\sum_{j=1}^{l}Y_{ij}/l, \quad \bar{Y_j}=\sum_{i=1}^{k}Y_{ij}/k \qquad (3.7) Yiˉ=j=1∑lYij/l,Yjˉ=i=1∑kYij/k(3.7)
由(3.7)(3.7)(3.7)知，Yjˉ\bar{Y_j}Yjˉ为μ+ai\mu+a_iμ+ai的一个无偏估计。于是得到aia_iai的一个无偏估计为
ai^=Yiˉ−Yˉ,i=1,⋯,k(3.8)\hat{a_i}=\bar{Y_i}-\bar{Y}, i=1,\cdots,k \qquad(3.8) ai^=Yiˉ−Yˉ,i=1,⋯,k(3.8)
同理，
bj^=Yjˉ−Yˉ,j=1,⋯,l(3.9)\hat{b_j}=\bar{Y_j}-\bar{Y}, j=1,\cdots,l \qquad(3.9) bj^=Yjˉ−Yˉ,j=1,⋯,l(3.9)
ai^,bj^\hat{a_i},\hat{b_j}ai^,bj^适合约束条件(3.3)(3.3)(3.3)。

下面进行方差分析，要设法把总平方和
SS=∑i=1k∑j=1l(Yij−Yˉ)2SS=\sum_{i=1}^{k}\sum_{j=1}^{l}(Y_{ij}-\bar{Y})^2 SS=i=1∑kj=1∑l(Yij−Yˉ)2
分解为三部分：SSA,SSB,SSeSS_A,SS_B,SS_eSSA,SSB,SSe，分别表示因素A,BA,BA,B和随机误差的影响。这种分解的主要目的是假设检验：
H0A:a1=⋯=ak=0(3.10)H_{0A}:a_1=\cdots=a_k=0 \qquad(3.10) H0A:a1=⋯=ak=0(3.10)
和
H0B:b1=⋯=bk=0(3.11)H_{0B}:b_1=\cdots=b_k=0 \qquad(3.11) H0B:b1=⋯=bk=0(3.11)
H0AH_0AH0A成立表示因素AAA对指标其实无影响。在实际问题中，绝对无影响的场合少见，但如影响甚小以致被随机误差所掩盖时，这种影响事实上等于没有。因此，拿SSASS_ASSA和SSeSS_eSSe的比作为检验统计量正符合这一想法．

接下来讲一下方差分解的小技巧：
Yij−Yˉ=(Yiˉ−Yˉ)+(Yjˉ−Yˉ)+(Yij−Yiˉ−Yjˉ+Yˉ)Y_{ij}-\bar{Y}=(\bar{Y_i}-\bar{Y}) + (\bar{Y_j}-\bar{Y})+(Y_{ij}-\bar{Y_i}-\bar{Y_j}+\bar{Y}) Yij−Yˉ=(Yiˉ−Yˉ)+(Yjˉ−Yˉ)+(Yij−Yiˉ−Yjˉ+Yˉ)
两边平方，对i,ji,ji,j求和，结合约束条件(3.3)，注意到
∑i=1k(Yiˉ−Yˉ)=0，∑j=1l(Yjˉ−Yˉ)=0，\sum_{i=1}^{k}(\bar{Y_{i}}-\bar{Y})=0， \sum_{j=1}^{l}(\bar{Y_{j}}-\bar{Y})=0， i=1∑k(Yiˉ−Yˉ)=0，j=1∑l(Yjˉ−Yˉ)=0，

∑i=1k(Yij−Yiˉ−Yjˉ+Yˉ)=∑j=1l(Yij−Yiˉ−Yjˉ+Yˉ)=0\sum_{i=1}^{k}(Y_{ij}-\bar{Y_i}-\bar{Y_j}+\bar{Y})=\sum_{j=1}^{l}(Y_{ij}-\bar{Y_i}-\bar{Y_j}+\bar{Y})=0 i=1∑k(Yij−Yiˉ−Yjˉ+Yˉ)=j=1∑l(Yij−Yiˉ−Yjˉ+Yˉ)=0

即知所有交叉积之和皆为0，而得到
SS=l∑i=1k(Yiˉ−Yˉ)2+k∑j=1l(Yjˉ−Yˉ)2+∑i=1k∑j=1l(Yij−Yiˉ−Yjˉ+Yˉ)2=SSA+SSB+SSe(3.12)\begin{aligned} SS&=l\sum_{i=1}^{k}(\bar{Y_{i}}-\bar{Y})^2+k\sum_{j=1}^{l}(\bar{Y_{j}}-\bar{Y})^2+\sum_{i=1}^{k}\sum_{j=1}^{l}(Y_{ij}-\bar{Y_i}-\bar{Y_j}+\bar{Y})^2 \\ &=SS_A + SS_B + SS_e \qquad(3.12) \end{aligned} SS=li=1∑k(Yiˉ−Yˉ)2+kj=1∑l(Yjˉ−Yˉ)2+i=1∑kj=1∑l(Yij−Yiˉ−Yjˉ+Yˉ)2=SSA+SSB+SSe(3.12)

第一个平方和可以作为因素AAA的影响的衡量，从前述Yiˉ−Yˉ\bar{Y_{i}}-\bar{Y}Yiˉ−Yˉ作为 aia_iai的估计可以理解第二个平方和同理。至于第三个平方和可作为随机误差的影响这一点，直接看不甚明显。可以从两个角度去理解：在SSSSSS中去掉SSASS_ASSA 和SSBSS_BSSB后，剩余下的再没有其他系统性因素的影响，故只能作为SSeSS_eSSe。另外，由模型(3.1)(3.1)(3.1)及约束条件(3.3)(3.3)(3.3)，易知
Yij−Yiˉ−Yjˉ+Yˉ=(μ+ai+bj+eij)−(μ+ai+eiˉ)−(μ+bj+ejˉ)+(μ+eˉ)=eij−eiˉ−ejˉ+eˉ(3.13)\begin{aligned} Y_{ij}-\bar{Y_i}-\bar{Y_j}+\bar{Y} &= (\mu + a_i + b_j + e_{ij}) - (\mu + a_i + \bar{e_{i}}) - (\mu + b_j + \bar{e_{j}} ) + (\mu + \bar{e}) \\ &=e_{ij}-\bar{e_i}-\bar{e_j}+\bar{e} \qquad(3.13) \end{aligned} Yij−Yiˉ−Yjˉ+Yˉ=(μ+ai+bj+eij)−(μ+ai+eiˉ)−(μ+bj+ejˉ)+(μ+eˉ)=eij−eiˉ−ejˉ+eˉ(3.13)
这里面已经毫无μ,ai,bj\mu,a_i,b_jμ,ai,bj的影响，而只含随机误差。

得到分解式(3.12)(3.12)(3.12)后，我们就可以像单囚素情况那样，写出下面的方差分析表：
SSA,SSBSS_A , SS_BSSA,SSB 自由度分别为其水平数减去1，这一点与单因素情况相同．总和自由度为全部观察值数目klklkl减去1．剩下的就是误差平方和自由度：
(kl−1)−(k−1)−(l−1)=(k−1)(l−1)(kl - 1) - (k - 1) - (l - 1) = (k - 1) (l - 1) (kl−1)−(k−1)−(l−1)=(k−1)(l−1)
表3.1 双因素方差分析表

项目	SSSSSS	自由度	MSMSMS	FFF比	显著性
AAA	SSASS_ASSA	k−1k-1k−1	MSAMS_AMSA	MSA/MSeMS_A / MS_eMSA/MSe	, *, 或无
BBB	SSBSS_BSSB	l−1l-1l−1	MSBMS_BMSB	MSB/MSeMS_B / MS_eMSB/MSe
误差	SSeSS_eSSe	(k−1)(l−1)(k - 1) (l - 1)(k−1)(l−1)	MSeMS_eMSe
总和	SSSSSS	kl−1kl-1kl−1

还有一点要注意：在采纳模型(3.1)(3.1)(3.1)时，我们事实上引进了一种假定，即两因素A,BA,BA,B对指标的效应是可以叠加的．换一种方式说：因素AAA的各水平的优劣比较，与因素BBB处在哪个水平无关，反之亦然．更一般的情况是：A,BA,BA,B两因子有“交互作用＂。这时在模型(5.13)中，还要加上表示交互作用的项cijc_{ij}cij．这时不仅统计分析复杂化了，尤其是分析结果的解释也复杂化了．本文档暂不讨论这种情况。在一个特定的问题中，交互作用是否需要考虑，在很大程度上取决于问题的实际背景和经验．有时，通过试验数据的分析也可以看出一些问题。例如，若误差方差σ2\sigma^2σ2的估计MSeMS_eMSe反常地大，则有可能是由于交互作用所致．因为可以证明：若交互作用确实存在而未加考虑，则它的影响进入随机误差而增大了MSeMS_eMSe。

2 代码实现

由于目前还没有接触过R，因此这部分的坑之后再来填，给出spss的实现过程和结果

最后

第一次参与了DW的组队学习，我觉得非常好，但是发现了自己理论和实践上都存在很大程度的不足，最近正好期末考试了，不能全身心投入学习和充电的工作，希望未来可以进一步完善自己该部分的内容，感谢组织！