上一篇文章说了方差和t检验的差异，这篇说说几种实用的方差分析方法和R语言实现。
一般情况下，基本的方差分析模型包含以下三类，三类下面会根据具体情况再进行细分，主要的三类为一元方差分析，协方差分析，多元方差分析。
1、一元方差分析
一元方差分为单因素、多因素两类（协方差单独分类），既然方差是检验各组差异的，那么从一个最简单的例子入手，探寻各类方差分析的适用条件和特点。
OK，正题开始，鉴于自己也算是酷爱篮球，就举个篮球运动员的例子吧，以下数据纯属瞎编，如有雷同纯属巧合。
有一天，詹姆斯和科比遇见了，鉴于科比已经退役，詹姆斯说老科，我们比投篮吧，看看是你厉害还是我厉害。科比一听这话顿时来了精神，就在我家球场，放学别走。
詹姆斯说，我的投篮命中率应该比你高，但是为了避免你蒙进去的多，我们比试5组，每一组投10次，看谁进的多，科比：Deal！
下面是相关代码
name<-c(rep('Kobe',5),rep('James',5))
good<-c(5,4,6,5,7,7,5,6,4,8)
a<-data.frame(name,good)
fit1<-aov(good~name,a)
summary(fit1)其中建立了一个数据框：a，也就是两个人的投篮比较：
name good
1 Kobe 5
2 Kobe 4
3 Kobe 6
4 Kobe 5
5 Kobe 7
6 James 7
7 James 5
8 James 6
9 James 4
10 James 8结果如下：
Df Sum Sq Mean Sq F value Pr(>F)
name 1 0.9 0.9 0.474 0.511
Residuals 8 15.2 1.9 可以看到，p=0.511，也就是说，二者的投篮没有显著差异（同样可以使用t检验）。
这时候，乔丹老流氓来了，看到他们在比投篮，说，要不我也参与下，好久没活动了。
随着老流氓的连续空心，科比和詹姆斯马上就要报警了，最后看看三人的命中数：
b<-rbind(a,data.frame(name=rep('jordan',5),good=c(10,9,10,10,10)))即
name good
1 Kobe 5
2 Kobe 4
3 Kobe 6
4 Kobe 5
5 Kobe 7
6 James 7
7 James 5
8 James 6
9 James 4
10 James 8
11 jordan 10
12 jordan 9
13 jordan 10
14 jordan 10
15 jordan 10接下来就做一个t检验做不了的事情：两组以上的方差分析：
fit2<-aov(good~name,b)
summary(fit2)结果是：
Df Sum Sq Mean Sq F value Pr(>F)
name 2 56.93 28.467 21.35 0.000111 ***
Residuals 12 16.00 1.333
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1我们看到P值是0.000111，也就是超级显著了，但是我们知道三者有显著差别，但是谁和谁有显著差别暂时还不知道（虽然能猜到，但还是装不知道），这时候需要借助Tukey或Duncan方法进行检验下，这里先使用Tukey方法检测：
TukeyHSD(fit2)结果为
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = good ~ name, data = b)
$name
diff lwr upr p adj
Kobe-James -0.6 -2.548332 1.348332 0.6973146
jordan-James 3.8 1.851668 5.748332 0.0005973
jordan-Kobe 4.4 2.451668 6.348332 0.0001637可以看到，Kobe和James没有差异，乔老爷子和两个人都有差异，但是和Kobe差异最大，这个用图形表示就是
plot(TukeyHSD(fit2))
图形中纵轴是均值的差异，包含0说明差异不明显，不包含0说明差异明显，可以看到乔丹的均值比詹姆斯和科比都大很多。
那么这个结果说明了乔丹比科比和詹姆斯投篮准，最后还需要进行正态性假设，因为我们需要知道这5组投篮是随机从他们职业生涯的千千万万个投篮中抽取的，这里需要检测下残差是否符合正态分布。
到这里说的有点远，因为残差需要进行线性拟合，但是有三个因素性变量，但是我不能把名字作为横轴吧，那么乔丹=？科比=？，一般在这种情况下处理就是两两比较，比如把詹姆斯的投篮命中作为因变量，乔丹作为自变量1，科比作为自变量2，那么其实舍弃了科比和乔丹的比较，比如载入car包中，看具体的因子编码：
library(car)
contrasts(b$name)结果：
Kobe jordan
James 0 0
Kobe 1 0
jordan 0 1可以看到james都是0，也就是因变量，然后Kobe是1，然后乔丹是1，也就是先拿科比作为1，然后保持科比不变，拿乔丹作为1，看和詹姆斯命中的影响，既然看残差符合正态分布，所以可以通过两个方式看qq图或正态性检测，先看第一个qqplot：
fit2.1<-lm(good~name,b)
summary(fit2.1)
qqPlot(fit2.1)先看看fit2.1的情况：
Call:
lm(formula = good ~ name, data = b)
Residuals:
Min 1Q Median 3Q Max
-2.0 -0.6 0.2 0.4 2.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.0000 0.5164 11.619 6.92e-08 ***
nameKobe -0.6000 0.7303 -0.822 0.427336
namejordan 3.8000 0.7303 5.203 0.000221 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.155 on 12 degrees of freedom
Multiple R-squared: 0.7806, Adjusted R-squared: 0.7441
F-statistic: 21.35 on 2 and 12 DF, p-value: 0.0001115
可以看到，namejordan那一列的p值是0.000221，说明保持科比不变的前提下二者有显著性差异，但是反而科比和詹姆斯没啥差异（p=0.427336）。然后看看qq图：
所有残差都在两个红虚线，95%置信区间内，如果不放心，我们再用第二种正态检测下，并画出图看看。
shapiro.test(rstudent(fit2.1))
第一个就是shapiro正态性检测，结果为
Shapiro-Wilk normality test
data: rstudent(fit2.1)
W = 0.97576, p-value = 0.9323
看到，有0.9323的概率是正态分布。然后看密度图：
plot(density(rstudent(fit2.1)))
比较符合正态分布。
本来想只说说单因素方差分析的，结果说了这么多。
下面继续刚才的话题，刚才乔丹加入了，我们看到三个人的投篮命中数：
name good
1 Kobe 5
2 Kobe 4
3 Kobe 6
4 Kobe 5
5 Kobe 7
6 James 7
7 James 5
8 James 6
9 James 4
10 James 8
11 jordan 10
12 jordan 9
13 jordan 10
14 jordan 10
15 jordan 10到现在为止，都是单因素方差分析，下面看看双因素分析；
詹姆斯和科比在看了这篇文章之后，觉得老流氓简直太厉害了，但是詹姆斯不服科比，我应该比你投篮准，但是因为这是你家后院，你家后院篮球场场地不好太晃眼，我们三个在每个人家的篮球场比比怎么样？
老流氓自然无所谓，科比倒是也觉得这样公平，那么把三个场地的每个人的投篮命中个数都统计下：
b1<-c(5,4,3,7,5,8,8,9,7,8,9,9,9,7,8)
b2<-c(6,6,7,6,7,8,7,6,7,8,9,10,10,8,9)
c<-cbind(b,b1,b2)
colnames(c)<-c('name','count.J','count.K','count.j')
library(reshape2)
d<-melt(c,variable_name = 'count')d数据框就是整形好的数据：
name count value
1 Kobe count.J 5
2 Kobe count.J 4
3 Kobe count.J 6
4 Kobe count.J 5
5 Kobe count.J 7
6 James count.J 7
7 James count.J 5
8 James count.J 6
9 James count.J 4
10 James count.J 8
11 jordan count.J 10
12 jordan count.J 9
13 jordan count.J 10
14 jordan count.J 10
15 jordan count.J 10
16 Kobe count.K 5
17 Kobe count.K 4
18 Kobe count.K 3
19 Kobe count.K 7
20 Kobe count.K 5
21 James count.K 8
22 James count.K 8
23 James count.K 9
24 James count.K 7
25 James count.K 8
26 jordan count.K 9
27 jordan count.K 9
28 jordan count.K 9
29 jordan count.K 7
30 jordan count.K 8
31 Kobe count.j 6
32 Kobe count.j 6
33 Kobe count.j 7
34 Kobe count.j 6
35 Kobe count.j 7
36 James count.j 8
37 James count.j 7
38 James count.j 6
39 James count.j 7
40 James count.j 8
41 jordan count.j 9
42 jordan count.j 10
43 jordan count.j 10
44 jordan count.j 8
45 jordan count.j 9count为场地
那么我们看到三个人在三块自己的主场上作战，在这里需要针对实际情况进行分析，在现有的情况下，我们看的是每个人在每块场地的投篮命中，也就是投篮命中和场地和人都有关系，所以应该看的是两个因素的
交互项
，也就是变量乘积：
fit3<-aov(value~name*count,data=d)
summary(fit3)看到fit3结果是
Df Sum Sq Mean Sq F value Pr(>F)
name 2 97.91 48.96 47.891 7.18e-11 ***
count 2 2.84 1.42 1.391 0.26181
name:count 4 18.76 4.69 4.587 0.00427 **
Residuals 36 36.80 1.02
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1不同姓名间的投篮命中确实有差别，同时，在不同场地的不同人之间也有显著差别，那么跟第一个分析一样，差别在哪里呢？
这里用HH包的交互图说明下：
library(HH)
interaction2wt(value~name*count,data=d)
图形如下：
先看左上图，横轴是三个人，纵轴是命中数量，颜色是三块场地，从这个图上看，交错比较明显，所以不同场地对所有人来说，差异不是很明显；和summary中的count p值为0.26相符；
看左下，横纵不变，颜色为三个人，可以明显看出，不管在那个场地，Jordan命中都比其他两人高；
看右上，横轴为三个场地，场地和命中差异不明显；
看右下，在任何的场地上，James表现比科比好，而在KOBE主场，这个差异最明显。
下面，看另一种情况，在前面我们假设的是场地和投篮命中有交互关系，那么如果三个人住在一起，在每个不同场地投篮的心态都一样，那么二者没有交互，我们的方差分析只看两个条件的差异：场地的差异和球员的差异：
fit4<-aov(value~name+count,data=d)
summary(fit4)结果：
Df Sum Sq Mean Sq F value Pr(>F)
name 2 97.91 48.96 35.248 1.49e-09 ***
count 2 2.84 1.42 1.024 0.368
Residuals 40 55.56 1.39
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1看到就是姓名之间命中率有差异，而场地之间无差异；
这次利用duncan.test进行不同姓名间比对：
duncan.test(fit4,'name',console=T)结果为
Study: fit4 ~ "name"
Duncan's new multiple range test
for value
Mean Square Error: 1.388889
name, means
value std r Min Max
James 7.066667 1.3345233 15 4 9
jordan 9.133333 0.9154754 15 7 10
Kobe 5.533333 1.2459458 15 3 7
alpha: 0.05 ; Df Error: 40
Critical Range
2 3
0.8697324 0.9144844
Means with the same letter are not significantly different.
Groups, Treatments and means
a jordan 9.133
b James 7.067
c Kobe 5.533 差异主要看后面最后三行，a/b/c，都是单个的值，则说明三个人都有差异，这次在我修改数值后，三个人都有了较大的差异，比如如果James和Kobe没差异的话，应该是bc，cb这样的表示。
2、多元方差分析
其实多元方差只是因变量多了一个，比如我在这个数据集后面加上篮板的参数（实在数据太多，随机生成吧），
set.seed(1000)
e<-cbind(d,rebound=abs(round(rnorm(45,10,3))))
perform<-cbind(e$value,e$rebound)最新的数据集是
name count value rebound
1 Kobe count.J 5 9
2 Kobe count.J 4 6
3 Kobe count.J 6 10
4 Kobe count.J 5 12
5 Kobe count.J 7 8
6 James count.J 7 9
7 James count.J 5 9
8 James count.J 6 12
9 James count.J 4 10
10 James count.J 8 6
11 jordan count.J 10 7
12 jordan count.J 9 8
13 jordan count.J 10 10
14 jordan count.J 10 10
15 jordan count.J 10 6
16 Kobe count.K 5 11
17 Kobe count.K 4 10
18 Kobe count.K 3 10
19 Kobe count.K 7 4
20 Kobe count.K 5 11
21 James count.K 8 18
22 James count.K 8 6
23 James count.K 9 13
24 James count.K 7 12
25 James count.K 8 8
26 jordan count.K 9 12
27 jordan count.K 9 5
28 jordan count.K 9 11
29 jordan count.K 7 12
30 jordan count.K 8 14
31 Kobe count.j 6 9
32 Kobe count.j 6 12
33 Kobe count.j 7 8
34 Kobe count.j 6 9
35 Kobe count.j 7 5
36 James count.j 8 11
37 James count.j 7 9
38 James count.j 6 13
39 James count.j 7 8
40 James count.j 8 6
41 jordan count.j 9 8
42 jordan count.j 10 14
43 jordan count.j 10 11
44 jordan count.j 8 10
45 jordan count.j 9 9上数据集的value就是投篮good，很郁闷reshape2包的value.name不好用了。。。
然后如果我们假设人和场地对篮球运动员的表现有显著地影响，那么：
fit5<-manova(perform~name*count,data=e)
summary(fit5)得出
Df Pillai approx F num Df den Df Pr(>F)
name 2 0.76600 11.1735 4 72 4.17e-07 ***
count 2 0.14389 1.3954 4 72 0.24419
name:count 4 0.40204 2.2644 8 72 0.03213 *
Residuals 36
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1可以看出，还是姓名对表现有很强的影响，以及人和场地的交互影响。
下面细分来看：
summary.aov(fit5)结果：
Response 1 :
Df Sum Sq Mean Sq F value Pr(>F)
name 2 97.911 48.956 47.8913 7.178e-11 ***
count 2 2.844 1.422 1.3913 0.261805
name:count 4 18.756 4.689 4.5870 0.004266 **
Residuals 36 36.800 1.022
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Response 2 :
Df Sum Sq Mean Sq F value Pr(>F)
name 2 9.644 4.8222 0.5865 0.5615
count 2 21.111 10.5556 1.2838 0.2894
name:count 4 14.222 3.5556 0.4324 0.7842
Residuals 36 296.000 8.2222 这个就和一元方差类似，response1就是投篮，2是篮板，可以看出，投篮跟之前一样有很大差异，但是篮板和人、场地以及交互都没有差异。
3.协方差分析
如果考虑一下这种情况：詹姆斯对乔丹不服，说你进联盟这么长时间了，所以经验丰富，科比经验第二丰富，如果我们剔除掉训练因素，看看究竟谁天赋更高？这就好像几个学生考试，小明考了100分，小美考了90分，但是小明天天学习24小时，小美天天学习1小时，如果我们比较小明和小美的本身天赋（及去掉努力的因素），这时候就该用协方差分析。
简单按照年龄大小排序，乔丹50，科比40，詹姆斯30，ok，建立数据框
f<-cbind(d,age=rep(c(rep(40,5),rep(30,5),rep(50,5)),3))建好的数据框如下：
name count value age
1 Kobe count.J 5 40
2 Kobe count.J 4 40
3 Kobe count.J 6 40
4 Kobe count.J 5 40
5 Kobe count.J 7 40
6 James count.J 7 30
7 James count.J 5 30
8 James count.J 6 30
9 James count.J 4 30
10 James count.J 8 30
11 jordan count.J 10 50
12 jordan count.J 9 50
13 jordan count.J 10 50
14 jordan count.J 10 50
15 jordan count.J 10 50
16 Kobe count.K 5 40
17 Kobe count.K 4 40
18 Kobe count.K 3 40
19 Kobe count.K 7 40
20 Kobe count.K 5 40
21 James count.K 8 30
22 James count.K 8 30
23 James count.K 9 30
24 James count.K 7 30
25 James count.K 8 30
26 jordan count.K 9 50
27 jordan count.K 9 50
28 jordan count.K 9 50
29 jordan count.K 7 50
30 jordan count.K 8 50
31 Kobe count.j 6 40
32 Kobe count.j 6 40
33 Kobe count.j 7 40
34 Kobe count.j 6 40
35 Kobe count.j 7 40
36 James count.j 8 30
37 James count.j 7 30
38 James count.j 6 30
39 James count.j 7 30
40 James count.j 8 30
41 jordan count.j 9 50
42 jordan count.j 10 50
43 jordan count.j 10 50
44 jordan count.j 8 50
45 jordan count.j 9 50先看看年龄的影响
fit6<-aov(value~age+name+count,f)
summary(fit6)结果为
Df Sum Sq Mean Sq F value Pr(>F)
age 1 32.03 32.03 23.064 2.22e-05 ***
name 1 65.88 65.88 47.432 2.69e-08 ***
count 2 2.84 1.42 1.024 0.368
Residuals 40 55.56 1.39
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1可以看到count依然没差异，但是age的差异比较大，先看看包含age的三个人的投篮命中均值：
aggregate(f$value,by=list(f$name),mean)得出
Group.1 x
1 James 7.066667
2 Kobe 5.533333
3 jordan 9.133333ok，暂时乔丹领先，那么看看调整后的均值：
library(effects)
effect('name',fit6)结果为：
name effect
name
James Kobe jordan
8.100000 5.533333 8.100000 看到詹姆斯和乔丹的命中相当，也就是说假以时日，詹姆斯的投篮能赶上乔丹。
OK，这篇文章看起来好长，虽然作为勒布朗的球迷给他贴了金，但是三者的原理有些相似但略有不同，据说线性回归和方差分析都是广义回归模型的特例，以后可能回头再看看有什么联系。