4. 方差分析

方差分析

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 推导过程

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

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

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

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

接下来，把$SS$分为两部分，一部分表示随机误差的影响，记为$S{S}_e$；另一部分表示因素$A$的各水平理论平均值$a_i$不同带来的影响，记为$S{S}_A$。

关于$S{S}_e$，先固定一个$i$，此时对应的所有观测值$Y_{i1},Y_{i2},\cdots ,Y_{in}$，他们之间的差异与每个水平的理论平均值不等无关，而是取决于随机误差，反映这些观察值差异程度的量是${\sum }_{j=1}^{n_i}{\left(Y_{ij}-\overline{Y_i}\right)}^2$，其中
Y i ˉ = ( Y i 1 + Y i 2 + ⋯ + Y i n ) / n i , i = 1 , 2 ,