R语言统计与绘图:组间差异的多重比较

前面我们一起学习了如何进行多组独立样本的差异比较,包括参数和非参数的检验方法。但备择假设为组间的差异不全相等,如拒绝原假设的话,通常需要进行组间差异多重比较,即两两比较

针对不同的设计和数据特征,我们可以选择合适的方法,如图1所示。

图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 LSD.test()函数设置Bonferroni法调整P值结果
图3 Groups and Range图示

图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 pairwise.t.test()函数运行结果(bonferroni检验)

结果如图4所示,可以看出,结果与图2、图3一致。

Bonferroni检验比较保守,在组数较多时需要比较的次数较多,校正后的检验水准过小,有时很难达到。

Holm是一种常用的修正Bonferroni过保守的方法,将p.adj设置为”holm”即可。

图5的结果显示,相对于校正之前的Bonferroni法,holm的结果没那么严格,A药物组与B药物组、B药物组与C药物组、C药物组与D药物组间的P值接近检验水准。

图5 pairwise.t.test()函数运行结果(holm法)

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 ##显示归类结果
图6 LSD-t检验结果

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 Dunnett-t检验结果

如图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)
图8 Tukey检验结果1
图9 效应之差的置信区间

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# 调整检验水准
图10 K-W H检验后的多重比较1

kruskalmc()函数以逻辑判断的方式表示是否有显著性差异,如图10所示,检验水准为0.05时,仅高剂量组和对照组间差异有统计学意义。采用更严格的检验水准后,各组间不存在显著性差异。

我们还可以用PMCMR包中的posthoc.kruskal.nemenyi.test()函数实现,代码和结果如下:

library(PMCMR)
posthoc.kruskal.nemenyi.test(data4$TNF,data4$group)
图11 K-W H检验后的多重比较2

如图11所示,结果与图10检验水准为0.05时的结果一致,高剂量组和对照组间差异有统计学意义。

2.2.2 多个相关样本

如果上一例中的设计变为随机区组设计,即如图12中矩阵数据(行为区组,列为处理组)。

图12 Friedman M检验的数据

如果在Friedman M检验后需要两两比较,可以用PMCMR包中的posthoc.friedman.nemenyi.test()函数实现,数据和代码如下:

library(PMCMR)
posthoc.friedman.nemenyi.test(data4.1)

小结

组间差异多重比较的方法较多,同一种方法在R中又有一到多种实现方式,因此,实际应用中我们需要根据特定的目的、设计和数据的特征来选择最适合的分析方法。那么,看完上面的介绍,对于手上的数据你有一些思路了吗?

参考来源:

  1. 孙振球,徐勇勇.医学统计学:第4版[M].北京:人民卫生出版社.2014.
  2. 方积乾,徐勇勇,陈峰,等.卫生统计学:第7版[M].北京:人民卫生出版社.2012.
  3. [美]Robert I. Kabacoff著.R语言实战(第2版) [M].王小宁等译.北京:人民邮电出出版社.2016.
统计与绘图

什么是区组随机化?

2020-6-23 20:04:01

统计与绘图

评价干预措施效果的指标:RR,RRR,ARR,NNT

2020-6-23 20:14:24

0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
购物车
优惠劵
今日签到
有新私信 私信列表
搜索