菜鸟学R语言（组间多重比较）

菜鸟学R语言（组间多重比较）

经过方差分析可以说明各总体均值间的差异是否显著，即只能说明均值不全相等，但不能具体说明哪几个均值之间有显著差异。此时多重比较就派上用场了，在科研中也是比较常用的方法。
具体的理论知识不再多讲，上代码之前先了解一下多重比较的类别。

1. LSD检验（最小显著差数检验法）：这也是我最常用的方法，基本上就是T检验的简单变形，T检验是对两组，而这个可以对多组间的均数做检验；
2. Dunnett检验：适用于多个试验组与一个对照组的比较，多对一；
3. Turkey检验：适用于组数大于6以上（不确定）；
4. SNK法（Student-Newman-Keuls）：最为流行的方法，广泛使用，但只告诉有无差异,不提供精确P值。
5. Duncan法：没用过，不过结果应该和LSD差不多；
6. Scheffe检验：适用于任何比较。
7. Bonferroni法：当比较次数不多时（小于5组），Bonferroni法的效果较好。

进行两两比较的时候，假如是验证性研究用Bonferroni或LSD比较好，Bonferroni法适用于小于5组的比较。假如是探索性研究且各组人数相同用TurKey法较好，其他的用scheffe较好。

我解释的不清楚，有兴趣的话可以看看下面这个附件。
WHAT???不能插入附件嘛？这是其中的一张图。

具体用什么哪个方法，还是根据实际情况吧。

选几个用代码试试，我拿自己的实验数据举栗子了。
就不展示进行正态分布和方差齐性检验的部分了。下图是我论文中的一幅图，每个柱子上边的的abc字母不同代表差异显著。可以认为是5个变量（5个深度），4个水平（4种植被，用1 2 3 4表示），或者反过来说也可以吧。。因变量是同位素值。比较各个深度下，不同的植被同位素值的组间差异。第一个和最后一个深度，方差分析显示差异不显著，其实也就没有必要做多重比较了，但我也标上字母了，只不过都一样。

扯了不少废话，上代码！！！

LSD检验

library(agricolae)
qiudexun <- aov(nb1~type1, data = soil)
out1 <- LSD.test(qiudexun, "type1", p.adj = "bonferroni" )
print(out1$groups)
qiudexun <- aov(nb2~type2, data = soil)
out1 <- LSD.test(qiudexun, "type2", p.adj = "bonferroni" )
print(out1$groups)
qiudexun <- aov(nb3~type3, data = soil)
out1 <- LSD.test(qiudexun, "type3", p.adj = "bonferroni" )
print(out1$groups)
qiudexun <- aov(nb4~type4, data = soil)
out1 <- LSD.test(qiudexun, "type4", p.adj = "bonferroni" )
print(out1$groups)
qiudexun <- aov(nb5~type5, data = soil)
out1 <- LSD.test(qiudexun, "type5", p.adj = "bonferroni" )
print(out1$groups)

结果如下：

> library(agricolae)
> qiudexun <- aov(nb1~type1, data = soil)
> out1 <- LSD.test(qiudexun, "type1", p.adj = "bonferroni" )
> print(out1$groups)
        nb1 groups
1 -7.683238      a
4 -7.932266      a
3 -8.271141      a
2 -8.523507      a
> qiudexun <- aov(nb2~type2, data = soil)
> out1 <- LSD.test(qiudexun, "type2", p.adj = "bonferroni" )
> print(out1$groups)
         nb2 groups
4  -9.775052      a
1 -10.147676     ab
3 -10.161522     ab
2 -11.173642      b
> qiudexun <- aov(nb3~type3, data = soil)
> out1 <- LSD.test(qiudexun, "type3", p.adj = "bonferroni" )
> print(out1$groups)
        nb3 groups
4 -10.24492      a
3 -10.53397      a
1 -11.53943     ab
2 -12.50955      b
> qiudexun <- aov(nb4~type4, data = soil)
> out1 <- LSD.test(qiudexun, "type4", p.adj = "bonferroni" )
> print(out1$groups)
         nb4 groups
3  -9.768635      a
2 -10.134467     ab
1 -10.336852     ab
4 -11.266095      b
> qiudexun <- aov(nb5~type5, data = soil)
> out1 <- LSD.test(qiudexun, "type5", p.adj = "bonferroni" )
> print(out1$groups)
         nb5 groups
3  -9.633797      a
2  -9.936403      a
1 -10.376280      a
4 -11.256957      a

代码应该可以简化，不过重要的方法和结果，写论文谁在乎你的过程啊（我胡说的）。
用的方法是LSD，但是注意p.adj = “bonferroni”，这是不是就等同于Bonferroni法？如果去掉这行代码呢？用nb2~type2检验一下。

> qiudexun <- aov(nb2~type2, data = soil)
> out1 <- LSD.test(qiudexun, "type2")
> print(out1$groups)
         nb2 groups
4  -9.775052      a
1 -10.147676      a
3 -10.161522      a
2 -11.173642      b

可以看到结果是不同的！分别是a, a, a, b和a, ab, ab, b！但哪个对呢？我也不清楚。。我做的时候加上了p.adj = “bonferroni”

Dunnett检验

这个适用于有对照组的数据，多个组分别于一个组进行比较。
其实我的数据也有对照组，PR、CK、AG、Cropland分别是刺槐、柠条、撂荒草地、农地。农地就是对照组，下面试试：

试了半天出错了，凸(艹皿艹 )，略过！

Tukey检验

上面提到Tukey检验适用于组数大于6以上，这个其实并不是判断依据，如果各组例数相等，Tukey法也很不错！而且是所有组之间的比较，我感觉在不考虑适用条件的情况下，Tukey检验可以代替Dunnett检验，只不过输出结果是P值的形式。

> qiudexun <- aov(nb1 ~ as.factor(type1), data = soil)
> tuk <- TukeyHSD(qiudexun)
> tuk
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = nb1 ~ as.factor(type1), data = soil)

$`as.factor(type1)`
         diff       lwr      upr     p adj
2-1 -0.840269 -4.288832 2.608294 0.9126657
3-1 -0.587903 -4.036466 2.860660 0.9673972
4-1 -0.249028 -3.697591 3.199535 0.9973443
3-2  0.252366 -3.196197 3.700929 0.9972372
4-2  0.591241 -2.857322 4.039804 0.9668702
4-3  0.338875 -3.109688 3.787438 0.9933945

> qiudexun <- aov(nb2 ~ as.factor(type2), data = soil)
> tuk <- TukeyHSD(qiudexun)
> tuk
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = nb2 ~ as.factor(type2), data = soil)

$`as.factor(type2)`
         diff         lwr         upr     p adj
2-1 -1.025966 -1.98726176 -0.06467024 0.0344268
3-1 -0.013846 -0.97514176  0.94744976 0.9999737
4-1  0.372624 -0.58867176  1.33391976 0.6893429
3-2  1.012120  0.05082424  1.97341576 0.0373150
4-2  1.398590  0.43729424  2.35988576 0.0036800
4-3  0.386470 -0.57482576  1.34776576 0.6651247

以nb2 ~ as.factor(type2)为例，注意这里要加as.factor，这是因为自变量用数字1234表示的，R可能识别不出来这是因子。看P值，小于0.05的有2-1,3-2,4-2，转换为字母表示（平均值由大到小为4, 1, 3, 2）：a, a, a, b。这与LSD方法结果一样，与bonferroni不一样。
写到这里，难道我论文中的结果，要修改嘛？？？

一言难尽，睡觉了，剩下的几种方法下次再补充。