R语言非参数检验多重比较

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。

之前介绍了多个样本均数的多重比较，今天说说kruskal-Wallis H检验后的多重比较，Friedman M检验后的多重比较。

也是和课本对照着来，孙振球，徐勇勇《医学统计学》第四版。

非参数检验后的多重比较，我们也是用这个宝藏R包：PMCMRplus。

kruskal-Wallis H检验及多重比较

使用课本例8-5的数据。

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 ...
1
2
3
4
5
6
7
8
9
10
11

进行kruskal-Wallis H 检验：

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
1
2
3
4
5
6
7

多重检验，课本上用的是Nemenyi检验，我们通过多重比较的全能R包PMCMRplus实现。

library(PMCMRplus)
1

也是提供两种输入方式，直接提供kruskal-Wallis H检验的结果，或者formula形式，都可以。

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

除此之外，还提供了kwAllPairsConoverTest()和kwAllPairsDunnTest()方法。

Friedman M检验及多重比较

使用课本本例8-9的数据，这个方式适用于随机区组设计资料的多样本比较。

进行Friedman M检验需要矩阵形式的数据(这个是R语言里为数不多的不支持formula形式的统计检验函数之一)，可以自己输入，也可以直接读取spss格式数据，然后变成矩阵即可。

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

进行Friedman M检验：

fit <- friedman.test(df)
fit
## 
## 	Friedman rank sum test
## 
## data:  df
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
1
2
3
4
5
6
7

使用q检验（quade test）进行多重比较:

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
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

可以看到非常简单，函数名称清晰易懂，结果也是非常直观，直接给出了两两比较的P值和统计量。

除此之外，还提供了非常多其他方法:

• frdAllPairsConoverTest()
• frdAllPairsExactTest()
• frdAllPairsConoverTest()
• frdAllPairsMillerTest()
• frdAllPairsSiegelTest()

大家选择自己需要的即可。

本文首发于公众号：医学和生信笔记，完美观看体验请至公众号查看本文。

医学和生信笔记，专注R语言在临床医学中的使用，R语言数据分析和可视化。