方差分析ANOVA:理论、推导与R语言实现
方差分析
1 概要
方差分析(Analysis of variance, ANOVA) 主要研究分类变量作为自变量时，对因变量的影响是否是显著的。
方差分析的方法是由20世纪的统计学家Ronald Aylmer Fisher在1918年到1925年之间提出并陆续完善起来的，该方法刚开始是用于解决田间实验的数据分析问题，因此，方差分析的学习是和实验设计、实验数据的分析密不可分的。
实验设计和方差分析都有自己相应的语言。因此，在这里我们通过一个焦虑症治疗的实例，先了解一些术语，并且思考一下，方差分析主要用于解决什么样的问题。
以焦虑症治疗为例，现有两种治疗方案：认知行为疗法（CBT）和眼动脱敏再加工法（EMDR）。我们招募10位焦虑症患者作为志愿者，随机分配一半的人接受为期五周的CBT，另外一半接受为期五周的EMDR，设计方案如表1-1所示。在治疗结束时，要求每位患者都填写状态特质焦虑问卷（STAI），也就是一份焦虑度测量的自我评测报告。
表1-1 单因素组间方差分析
| CBT | EMDR |
|---|---|
| s1 | s6 |
| s2 | s7 |
| s3 | s8 |
| s4 | s9 |
| s5 | s10 |
在这个实验设计中，治疗方案是两水平（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 | ||
| s7 | ||
| s8 | ||
| s9 | ||
| s10 |
此时，时间（time）是两水平（五周、六个月）的组内因子。因为每位患者在所有水平下都进行了测量，所以这种统计设计称单因素组内方差分析；又由于每个受试者都不止一次被测量，也称作重复测量方差分析。当时间的F检验显著时，说明患者的STAI得分均值在五周和六个月间发生了改变。
现假设你对治疗方案差异和它随时间的改变都感兴趣，则将两个设计结合起来即可：随机分配五位患者到CBT，另外五位到EMDR，在五周和六个月后分别评价他们的STAI结果（见表1-3）。
表1-3 含组间和组内因子的双因素方差分析
| 时间 | |||
|---|---|---|---|
| 疗法 | 患者 | 5周 | 6个月 |
| CBT | s1 | ||
| s2 | |||
| s3 | |||
| s4 | |||
| s5 | |||
| EMDR | s6 | ||
| s7 | |||
| s8 | |||
| s9 | |||
| s10 |
疗法（therapy）和时间（time）都作为因子时，我们既可分析疗法的影响（时间跨度上的平均）和时间的影响（疗法类型跨度上的平均），又可分析疗法和时间的交互影响。前两个称作主效应，交互部分称作交互效应。
当设计包含两个甚至更多的因子时，便是因素方差分析设计，比如两因子时称作双因素方差分析，三因子时称作三因素方差分析，以此类推。若因子设计包括组内和组间因子，又称作混合模型方差分析，当前的例子就是典型的双因素混合模型方差分析。
本例中，你将做三次F检验：疗法因素一次，时间因素一次，两者交互因素一次。若疗法结果显著，说明CBT和EMDR对焦虑症的治疗效果不同；若时间结果显著，说明焦虑度从五周到六个月发生了变化；若两者交互效应显著，说明两种疗法随着时间变化对焦虑症治疗影响不同（也就是说，焦虑度从五周到六个月的改变程度在两种疗法间是不同的）。
现在，我们对上面的实验设计稍微做些扩展。众所周知，抑郁症对病症治疗有影响，而且抑郁症和焦虑症常常同时出现。即使受试者被随机分配到不同的治疗方案中，在研究开始时，两组疗法中的患者抑郁水平就可能不同，任何治疗后的差异都有可能是最初的抑郁水平不同导致的，而不是由于实验的操作问题。抑郁症也可以解释因变量的组间差异，因此它常称为混淆因素（confounding factor）。由于你对抑郁症不感兴趣，它也被称作干扰变数（nuisance variable）。
假设招募患者时使用抑郁症的自我评测报告，比如白氏抑郁症量表（BDI），记录了他们的抑郁水平，那么你可以在评测疗法类型的影响前，对任何抑郁水平的组间差异进行统计性调整。本案例中，BDI为协变量，该设计为协方差分析（ANCOVA）。
以上设计只记录了单个因变量情况（STAI），为增强研究的有效性，可以对焦虑症进行其他的测量（比如家庭评分、医师评分，以及焦虑症对日常行为的影响评价）。当因变量不止一个时，设计被称作多元方差分析（MANOVA）， 若协变量也存在， 那么就叫多元协方差分析（MANCOVA）。
下面我们主要介绍单因素方差分析与双因素方差分析的原理与实现。
2 单因素方差分析
2.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−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)i=1∑kj=1∑ni(Yiˉ−Yˉ)2−i=1∑kj=1∑ni(Yij−Yiˉ)2
KaTeX：对齐由原来的\begin{align}和\end{align} 目前改为 \begin{aligned}和\end{aligned}
因为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.2 代码实例
单因素方差分析的R语言实现
单因素方差分析中，你感兴趣的是比较分类因子定义的两个或多个组别中的因变量均值。以multcomp包中的cholesterol数据集为例，50个患者均接受降低胆固醇药物治疗（trt）五种疗法中的一种疗法。其中三种治疗条件使用药物相同，分别是20mg一天一次（1time）、10mg一天两次（2times）和5mg一天四次（4times）。剩下的两种方式（drugD和drugE）代表候选药物。
> library(multcomp)
> attach(cholesterol)
>
> # 统计各组样本大小
> table(trt)
trt1time 2times 4times drugD drugE 10 10 10 10 10
>
> # 各组均值
> aggregate(response, by=list(trt), FUN=mean)Group.1 x
1 1time 5.78197
2 2times 9.22497
3 4times 12.37478
4 drugD 15.36117
5 drugE 20.94752
>
> # 各组标准差
> aggregate(response, by=list(trt), FUN=sd)Group.1 x
1 1time 2.878113
2 2times 3.483054
3 4times 2.923119
4 drugD 3.454636
5 drugE 3.345003
>
> # 进行方差分析
> fit <- aov(response ~ trt)
> summary(fit)Df Sum Sq Mean Sq F value Pr(>F)
trt 4 1351.4 337.8 32.43 9.82e-13 ***
Residuals 45 468.8 10.4
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
方差分析的结果中，各项数字的含义可以参照表2-1。
查看各水平对应的组均值的差异
gplots包中的plotmeans()可以用来绘制带有置信区间的组均值图形。如图9-1所示，图形展示了带有95%的置信区间的各疗法均值，可以清楚看到它们之间的差异。
library(gplots)
plotmeans(response ~ trt, xlab="Treatment", ylab="Response",main="Mean Plot\nwith 95% CI")
detach(cholesterol)
2-1 五种降低胆固醇药物疗法的均值，含95%的置信区间
多重比较
虽然ANOVA对各疗法的F检验表明五种药物疗法效果不同，但是并没有告诉你哪种疗法与其他疗法不同。多重比较可以解决这个问题。例如，TukeyHSD()函数提供了对各组均值差异的成对检验。
> TukeyHSD(fit)Tukey multiple comparisons of means95% family-wise confidence levelFit: aov(formula = response ~ trt)$trtdiff lwr upr p adj
2times-1time 3.44300 -0.6582817 7.544282 0.1380949
4times-1time 6.59281 2.4915283 10.694092 0.0003542
drugD-1time 9.57920 5.4779183 13.680482 0.0000003
drugE-1time 15.16555 11.0642683 19.266832 0.0000000
4times-2times 3.14981 -0.9514717 7.251092 0.2050382
drugD-2times 6.13620 2.0349183 10.237482 0.0009611
drugE-2times 11.72255 7.6212683 15.823832 0.0000000
drugD-4times 2.98639 -1.1148917 7.087672 0.2512446
drugE-4times 8.57274 4.4714583 12.674022 0.0000037
drugE-drugD 5.58635 1.4850683 9.687632 0.0030633> par(las=2)
> par(mar=c(5,8,4,2))
> plot(TukeyHSD(fit))
成对比较图形如图2-2所示。图形中置信区间包含0的疗法说明差异不显著（p>0.05）。
图2-2 Tukey HSD均值成对比较图
评估检验的假设条件
根据2.1中我们讲的关于方差分析的推导中，我们知道，方差分析结果的有效性是建立在一系列假设条件之上的，因此，在我们使用方差分析模型时，需要评估进行方差分析的数据，是否符合模型使用的假设条件。
正态性检验
第一，在建立模型时，我们假设因变量是服从正态分布的，需要进行正态性检验。
正态性检验的方法有两种，一是通过QQ图进行检验。
# QQ plot
library(car)
qqPlot(lm(response ~ trt, data=cholesterol),simulate=TRUE, main="Q-Q Plot", labels=FALSE)
除此之外，R里面也提供了一些package来进行正态性检验。
K-S test
统计学里, Kolmogorov–Smirnov 检验(亦称：K–S 检验)是用来检验数据是否符合某种分布的一种非参数检验。其原假设H0H_0H0:两个数据分布一致或者数据符合理论分布。在R语言里，我们可以使用ks.test(x, pnorm)进行正态性检验，若结果中的p值大于0.05，则数据符合正态分布。
Anderson–Darling test
Anderson–Darling检验是一种用来检验给定的样本是否来自于某个确定的概率分布的统计检验方法。在R语言中，我们可以从nortest包中的ad.test()进行检验。若结果中的p值大于0.05，则数据符合正态分布。
Shapiro-Wilk test
Shapiro-Wilk检验在小样本情况下，是很普通的正态性检验方法，Shapiro.test()在默认安装的stats包中。原假设H0H_0H0: 数据符合正态分布。
Lilliefor test
Lilliefor test是基于Kolmogorov–Smirnov test的一种正态性检验。原假设H0H_0H0: 数据符合正态分布，lillie.test()也在nortest包中。
方差齐性检验
因为方差分析的实质是检验多个水平的均值是否有显著差异，如果各个水平的观察值方差差异太大，只检验均值之间的差异就没有意义了，所以要进行方差齐性检验。
Bartlett test可以用来检验数据的方差齐性。
> bartlett.test(response ~ trt, data=cholesterol)Bartlett test of homogeneity of variancesdata: response by trt
Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
Bartlett检验表明五组的方差并没有显著不同（p=0.97）。其他检验如Fligner-Killeen检验
（fligner.test()函数）和Brown-Forsythe检验（HH包中的hov()函数）此处没有做演示，但它们获得的结果与Bartlett检验相同。
不过，方差齐性分析对离群点非常敏感。可利用car包中的outlierTest()函数来检测离群点：
outlierTest(fit)
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:rstudent unadjusted p-value Bonferroni p
19 2.251149 0.029422 NA
从输出结果来看，并没有证据说明胆固醇数据中含有离群点（当p>1时将产生NA）。因此根据正态性检验、方差齐性检验和离群点检验，该数据似乎可以用ANOVA模型拟合得很好。这些方法反过来增强了我们对于所得结果的信心。
3 双因素方差分析
3.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。
3.2 代码实例
在双因素方差分析中，受试者被分配到两因子的交叉类别组中。以基础安装中的ToothGrowth数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗坏血酸含量有三种水平（0.5mg、1mg或2mg），每种处理方式组合都被分配10只豚鼠，牙齿长度为因变量。
> attach(ToothGrowth)
> table(supp, dose)dose
supp 0.5 1 2OJ 10 10 10VC 10 10 10
>
> aggregate(len, by=list(supp, dose), FUN=mean)Group.1 Group.2 x
1 OJ 0.5 13.23
2 VC 0.5 7.98
3 OJ 1.0 22.70
4 VC 1.0 16.77
5 OJ 2.0 26.06
6 VC 2.0 26.14
>
> aggregate(len, by=list(supp, dose), FUN=sd)Group.1 Group.2 x
1 OJ 0.5 4.459709
2 VC 0.5 2.746634
3 OJ 1.0 3.910953
4 VC 1.0 2.515309
5 OJ 2.0 2.655058
6 VC 2.0 4.797731
>
> dose <- factor(dose)
#dose变量被转换为因子变量，这样aov()函数就会将它当做一个分组变量，而不是一个数值型协变量
> # condider interactive factor
> fit <- aov(len ~ supp*dose)
> summary(fit)Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1计算结果表明，主效应和交互效应都是显著的。
有多种方式对结果进行可视化处理。此处可用interaction.plot()函数来展示双因素方
差分析的交互效应。
# interactive effect
interaction.plot(dose, supp, len, type="b",col=c("red","blue"), pch=c(16, 18),main = "Interaction between Dose and Supplement Type")
图3-1 各种剂量喂食下豚鼠牙齿长度的均值（interaction.plot()函数绘制）
还可以用gplots包中的plotmeans()函数来展示交互效应。
图3-2 喂食方法和剂量对牙齿生长的交互作用。用plotmeans()函数绘制的95%的置
信区间的牙齿长度均值
图形展示了均值、误差棒（95%的置信区间）和样本大小。
最后，你还能用HH包中的interaction2wt()函数来可视化结果，图形对任意顺序的因子设计的主效应和交互效应都会进行展示（图3-3）。
library(HH)
interaction2wt(len~supp*dose)
>
> detach(ToothGrowth)
图3-3 ToothGrowth数据集的主效应和交互效应。图形由interaction2wt()函数创建
以上三幅图形都表明随着橙汁和维生素C中的抗坏血酸剂量的增加，牙齿长度变长。对于0.5mg和1mg剂量，橙汁比维生素C更能促进牙齿生长；对于2mg剂量的抗坏血酸，两种喂食方法下牙齿长度增长相同。
参考书目：
- 陈希孺，概率论与数理统计
- Robert I. Kabacoff, R in Action.
来源：https://github.com/datawhalechina/team-learning-data-mining/tree/master/ProbabilityStatistics
方差分析ANOVA:理论、推导与R语言实现相关推荐
- r语言算巢式设计方差分析_应用统计学与R语言实现学习笔记（八）——方差分析...
Chapter 8 ANOVA 本篇是第八章,内容是方差分析.前一段考试,汇报,作业.忙不过来,停更了一段时间,现在重新开始更这一部分内容.方差分析是很多实验的基础以及很重要的分析手段,这一章内容相比 ...
- 【视频】极值理论EVT与R语言应用：GPD模型火灾损失分布分析
最近我们被客户要求撰写关于极值理论的研究报告,包括一些图形和统计输出.正态分布属于统计学里的知识,对于我们科研来说在数据处理时常常用到所以需要学习相关的知识. 正态分布在自然界中是一种最常见的分布.例 ...
- R语言单因素重复测量方差分析（one-way repeated measures ANOVA）实战
R语言单因素重复测量方差分析(one-way repeated measures ANOVA)实战 目录 R语言单因素重复测量方差分析(one-way repeated measures ANOVA) ...
- R语言嵌套方差分析（Nested ANOVA）实战
R语言嵌套方差分析(Nested ANOVA)实战 目录 R语言嵌套方差分析(Nested ANOVA)实战 #嵌套方差分析(Nested ANOVA)
- R语言方差分析的注意事项
本文首发于公众号:医学和生信笔记,完美观看体验请至公众号查看本文. 文章目录 均衡设计和非均衡设计 方差分析的3种类型 示例 one-way anova two-way anova 协方差分析 R语言 ...
- R语言极值理论 EVT、POT超阈值、GARCH 模型分析股票指数VaR、条件CVaR：多元化投资组合预测风险测度分析
最近我们被客户要求撰写关于极值理论的研究报告,包括一些图形和统计输出. 相关视频: 极值理论EVT与R语言应用:GPD模型火灾损失分布分析 R语言极值理论EVT:基于GPD模型的火灾损失分布分析 相关 ...
- R语言差异检验：单因素方差分析
文章目录 @[toc] 方差分析介绍 适用条件 分类 R语言 单因素方差分析示例 数据集 示例 多重比较 评估检验的假设条件 t检验可以解决单样本.双样本时的均数比较.当要比较的组多于两个时,t检验方 ...
- 数据分析-R语言资料整理
独家分享--48页PPT解密数据可视化! Excel图表快捷操作小技巧 基于随机森林的分类与回归 R语言制作网页 ggplot2:可视化设计师的神器,了解一下 [译]R包介绍:Online Rando ...
- 精心整理 | R语言中文社区历史文章整理（类型篇）
2018年过去一半了~又到了盘点的时间~感谢长时间来各位好友的关注,我们的成长与你们的爱护是分不开的.更感谢各位老师的投稿,支撑起了我们的这个社区,让更多R语言的爱好者和从业者获得最棒的知识!本文选取 ...
最新文章
- Eclipse Theme
- python使用教程cmd啥意思-对python中执行DOS命令的3种方法总结
- ASP.NET MVC编程——视图
- 人工智能到底威胁人类还是造福人类？
- hue安装及基本测试-笔记
- 移动开发之我见--“Android开发生涯”
- matlab 字符串处理（单引号、拼接、char）
- 小程序web开发框架－weweb介绍 1
- 团队项目:VS2013和SQL Server2012的连接使用
- 3.通信原理——随机过程（第七版 樊昌信 曹丽娜编著）
- 高通QFIL 导出所有分区
- 图解机器学习：人人都能懂的算法原理
- sobol灵敏度分析matlab_灵敏度分析 使用MATLAB编写
- 决策树(2)——CART算法
- 【游戏逆向】老飞飞怀恋魅力爱玩等老飞飞瞬移分析代码
- 月薪4万是一种什么样的感受？
- Shell循环语句（for、while、until）及echo、IFS
- Transact-SQL参考：sp_who
- 牛客多校第十场F-Popping Balloons
- 域名可以转让注册人吗_别人帮注册的域名怎么过户