前面我们一起学习了如何进行多组独立样本的差异比较，包括参数和非参数的检验方法。但备择假设为组间的差异不全相等，如拒绝原假设的话，通常需要进行组间差异多重比较，即两两比较。
针对不同的设计和数据特征，我们可以选择合适的方法，如图1所示。
下面，我们将接着方差分析和非参数检验的例子进一步介绍在R中实现组间差异的多重比较。
目 录
- 1. 参数检验
- 1.1 前文回顾
- 1.2 均数间的两两比较
- 1.2.1 Bonferroni检验及校正
- 1.2.2 LSD-t检验
- 1.2.3 Dunnett-t检验
- 1.2.4 Tukey HSD检验
- 1.2.5 Scheffe检验
- 1.2.6 SNK-q检验
- 2. 非参数检验
- 2.1 前文回顾
- 2.2 组间差异两两比较
- 2.2.1 多个独立样本
- 2.2.2 多个相关样本
- 小结
1. 参数检验
1.1 前文回顾
因为我们要做的是多组独立样本参数检验后的均数间的两两比较，我们先简单回顾完全随机设计的方差分析这个例子，详情可见《方差分析》。
某临床实验中心想测试五种药物降低胆固醇的疗效，将50名参与者随机分成5组，每组分别接受一种药物(drugA、drugB、drugC、drugD和drugE)，数据资料如表1所示，变量包括使用的药物类型(trt，分组变量)和胆固醇的降低值(response，反应变量)。
前文进行方差分析，结果显示P＜0.05，拒绝原假设，不能认为五个组的均数是相等。
1.2 均数间的两两比较
因此，我们需要进一步做均数间的两两比较。均数间的两两比较检验方法有的需要直接使用数据，有的则使用模型方差。
下面为大家介绍几种常见的检验方法。
1.2.1 Bonferroni检验及校正
Bonferroni检验用途最广，几乎可用于任何多重比较的情形，可用agricolae包中LSD.test()函数或基础包中pairwise.t.test()函数实现。
其基本思想是采用α’=α/n作为新的检验水准，其中，n为两两比较次数, α为累积I类错误的概率。
我们首先使用agricolae包中LSD.test()函数，注意将P值的调整方法p.adj设置为bonferroni。代码和结果如下：
library(agricolae) # 加载包
bon <- LSD.test(aov,"trt", p.adj="bonferroni")
bon$groups # 显示归类结果
plot(bon)
图2和图3显示了比较的结果，两组间有相同的字母表示差异不显著，两组间字母不同表示差异显著。
例如本例，E药物组只有单独字母a，与其他4组皆不同，代表E药物组分别与其他4个药物组的疗效有显著差异。而D药物组与C药物组有相同的b，与B药物组和A药物组无相同字母，表示D药物组与C药物组差异不显著，D药物组分别与B药物组和A药物组间有显著差异。
我们还可以用基础包中的pairwise.t.test()函数做Bonferroni检验。代码和结果如下：
attach(a2)
pairwise.t.test(response,trt, p.adj = "bonferroni")
detach()
结果如图4所示，可以看出，结果与图2、图3一致。
Bonferroni检验比较保守，在组数较多时需要比较的次数较多，校正后的检验水准过小，有时很难达到。
Holm是一种常用的修正Bonferroni过保守的方法，将p.adj设置为”holm”即可。
图5的结果显示，相对于校正之前的Bonferroni法，holm的结果没那么严格，A药物组与B药物组、B药物组与C药物组、C药物组与D药物组间的P值接近检验水准。
1.2.2 LSD-t检验
LSD-t检验，也叫最小显著性差异，由agricolae包中的LSD.test()函数实现。依据t分布，是t检验的简单变形。
LSD检验适用于在专业上有特殊意义的样本均数间的比较，是在设计之初，就已明确要比较某几个组均数间是否有差异。LSD法侧重于减小II类错误，但有增大I类错误(假阳性)的可能。代码和结果如下：
library(agricolae)
L <- LSD.test(aov,"trt", p.adj = "none")
L$groups ##显示归类结果
LSD-t检验不像Bonferroni检验那么保守，所以结果有所不同。如图6所示，各药物组对应的指示字母均不同，表明各药物组的疗效有显著差异。
1.2.3 Dunnett-t检验
Dunnett-t检验由multcomp包中glht()函数实现，适用于k-1个试验组与一个对照组均数差异的多重比较。这种比较也是在设计阶段就确定了。
library(multcomp)
D<-glht(aov, linfct = mcp(trt = 'Dunnett'), alternative = 'two.side')
summary(D)
如图7所示，模型默认第一组为对照组，并依次对比了其他组和它的差异，发现B组与A组的疗效差异不显著，而C、D、E药物组分别与A药物组的疗效有显著差异。
1.2.4 Tukey HSD检验
Tukey HSD检验由基础包中的Tukey HSD()函数或multcomp包中的glht()函数实现，用于各组样本均数间的全面比较，适用于组间例数相等的情况。代码和结果如下：
TukeyHSD(aov)
plot(TukeyHSD(aov))
###或者##
library(multcomp)
T<-glht(aov, linfct = mcp(trt = "Tukey"))
summary(T)
plot(T)
TukeyHSD()和glht()函数的运行结果一致，图8和图9结果显示，除了B-A、C-B、D-C外，其他各组间的差异均具有统计学意义。
1.2.5 Scheffe检验
Scheffe检验由agricolae包中的scheffe.test()函数实现。各组样本数相等或不等均可以使用，但是以各组样本数不相等时使用较多。Scheffe也是通过指示字母的相同与否判定二者间差异是否有统计学意义，不显示两两比较的P值。代码如下：
library(agricolae)
scheffe.test(aov,"trt", console=TRUE)
1.2.6 SNK-q检验
SNK-q检验由agricolae包中的SNK.test()函数实现，适用于多个样本均数两两之间的全面比较，与LSD-t检验相似，可能存在假阳性。代码如下：
library(agricolae)
S<-SNK.test(aov,"trt")#console=T表示输出结果
2. 非参数检验
2.1 前文回顾
在数据不满足正态分布等条件的情况下，组间差异性比较通常会选择非参数检验。
同样的，多组样本的非参数检验也存在无法进行组间两两比较的问题。在前文《非参数检验》中，我们分析了这样一个案例：
15只家兔按体重随机分为三组，分别注入无菌生理盐水，0.05μg/g和0.10μg/g剂量的DON进行实验处理，实验期满后测定关节冲洗液中肿瘤坏死因子(TNF-α)的水平(μg/L)，数据如表2，比较这三组家兔关节冲洗液TNF-α测定结果差异是否具有统计学意义。
2.2 组间差异两两比较
我们采用多组独立样本数据的K-W平均秩检验，结果显示拒绝原假设，可以认为这三组家兔关节冲洗液TNF-α测定结果差异有统计学意义。
2.2.1 多个独立样本
多组独立样本非参数检验后组间差异的两两比较可以用pgirmess包中的kruskalmc()函数或PMCMR包中的posthoc.kruskal.nemenyi.test()函数实现。
我们先看kruskalmc()函数。代码和结果如下：
library(pgirmess)
kruskalmc(data4$TNF,data4$group)
kruskalmc(data4$TNF,data4$group, probs=0.001) # 调整检验水准
kruskalmc()函数以逻辑判断的方式表示是否有显著性差异，如图10所示，检验水准为0.05时，仅高剂量组和对照组间差异有统计学意义。采用更严格的检验水准后，各组间不存在显著性差异。
我们还可以用PMCMR包中的posthoc.kruskal.nemenyi.test()函数实现，代码和结果如下：
library(PMCMR)
posthoc.kruskal.nemenyi.test(data4$TNF,data4$group)
如图11所示，结果与图10检验水准为0.05时的结果一致，高剂量组和对照组间差异有统计学意义。
2.2.2 多个相关样本
如果上一例中的设计变为随机区组设计，即如图12中矩阵数据(行为区组，列为处理组)。
如果在Friedman M检验后需要两两比较，可以用PMCMR包中的posthoc.friedman.nemenyi.test()函数实现，数据和代码如下：
library(PMCMR)
posthoc.friedman.nemenyi.test(data4.1)
小结
组间差异多重比较的方法较多，同一种方法在R中又有一到多种实现方式，因此，实际应用中我们需要根据特定的目的、设计和数据的特征来选择最适合的分析方法。那么，看完上面的介绍，对于手上的数据你有一些思路了吗？
参考来源：
-
孙振球,徐勇勇.医学统计学:第4版[M].北京:人民卫生出版社.2014. -
方积乾,徐勇勇,陈峰,等.卫生统计学:第7版[M].北京:人民卫生出版社.2012. -
[美]Robert I. Kabacoff著．R语言实战(第2版) [M].王小宁等译．北京:人民邮电出出版社.2016.
