方差分析
什么是方差分析
方差分析（Analysis of Variance，ANOVA）由R.A.Fisher发明.听起来好像是研究参数方差，其实用于两个及两个以上样本均数差别的显著性检验。
为什么要用方差分析
课堂中老师举了一个例子：一个律师怀疑自己的代理人的陪审团中的女性所在比例太少。为了检验这个猜想， 律师找了本案法官和其他6名法官近几次来开庭时陪审团成员信息，见下图：
虽然一看就知道，律师的代理人陪审团成员中女性的比例的确少于其他法官，但是毕竟是抽样调查，不做点统计检验没有说服力。所以他需要：
证明代理人的陪审团女性成员与其他相比显著性的少，要进行6次T检验。
证明其他法官中陪审团女性成员比例间没有显著性差异，要进行次T检验。
由于每次统计检验中都有一定概率犯I型错误，所以当检验次数增加时，就不可避免就会提高I型错误概率。如果想通过P值修正，降低I型错误，那么II型错误就会增加，功效机会降低。
有没有办法能够一次性进行多个样本间的比较呢？方法就是本次文章要介绍的方差分析。
方差分析的前提
不过方差分析不是随随便便就能用的，至少要满足以下几个条件：
独立方差性：样本必须来自于一个正态分布总体，并且样本间是相互独立的
等方差性或方差齐性：抽样的总体必须是等方差的。
方差分析的基本思想
引起观测值不同（波动）的原因主要有两类：一类是实验过程中随机因素的干扰或观测误差引起的不可控制的波动；另一类是由于试验中处理方式不同或试验条件不同引起的可以控制的波动。
即
方差分析的主要工作就是讲观测数据的总变异（波动）按照变异的原因的不同分解为因子效应与试验误差，并对其做出数量分析，比较各种原因在总体变异种所占的重要程度，以此作为进一步统计推断的依据。
方差分析的假设
空假设（Null hypothesis): 所有比较组的总体均值相同
备择假设（alternative hypothesis): 至少有一个均值与其他不同
方差分析的术语
目前有两种方法治疗焦虑症：认知行为疗法（CBT)和眼动脱敏再加工法（EMDR)。为了验证哪类资料方法效果好，一共找了10位焦虑症志愿者，随机分为两组，各接受五周的实验。然后要求每位患者填写状态特质焦虑问卷（STAI)。
实验设计表如下：
这次实验中治疗方案是自变量(independent variable)，STAI是应变量(dependent variable or response variable).
治疗方案是两水平（CBT,EMDR)的组间因子（组间因子,就是每位患者都只接受一种治疗方案）
此外每种治疗方案的观测数量相等，所以是均衡设计，观测数不同就是非均衡设计了。
因为均有一个类别型变量，也就是指只有治疗方案一个变量，所以就可以称之为单因素方差分析,one-way ANOVA.
同样是10个志愿者，如果他们都接受CBT治疗，但是需要在5周和6个月后填写STAI，那么这个就是单因素组内方差分析。因为存在组内因子—时间，即每位患者都要有5周和6个月的资料报告，这与之前患者要么CBT,要么EMDR不同。又因为每个志愿者要进行两次测量，也可以称之为重复测量方差分析
如果需要同时研究治疗方案和时间对资料效果的影响，那么就是双因素(two-way ANOVA)。
实验设计包含两个及以上因子时，就称为因素方差分析设计。如果因子设计包括组内和组间因子，称为混合模式方差分析。
此外，我们还发现抑郁症对病症的治疗也有影响，并且抑郁症和焦虑症常常同时发作。所以尽管志愿者是随机分组，但是不同的抑郁水平也会导致最终的治疗效果存在差异。这里的抑郁症就被称为混淆因素(confounding factor），由于本实验不是为了研究抑郁症，所以它称为干扰变数（nuisance variable)
进一步，假设志愿者使用了抑郁症的自我评价报告，比如XX，记录了他们的抑郁水平，那么你可以在评价疗法类型的影响前，对任何影响抑郁水平的组内差异进行统计学调整。XX为协变量，该设计为协变量方差分析
再进一步，以上设计只记录了单个因变量的情况，为了增强研究的有效性，我们还可以对焦虑症进行其他的测量（如家庭评分，焦虑症对日常行为的影响评价）。所以因变量也就增加了。这个时候实验设计就被称为多元方差分析（MANOVA），如果协变量也存在，就叫多元协变量方差分析（MANCOVA)。
小结一下之前出现的名词：
自变量，因变量，协变量，混淆因素
均衡设计，非均衡设计
组内变量，组间变量
单因素方差分析，双因素方差分析，重复测量方差分析，因素方差分析统计，混合模型方差分析
混淆因素， 干扰变数
协变量方差分析，多元方差分析，多元协变量方差分析
好了，你可以拿这些名词去吓人了。
如何用R进行方差分析
课堂或者课本对于最简单的单因素方差检验，都会介绍如何计算SS,df, MS, 和F比例，最后得出一个ANOVA表：
但我们只需要知道原理就好了，毕竟有很多软件能够进行ANOVA分析了，这里就是介绍如何用R进行方差分析。
R语言主要使用aov()函数进行方差分析，语法为aov(formula, data=dataframe).下表来自于《R语言实战》的总结：
单因素方差分析
以multcomp包中cholesterol数据作为案例进行演示。原始数据是用来研究不同药物处理下，降低胆固醇的量。不过你也可以理解为不同施肥下的产量，不同处理下，同一个基因的表达量等等。
注:>表明输入，没有> 结果
加载数据
> install.packages("multcomp")
> library("multcomp")# 加载数据集到环境变量中，就可以不用$选择某一列数据了> attach(cholesterol)探索性分析
简单对数据进行观察
> str(cholesterol)'data.frame': 50 obs. of 2 variables:
$ trt : Factor w/ 5 levels "1time","2times",..: 1 1 1 1 1 1 1 1 1 1 ...
$ response: num 3.86 10.39 5.91 3.06 7.72 ...说明’choleserol’是一个数据框，包括一个因子’trt’和数值变量’response’
接下来，具体看下这两个变量的描述性统计结果
> summary(cholesterol)
trt response
1time :10 Min. : 2.304
2times:10 1st Qu.: 8.409
4times:10 Median :12.605
drugD :10 Mean :12.738
drugE :10 3rd Qu.:17.519
Max. :27.244不难发现’trt’分为5个处理，每个处理都是10个样本。不过’response’只是反映了总体的情况，最好是能知道不同处理的四分位值和方差等信息，这些可以使用aggregate()根据因子计算得到。
# aggregate> aggregate(response,by=list(trt),summary)
Group.1 x.Min. x.1st Qu. x.Median x.Mean x.3rd Qu. x.Max.1 1time 2.304 3.261 5.415 5.782 7.673 10.3902 2times 3.505 7.914 8.615 9.225 10.110 15.8303 4times 8.053 10.550 12.610 12.370 13.950 18.1804 drugD 8.819 13.870 16.220 15.360 17.620 19.9805 drugE 15.770 18.830 21.020 20.950 22.560 27.240# 其实做一个箱图更加直观> boxplot(response~trt)# 了解一下标准差> aggregate(response,by=list(trt),FUN=sd)
Group.1 x1 1time 2.8781132 2times 3.4830543 4times 2.9231194 drugD 3.4546365 drugE 3.345003基本信息了解的差不多了，就可以进行方差分析了。在方差分析之前，先要检验一下是否符合方差分析的前提：独立正态性，同方差性。
其中最重要的就是同方差性，常用两个方法
Bartlett检验
Levene检验（由car包提供）
# Bartlett检验> bartlett.test(response ~ trt,data = cholesterol)
Bartlett test of homogeneity of variances
data: response by trt
Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
# Levene检验,对原始数据的正态性不敏感
> library(car)
> leveneTest(response ~ trt,data = cholesterol)
Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 4 0.0755 0.9893
45无论Bartlett检验还是Levene检验，两者的P值都大于0.05，因此接受原假设：样本之间的方差是相同的。因此可以接着做方差分析了。
注意：如果发现不是样本件非等方差，那么就需要用Welch的ANOVA分析法.函数为oneway.test()，参数为var.equal=FALSE
方差分析
# 使用aov()进行方差分析> 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看这个P值就知道样本之间的均值存在显著性差异的。那么问题来了，我们需要进一步找到显著性差异的那个样本。比如说前面案例中就是证明代理人的法官的陪审团女性比例的确少了。
如果你想判断出是不是某个特定药物，特定处理或者样本是否和其他样本有所差别，比如说案例中律师可能只需要知道代理人的法官女性成员比其他法官有显著性少，那么只需要
排除特定样本，对剩余样本进行ANOVA，观察是否有差异
将除特定样本外的其余样本放在一起，然后和特定样本进行比较
多重比较
但实际上，我们是没有特定目标的，比如想确定到底哪个药物效果或者哪种治疗方案最好。因此需要对均值做多重比较。这里介绍两种方法：
多重t检验：
说明：多重t检验使用方便，但当多次重复使用t检验时会增加I类错误的概率，从而使得“有显著差异”的结论不一定可靠，所以在进行较多次重复比较时候，要对P值进行调整。
# 查看调整的所用方法> p.adjust.methods
[1] "holm" "hochberg" "hommel" "bonferroni" "BH" "BY" "fdr" "none"# 多重比较# 不修正> pairwise.t.test(response,trt,p.adjust.method = "none")
Pairwise comparisons using t tests with pooled SD
data: response and trt
1time 2times 4times drugD
2times 0.02133 - - -
4times 3.8e-05 0.03435 - -
drugD 3.5e-08 0.00011 0.04432 -
drugE 1.1e-13 2.3e-10 3.8e-07 0.00035P value adjustment method: none#使用fdr修正> pairwise.t.test(response,trt,p.adjust.method = "fdr")
Pairwise comparisons using t tests with pooled SD
data: response and trt
1time 2times 4times drugD
2times 0.02667 - - -
4times 7.6e-05 0.03817 - -
drugD 1.2e-07 0.00018 0.04432 -
drugE 1.1e-12 1.1e-09 9.6e-07 0.00050P value adjustment method: fdr同时置信区间：Tukey法
Tukey区间计算比较复杂，但是克服了多重t检验不修正时存在过多显著性差异，修正后可能又太少的问题。使用方法很简单
# 使用方法TuekyHSD(aov(response~trt))# 可视化plot(TukeyHSD(aov(response~trt)))这个图主要看虚线是否在成对比较的置信区间。置信区间在线外表明该成对比较具有显著性差异。
此外可用使用multcomp包的glht()函数。
PS： 以上都是对于正态分布，并且是可度量的变量，如果遇到名义型变量(ordinal)，那么就需要用到非参数ANOVA—Kruskal-Wallis。和上述过程一样，只不过使用
kruskal.test进行方差分析，使用pairwise.wilcox.test()进行多重比较
单因素方差分析小结
双因素方差分析
双因素的方差分析，其基本思想和方法和单因素的方差分析相似。差异在于双因素方差分析中可能会出现交互作用。
双因素ANOVA（有重复）,可研究交互作用
双因素ANOVA数学模型如下：
也就是说总变异可以分解为因素A的变异，因素B的变异，A和B交互导致的变异和误差导致的变异。即
检验的原假设（空假设）为：
相应的双因素ANOVA表：
以R语言自带的ToothGrowth数据集为例，随机分配60只豚鼠，分别采用两种喂食方法（橙汁或维生素C），各喂食方法中抗血酸含量有三种水平（0.5 mg, 1mg或 2mg），每种处理方式组合都被分配10中豚鼠。牙齿长度为因变量。
方差分析
# 加载数据attach(ToothGrowth)# 探索性分析，均衡设计？每组的均值和方差是多少？table(supp, dose)
aggregate(len, by=list(supp,dose), FUN=mean)
aggregate(len, by=list(supp,dose), FUN=sd)## 方差齐性检验bartlett.test(len ~ supp)
bartlett.test(len ~ dose)# 双因素ANOVA> 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从方差分析表中可以发现，主效应（supp和dose）和交互效应都非常显著。
探索交互作用
进一步，我们可以用interaction.plot可视化观察下交互效应
interaction.plot(dose, supp, len, type="b",
col=c("red","blue"),pch=c(16,18),
main="Interaction between Dose and Supplement Type")gplots包中的plotmeans()函数也可以展示交互效应，使用help研究吧，这里不赘述。
但是HH包的interaction2wt()函数一定要在这里展示一下,因为他能对任意顺序的因子设计的主效应和交互效应都会进行展示
library(HH)
interaction2wt(len~supp*dose)通过可视化可以直观的发现，牙齿长度会随着橙汁和维生素C的抗坏血酸的增加而变长。
多重比较
下一步也就是多重比较，也非常简单就是用TukeyHSD()函数, 寻找
TukeyHSD(fit)
plot(TukeyHSD(fit))双因素方差分析，无重复
没有重复的后果就是不能研究因素间的交互作用，有时候实验经费不够的话也就只能这样了。
数学模型为：
原假设为：
R语言中aov函数为：
aov(y ~ B + A, data)注意这里的顺序很重要，B表示为区组因子。假设你要研究不同施肥对作物产量的影响，一般而言你种植地的土壤本身就不同，那么不同突然就是这里的区间因子。或者你要研究不同药物不同剂量的效果，那么不同剂量就是区间因子。
后面流程基本和双因素方差分析（含重复）一致，不多说。
最后注意
以上都是基于固定效应模型进行ANOVA分析，随机效应和混合效应没有讨论，等待我下次介绍这些区别吧。可能还会介绍下ANOVA分析，回归和相关之间的关系。
参考资料
生物统计学课程PPT from Xingming Zhao
R语言实战
R语言与统计分析