<-c(60,142,195,80,242,220,190,25,198,38,236,95)
test1<-c(76,152,243,82,240,220,205,38,243,44,190,100) test2
5 秩转换的非参数检验
5.1 配对样本比较的Wilcoxon符号秩检验
英文名字:wilcoxon signed-rank test,即符号秩检验。
使用课本例8-1的数据,自己手动摘录:
两列数据,和配对t检验的数据结果完全一样。
简单看一下数据情况:
boxplot(test1,test2)
进行秩和检验:
wilcox.test(test1,test2,paired = T,alternative = "two.sided",exact = F)
##
## Wilcoxon signed rank test with continuity correction
##
## data: test1 and test2
## V = 11.5, p-value = 0.06175
## alternative hypothesis: true location shift is not equal to 0
结果和课本一致!
wilcox.test
在R4.4.0及以后不再支持在公式形式中使用paired
参数,所以如果你的R版本在4.4.0及以后的,在公式形式中使用paired
参数会得到以下报错:cannot use 'paired' in formula method
。逆天更新!!
5.2 两独立样本比较的Wilcoxon秩和检验
英文名字:wilcoxon rank sum test,即秩和检验。
根据课本上的说法,wilcoxon rank sum test 和 mann-whitney U test 并不完全是一个东西。
和两样本t检验的数据格式完全一样!
使用课本例8-3的数据,自己手动摘录。
<-c(2.78,3.23,4.20,4.87,5.12,6.21,7.18,8.05,8.56,9.60)
RD1<-c(3.23,3.50,4.04,4.15,4.28,4.34,4.47,4.64,4.75,4.82,4.95,5.10) RD2
进行两独立样本比较的Wilcoxon符号秩检验:
wilcox.test(RD1,RD2,paired = F, correct = F, exact = F)
##
## Wilcoxon rank sum test
##
## data: RD1 and RD2
## W = 86.5, p-value = 0.08049
## alternative hypothesis: true location shift is not equal to 0
结果取单侧检验,还是和课本一致!
5.3 完全随机设计多个样本比较的 Kruskal-Wallis H 检验
5.3.1 多样本比较的kruskal-wallis H检验
使用课本例8-5的数据,手动摘录:
rm(list = ls())
<- c(32.5,35.5,40.5,46,49,16,20.5,22.5,29,36,6.5,
death_rate 9.0,12.5,18,24)
<- rep(c("Drug_A","drug_B","drug_C"),each=5)
drug <- data.frame(death_rate,drug)
mydata
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 : chr "Drug_A" "Drug_A" "Drug_A" "Drug_A" ...
数据一共2列,第1列是死亡率,第2列是药物(3种)。
简单看下数据:
boxplot(death_rate ~ drug, data = mydata)
进行 Kruskal-Wallis H 检验:
kruskal.test(death_rate ~ drug, data = mydata)
##
## Kruskal-Wallis rank sum test
##
## data: death_rate by drug
## Kruskal-Wallis chi-squared = 9.74, df = 2, p-value = 0.007673
算出来结果和课本一致!
使用课本例8-6的数据,手动摘录:
<- data.frame(days=c(2,2,2,3,4,4,4,5,7,7,
data8_6 5,5,6,6,6,7,8,10,12,
3,5,6,6,6,7,7,9,10,11,11),
type=c(rep("9D",10),rep("11C",9),rep("DSC",11)))
head(data8_6)
## days type
## 1 2 9D
## 2 2 9D
## 3 2 9D
## 4 3 9D
## 5 4 9D
## 6 4 9D
进行 Kruskal-Wallis H 检验:
kruskal.test(days ~ type, data = data8_6)
##
## Kruskal-Wallis rank sum test
##
## data: days by type
## Kruskal-Wallis chi-squared = 9.9405, df = 2, p-value = 0.006941
结果和课本一致(课本中Hc=9.97)。
5.3.2 kruskal-Wallis H检验后的多重比较
例8-8使用的是例8-6的数据。
课本上是使用的 Nemenyi检验,我们使用非参数检验的全能R包:PMCMRplus
实现。
library(PMCMRplus)
# 首先要把分类变量因子化
$type <- factor(data8_6$type) data8_6
下面就可以使用 Nemenyi检验了。
# 也可以把kwh检验的结果作为输入
<- kwAllPairsNemenyiTest(days ~ type, data = data8_6)
res ## Warning in kwAllPairsNemenyiTest.default(c(2, 2, 2, 3, 4, 4, 4, 5, 7, 7, : Ties
## are present, p-values are not corrected.
summary(res)
## q value Pr(>|q|)
## 9D - 11C == 0 3.628 0.027794 *
## DSC - 11C == 0 0.177 0.991411
## DSC - 9D == 0 3.998 0.013057 *
得到的数值和课本有差别,但是结论是一样的。
5.4 随记区组设计多个样本比较的Friedman M检验
5.4.1 多个相关样本比较的Friedman M检验
使用课本例8-9的数据:
<- foreign::read.spss("datasets/例08-09.sav", to.data.frame = T)
df
str(df)
## 'data.frame': 8 obs. of 4 variables:
## $ a: num 8.4 11.6 9.4 9.8 8.3 8.6 8.9 7.8
## $ b: num 9.6 12.7 9.1 8.7 8 9.8 9 8.2
## $ c: num 9.8 11.8 10.4 9.9 8.6 9.6 10.6 8.5
## $ d: num 11.7 12 9.8 12 8.6 10.6 11.4 10.8
## - attr(*, "codepage")= int 65001
数据一共4列,分别是4中不同频率下的反应率。
简单看下数据:
boxplot(df$a,df$b,df$c,df$d)
进行 Friedman M 检验前先把数据格式转换一下:
<- as.matrix(df) # 变成矩阵 M
进行 Friedman M 检验:
friedman.test(M)
##
## Friedman rank sum test
##
## data: M
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
结果和课本一致!
5.4.2 多个相关样本两两比较的q检验
P126页,多个相关样本两两比较的q检验。课本上说的这个q检验,应该是quade test。
接下来就是使用R语言实现quade-test
。但是自带的quade.test()
函数不能进行两两比较,还是要借助第三方包。
# 准备数据,也是用的课本例8-9的数据
<- matrix(
df 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检验看看:
friedman.test(df)
##
## Friedman rank sum test
##
## data: df
## Friedman chi-squared = 15.152, df = 3, p-value = 0.001691
接下来进行quade检验:
library(PMCMRplus)
quadeAllPairsTest(df, dist = "Normal")
## A B C
## B 0.2200 - -
## C 0.0017 0.0644 -
## D 1.7e-07 7.7e-05 0.0860
当然也可以有更加详细的结果:
<- quadeAllPairsTest(df,dist = "Normal")
res toTidy(res)
## group1 group2 statistic p.value alternative
## 1 B A 1.226488 2.200150e-01 two.sided
## 2 C A 3.526154 1.686568e-03 two.sided
## 3 C B 2.299666 6.440153e-02 two.sided
## 4 D A 5.549859 1.715396e-07 two.sided
## 5 D B 4.323371 7.683144e-05 two.sided
## 6 D C 2.023706 8.600089e-02 two.sided
## method distribution p.adjust.method
## 1 Quade's testwith standard-normal approximation z holm
## 2 Quade's testwith standard-normal approximation z holm
## 3 Quade's testwith standard-normal approximation z holm
## 4 Quade's testwith standard-normal approximation z holm
## 5 Quade's testwith standard-normal approximation z holm
## 6 Quade's testwith standard-normal approximation z holm
这个结果和课本上也不是完全一样,不过不影响结果。还有很多其他的方法可以选择,除了这个quade检验,还可以用Nemenyi等检验方法。