菜鸟学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???不能插入附件嘛?这是其中的一张图。

R语言 比较两个数据集的差异 r语言比较两组数据大小_ide


具体用什么哪个方法,还是根据实际情况吧。选几个用代码试试,我拿自己的实验数据举栗子了。

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

R语言 比较两个数据集的差异 r语言比较两组数据大小_R语言_02


扯了不少废话,上代码!!!

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不一样。
写到这里,难道我论文中的结果,要修改嘛???

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