之前介绍了多个样本均数的多重比较，今天说说kruskal-Wallis H检验后的多重比较，Friedman M检验后的多重比较。
也是和课本对照着来，孙振球，徐勇勇《医学统计学》第四版。本书电子版已上传到qq群中，大家加群即可免费获取。
非参数检验后的多重比较，我们也是用这个宝藏R包：PMCMRplus。
kruskal-Wallis H检验及多重比较
使用课本例8-5的数据。
代码语言：javascript复制rm(list = ls())
death_rate <- c(32.5,35.5,40.5,46,49,16,20.5,22.5,29,36,6.5,9.0,12.5,18,24)
drug <- rep(c("Drug_A","drug_B","drug_C"),each=5)
mydata <- data.frame(death_rate,drug)
# 分类变量因子化
mydata$drug <- factor(mydata$drug)
str(mydata)
## 'data.frame': 15 obs. of 2 variables:
## $ death_rate: num 32.5 35.5 40.5 46 49 16 20.5 22.5 29 36 ...
## $ drug : Factor w/ 3 levels "Drug_A","drug_B",..: 1 1 1 1 1 2 2 2 2 2 ...
进行kruskal-Wallis H 检验：
代码语言：javascript复制fit <- kruskal.test(death_rate ~ drug, data = mydata)
fit
##
## Kruskal-Wallis rank sum test
##
## data: death_rate by drug
## Kruskal-Wallis chi-squared = 9.74, df = 2, p-value = 0.007673
多重检验，课本上用的是Nemenyi检验，我们通过多重比较的全能R包PMCMRplus实现。
library(PMCMRplus)
也是提供两种输入方式，直接提供kruskal-Wallis H检验的结果，或者formula形式，都可以。
代码语言：javascript复制res <- kwAllPairsNemenyiTest(fit)
res <- kwAllPairsNemenyiTest(death_rate ~ drug, data = mydata)
summary(res)
##
## Pairwise comparisons using Tukey-Kramer-Nemenyi all-pairs test with Tukey-Dist approximation
## data: death_rate by drug
## alternative hypothesis: two.sided
## P value adjustment method: single-step
## H0
## q value Pr(>|q|)
## drug_B - Drug_A == 0 2.5 0.1805089
## drug_C - Drug_A == 0 4.4 0.0052932 **
## drug_C - drug_B == 0 1.9 0.3710425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
除此之外，还提供了kwAllPairsConoverTest()和kwAllPairsDunnTest()方法。
Friedman M检验及多重比较
使用课本本例8-9的数据，这个方式适用于随机区组设计资料的多样本比较。
进行Friedman M检验需要矩阵形式的数据(这个是R语言里为数不多的不支持formula形式的统计检验函数之一)，可以自己输入，也可以直接读取spss格式数据，然后变成矩阵即可。
代码语言：javascript复制df <- matrix(
c(8.4, 11.6, 9.4, 9.8, 8.3, 8.6, 8.9, 7.8,
9.6, 12.7, 9.1, 8.7, 8, 9.8, 9, 8.2,
9.8, 11.8, 10.4, 9.9, 8.6, 9.6, 10.6, 8.5,
11.7, 12, 9.8, 12, 8.6, 10.6, 11.4, 10.8
),
byrow = F, nrow = 8,
dimnames = list(1:8,LETTERS[1:4])
)
print(df)
## A B C D
## 1 8.4 9.6 9.8 11.7
## 2 11.6 12.7 11.8 12.0
## 3 9.4 9.1 10.4 9.8
## 4 9.8 8.7 9.9 12.0
## 5 8.3 8.0 8.6 8.6
## 6 8.6 9.8 9.6 10.6
## 7 8.9 9.0 10.6 11.4
## 8 7.8 8.2 8.5 10.8
进行Friedman M检验：
代码语言：javascript复制fit <- friedman.test(df)
fit
##
## Friedman rank sum test
##
## data: df
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
使用q检验（quade test）进行多重比较:
代码语言：javascript复制res <- quadeAllPairsTest(df)
summary(res)
##
## Pairwise comparisons using Quade's test with TDist approximation
## data: df
## P value adjustment method: holm
## H0
## t value Pr(>|t|)
## B - A == 0 1.126 0.27307148
## C - A == 0 3.236 0.01583866 *
## D - A == 0 5.093 0.00028885 ***
## C - B == 0 2.110 0.14099045
## D - B == 0 3.967 0.00351141 **
## D - C == 0 1.857 0.15475978
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
可以看到非常简单，函数名称清晰易懂，结果也是非常直观，直接给出了两两比较的P值和统计量。
除此之外，还提供了非常多其他方法:
frdAllPairsConoverTest()frdAllPairsExactTest()frdAllPairsMillerTest()frdAllPairsSiegelTest()
大家选择自己需要的即可。
