首页 > 代码库 > R语言之方差分析

R语言之方差分析

一、单因素方差分析

单因素方差分析只有一个分组变量,因此数据看起来像一个多列的数据框,如
 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=http://www.mamicode.com/bfs)
Call:
   aov(formula = count ~ site, data = http://www.mamicode.com/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=http://www.mamicode.com/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 = http://www.mamicode.com/pw)
(2)> aov(height~plant*water,data = http://www.mamicode.com/pw)
(3)> aov(height~plant+water+plant:water,data = http://www.mamicode.com/pw)

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

> aov1<-aov(height~plant+water,data = http://www.mamicode.com/pw)
> aov2<-aov(height~plant*water,data = http://www.mamicode.com/pw)
> aov3<-aov(height~plant+water+plant:water,data = http://www.mamicode.com/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为指定的因变量。

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

R语言之方差分析