R语言系列第四期：②R语言多组样本方差分析与KW检验
本文首发于 ”百味科研芝士“ 微信公众号，转载请注明：百味科研芝士，Focus科研人的百味需求
R语言系列四的第二个部分是对多组连续性数据的处理，分组往往是三组或者三组以上，当然两组数据也可以利用方差分析，但是两组数据还是建议使用t检验。同样多组数据的比较也分为参数法和非参数法，包括这个部分介绍的重点参数法方差分析，以及非参数方法kruskal—Wallis检验。
A. 单因素方差分析
我们先从一个简单的单因素开始，单因素可以理解为各个组间的差别只有一个因素，而我们研究的就是这个因素的影响。以此类推，两因素方差分析，就是研究组间两个因素(及其交互作用)的影响，这个之后再谈。
我们利用奥特曼（Altman，1991）收集的“红细胞叶酸盐”数据。首先我们来看一下数据情况：
> library(ISwR)
> attach(red.cell.folate)
> red.cell.folate
folate ventilation
1 243 N2O+O2,24h
...
8 392 N2O+O2,24h
9 206 N2O+O2,op
...
17 309 N2O+O2,op
18 241 O2,24h
...
22 328 O2,24h
> summary(red.cell.folate)
folate ventilation
Min. :206.0 N2O+O2,24h:8
1st Qu.:249.5 N2O+O2,op :9
Median :274.0 O2,24h :5
Mean :283.2
3rd Qu.:305.5
Max. :392.0
#Tips: 可以看出来这个数据集的数据是测量值和分组情况分别放在两个变量里，同时数据是分成三组的，它们分别是“24小时内的O2和N2O含量”“手术中O2和N2O含量”“24小时内O2含量”。另外，可以通过summary()函数可以看出，数值向量与属性变量的汇总格式是不一样的。
下面我们就利用anova()函数将方差分析表取出来：
> anova(lm(folate~ventilation))
Analysis of Variance Table
Response: folate
Df Sum Sq Mean Sq F value Pr(>F)
ventilation 2 15516 7757.9 3.7113 0.04359 *
Residuals 19 39716 2090.3
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：在公式里有一个lm()函数，它是一个线性模型函数（linear model），我们会在相关和回归里重点介绍它，这里只需要记住写法。同样，和t检验和wilcoxon检验一样这里，这里有“~”，而“~”之前的变量是数值变量，之后是分组变量。
在统计教材里，平方和一般都被分为“组间”和“组内”。在R语言中组间方差的平方和利用分组属性变量的名字（ventilation）来称呼，而组内方差直接利用Residual来标注。我们可以看一下统计量F value的值是3.7113，pr代表p值结果为0.04359*<0.05,因此我们可以认为三组的均值不全相等。方差分析表能被用在很多的统计模型上，但统计模型的格式与比较各组数据差异这个特定问题并不相同，用起来也很方便。
当然，有些差别还是要点出来，因为有的时候因为数据的格式问题，可能导致错误的结果。比如juul数据集的例子。这个数据中的变量tanner是个数值向量，而不是属性向量。对于列出的表格没有任何影响，但是在做方差分析时就会出现严重错误。
> attach(juul)
> tail(juul)
> summary(juul)
age menarche sex
Min. : 0.170 Min. :1.000 Min. :1.000
1st Qu.: 9.053 1st Qu.:1.000 1st Qu.:1.000
Median :12.560 Median :1.000 Median :2.000
Mean :15.095 Mean :1.476 Mean :1.534
3rd Qu.:16.855 3rd Qu.:2.000 3rd Qu.:2.000
Max. :83.000 Max. :2.000 Max. :2.000
NA's :5 NA's :635 NA's :5
igf1 tanner testvol
Min. : 25.0 Min. :1.00 Min. : 1.000
1st Qu.:202.2 1st Qu.:1.00 1st Qu.: 1.000
Median :313.5 Median :2.00 Median : 3.000
Mean :340.2 Mean :2.64 Mean : 7.896
3rd Qu.:462.8 3rd Qu.:5.00 3rd Qu.:15.000
Max. :915.0 Max. :5.00 Max. :30.000
NA's :321 NA's :240 NA's :859
> anova(lm(igf1~tanner))
Analysis of Variance Table
Response: igf1
Df Sum Sq Mean Sq F value Pr(>F)
tanner 1 10985605 10985605 686.07 < 2.2e-16 ***
Residuals 790 12649728 16012
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：这里并不是真正的方差分析，而是一个对于分组编号的线性回归！因为变量tanner的自由度为1，而真正方差分析的自由度是组数-1，应该是5-1=4.
我们可以做如下的更改：
> juul$tanner<-factor(tanner,labels=c("I","II","III","IV","V"))
> detach(juul)
> attach(juul)
> summary(tanner)
I II III IV V NA's
515 103 72 81 328 240
> anova(lm(igf1~tanner))
Analysis of Variance Table
Response: igf1
Df Sum Sq Mean Sq F value Pr(>F)
tanner 4 12696217 3174054 228.35 < 2.2e-16 ***
Residuals 787 10939116 13900
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Tips：把连续数据转换成因子变量就可以解决这个问题，R就会明白做的不是回归，而是方差分析。另外，因为更改了数据内的值，需要重新绑定数据集juul。
我们可以通过df的值，来查看我们的计算是否正确，这个例子就是告诉我们如果要做方差分析，分组的变量必须是属性变量或者因子。
#Tips：anova()不能处理数据分组盛放的情况，必须有一个变量是存放分组的。
B. 多重比较
前面的F检验提示我们组间有差异，那么问题马上变成差异到底在哪里。这时候就需要进行组与组之间的两两比较了。
如果我们比较所有的组别，应该进行多重检验的修正。进行多次检验，会增加其中出现一个显著结果的概率；也就是说，这个p值会变得夸张。一个常用的调整方法是bonferroni修正法。bonferroni修正法进行任意两组间的比较。
函数pairwise.t.test()能够计算所有的两组比较。他能够针对多重检验做调整：
> pairwise.t.test(folate,ventilation,p.adj="bonferroni")
Pairwise comparisons using t tests with pooled SD
data: folate and ventilation
N2O+O2,24h N2O+O2,op
N2O+O2,op 0.042 -
O2,24h 0.464 1.000
P value adjustment method: bonferroni
输出结果是一个成对比较的p值表。我们通过结果可以发现只有N2O+O2,24h和N2O+O2,op这两组之间的p值是<0.05的。说明这两组之间差异有显著性意义，其他组之间可认为无差别。
C. 图像显示
当然，现在有很多种方法来展示分组数据。这里我们展示一个精妙复杂的图形，其中原始数据用条形图画出来，然后再叠加上均值与标准误。
> xbar=tapply(folate,ventilation,mean)
> s=tapply(folate,ventilation,sd)
> n<-tapply(folate,ventilation,length)
> sem=s/sqrt(n)
> stripchart(folate~ventilation,method="jitter",jitter=0.05,pch=16,vert=T)
> arrows(1:3,xbar+sem,1:3,xbar-sem,angle=90,code=3,length=0.1)
> lines(1:3,xbar,pch=4,type="b",cex=2)
#Tips：前四列是基本信息的赋值，stripchart()函数是绘制图形的函数，里面的参数method=“jitter”是我们之前使用过的，代表数据在绘制方向上垂直波动，避免重叠，jitter=0.05代表波动的尺度大小。pch=16表示小圆点，vert=T表示把方向改成垂直方向，而不是水平方向。arrows()函数能在图形上添加箭头。我们稍微灵活地利用箭头的头部可调整这一特性，在两端都加上一个交叉图像。前四个参数表示端点；参数angle指的是箭头和剑柄之间的角度，这里设置为90度；参数length指的是箭头的长度。最后code=3表示两端都有箭头。lines()函数不解释了。每次使用基本上需要更改的地方只有tapply和stripchart()的前两个参数和arrows()和lines()中的1:3,3改成自己的组数就可以了。
D. Bartlett检验
不知道大家还记不记得t检验里的方差齐性检验，var.test(),这里因为方差分析也是参数检验，所以也需要数据满足方差齐性，但在方差分析里有专门的检验方差齐性检验方法bartlett检验。
> bartlett.test(folate~ventilation)
Bartlett test of homogeneity of variances
data: folate by ventilation
Bartlett's K-squared = 2.0951, df = 2, p-value = 0.3508
这里可见p>0.05，所以可以认为三组的方差是齐的。可以进行方差分析。
E. Kruskal—Walis检验
它是方差分析的非参数版本。当数据不满足正态分布，或者数据类型不适合做方差分析的时候可以考虑KW检验，它同样比较的是数值的秩次而不是数值本身，这里不做过多的赘述。
> kruskal.test(folate~ventilation)
Kruskal-Wallis rank sum test
data: folate by ventilation
Kruskal-Wallis chi-squared = 4.1852, df = 2, p-value = 0.1234
这里的结果就与参数检验方差分析的结果相悖，这里显示的是3组数据间无显著性差异。足以看出参数检验更能检验出阳性结果，检验效能高于非参数检验。
F. 双因素方差分析
单因素方差分析处理的是依据单因素分类的数据。我们也能够分析依据不同的准则交叉分类的数据。双因素方差分析需要将数据放在一个向量里，以及与其平行的两个分类属性。我们以一个使用依那普利拉之后的心率数据（Altman，1991）作为例子。数据集heart.rate中的数据是这样的形式：
> attach(heart.rate)
> heart.rate
hr subj time
1 96 1 0
...
9 100 9 0
10 92 1 30
...
18 106 9 30
19 86 1 60
...
27 104 9 60
28 92 1 120
...
36 102 9 120
#Tips: 这里虽然分组变量subj和time变量都是数值变量，但是这个数据集在数据框中已经被定义成为因子。这个是原始数据就完成的，所以我们不需要再转化成因子，但是使用任何数据做方差分析前一定要确保数据的分组变量是正确的格式。
数据一旦被定义好了，双因素方差分析就可以简单地如下计算：
> anova(lm(hr~subj+time))
Analysis of Variance Table
Response: hr
Df Sum Sq Mean Sq F value Pr(>F)
subj 8 8966.6 1120.82 90.6391 4.863e-16 ***
time 3 151.0 50.32 4.0696 0.01802 *
Residuals 24 296.8 12.37
通过数据的结果显示，两个因素的p值都是小于0.05的，因此两个因素对结果的影响都具有显著性意义。
#Tips：在模型方程中交换subj和time，除了方差分析表中两行的顺序有变化，产生一模一样的分析结果（如果是不平衡设计的话，属性的顺序会有很大影响）。当然这里是没有交互效应的结果，我们这里只介绍到单独因素。
以上就是对方差分析在R语言中的简单利用，当然还有很多高深的方法，感兴趣的朋友可以网络上搜索。
t检验和方差分析都是对连续型数据的分析方法，当遇到离散数据或者分类数据时就需要改变方法了，之后我们就会介绍分类数据的处理方法。我们下期再见。
参考资料：
1. 《R语言统计入门（第二版）》 人民邮电出版社 Peter Dalgaard著
2. 《R语言初学者指南》 人民邮电出版社 Brian Dennis著
3. Vicky的小笔记本《blooming for you》by Vicky
