R语言之方差分析

Posted 2020-09-03

一、单因素方差分析

单因素方差分析只有一个分组变量，因此数据看起来像一个多列的数据框，如
 Grass Heath Arable
1      3     6     19
2      4     7      3
3      3     8      8
4      5     8      8
5      6     9      9
6     12    11     11
7     21    12     12
8      4    11     11
9      5    NA      9
10     4    NA     NA
11     7    NA     NA
12     8    NA     NA

方差分析的基本命令是aoc()，并且需要使用公式语法，并且数据结构也要为预测变量+因子的形式，因此对于前面的数据，如果直接使用将出现错误，我们使用stack()命令将其转化，转化后的数据形式如下：
  values    ind
1       3  Grass
2       4  Grass
3       3  Grass
4       5  Grass
5       6  Grass
6      12  Grass
7      21  Grass
8       4  Grass
9       5  Grass
10      4  Grass
11      7  Grass
12      8  Grass
13      6  Heath
14      7  Heath
15      8  Heath
16      8  Heath
17      9  Heath
18     11  Heath
19     12  Heath
20     11  Heath
21     NA  Heath
22     NA  Heath
23     NA  Heath
24     NA  Heath
25     19 Arable
26      3 Arable
27      8 Arable
28      8 Arable
29      9 Arable
30     11 Arable
31     12 Arable
32     11 Arable
33      9 Arable
34     NA Arable
35     NA Arable
36     NA Arable
原数据中带有NA项，如果想将其移除，可使用na.omit()。

将数据转化成预测变量+因子的形式之后，我们使用aov()命令进行方差分析，如：
> aov(count~site,data=bfs)
Call:
   aov(formula = count ~ site, data = bfs)

Terms:
                    site Residuals
Sum of Squares   55.3678  467.6667
Deg. of Freedom        2        26

Residual standard error: 4.24113
Estimated effects may be unbalanced

我们对结果使用summary()命令，将其呈现为经典的方差分析表格，结果如下：
> summary(aov(count~site,data=bfs))
            Df Sum Sq Mean Sq F value Pr(>F)
site         2   55.4   27.68   1.539  0.233
Residuals   26  467.7   17.99    
可以看出，上面的结果包含了F值及显著性。

方差分析经常涉及事后检验，我们可以使用TukeyHSD()命令来执行Tukey Honest检验。

二、多因素方差分析

多因素方差分析需要考虑因素间的交互作用，比单因素复杂一些，如有数据如下：
  height    plant water
1       9 vulgaris    lo
2      11 vulgaris    lo
3       6 vulgaris    lo
4      14 vulgaris   mid
5      17 vulgaris   mid
6      19 vulgaris   mid
7      28 vulgaris    hi
8      31 vulgaris    hi
9      32 vulgaris    hi
10      7   sativa    lo
11      6   sativa    lo
12      5   sativa    lo
13     14   sativa   mid
14     17   sativa   mid
15     15   sativa   mid
16     44   sativa    hi
17     38   sativa    hi
18     37   sativa    hi

aov()命令可按照如下设置：
（1）> aov(height~plant+water,data = pw)
（2）> aov(height~plant*water,data = pw)
（3）> aov(height~plant+water+plant:water,data = pw)

（1）表示分析~右侧的因素，（2）表示分析~右侧的因素和它们之间的交互作用，（3）表示分析~右侧的因素和指定两个因素的交互作用，当只有两个因素时，（2）和（3）是等价的。
对于这三个方差分析结果，我们可以进一步使用anova()命令对他们进行对比，如：

> aov1<-aov(height~plant+water,data = pw)
> aov2<-aov(height~plant*water,data = pw)
> aov3<-aov(height~plant+water+plant:water,data = pw)
> anova(aov1,aov2,aov3)
Analysis of Variance Table

Model 1: height ~ plant + water
Model 2: height ~ plant * water
Model 3: height ~ plant + water + plant:water
  Res.Df     RSS Df Sum of Sq      F   Pr(>F)   
1     14 199.111                                
2     12  69.333  2    129.78 11.231 0.001783 **
3     12  69.333  0      0.00                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以看出，aov1和aov2有明显不同，这也说明因素间的交互作用对结果产生影响，而aov2和aov3结果是完全相同的。

由于存在交互作用，因此在进行事后检验的时候，结果会非常多，我们可以通过TukeyHSD()命令的which选项指定想查看的因素，如：
> TukeyHSD(aov1,which = "water")
我们只查看water因素的两两对比。


在分析交互作用时，通常要做交互作用图，可以通过interaction.plot()命令实现，基础命令为：
interaction.plot(x.factor,trace.factor,response...)
x.factor判定交互作用如何划分
trace.factor为x.factor中划分类别的方式，实际上就是组合x.factor与trace.factor来显示交互作用。
response为指定的因变量。

三、更复杂的方差分析模型