非参数检验(non-parametric test)是相对于参数检验(parametric test)而言的。如果总体分布为已知的数学形式,用参数检验,反之用非参数检验。当总体分布不能由已知的数学形式表达,没有总体参数时,就无法用参数检验,两个或多个正态总体方差不等,也不能用t检验或F检验的参数检验。对于不满足参数检验条件的数据,一是进行变量变换,使其满足参数检验条件,另外就是用非参数检验。
非参检验对总体分布不作严格假定,又称任意分布检验(distribution-free test),《医学统计学》(第三版,孙振球)书中采用的是秩转换的非参数检验,即将数值变量从小到大排列,再计算检验统计量。
目录
配对Wilcoxon符号秩检验
独立样本Wilcoxon检验(Mann-Whitney U检验)
Kruskal-Wallis H检验(完全随机)
多个独立样本两两比较的Nemenyi检验法
Friedman M检验法
多个相关样本两两比较的q检验法(不会)
配对Wilcoxon符号秩检验
例8-1 对12份血清分别用原方法(检测时间20分钟)和新方法(检测时间10分钟)测谷-丙转氨酶,结果见表8-1的(2)、(3)栏。问两法所得结果有无差别?
# 源代码:
old <-c(60,142,195,80,242,220,190,25,198,38,236,95)
new <-c(76,152,243,82,240,220,205,38,243,44,190,100)
wilcox.test(old,new,paired =TRUE)# 运行结果
>wilcox.test(old,new,paired =TRUE) #配对wilcoxon检验
Wilcoxon signed rank test with continuity correction
data: old andnew
V =11.5, p-value =0.06175
alternative hypothesis:true location shift isnotequal to 0
Warning messages:
1:In wilcox.test.default(old, new, paired =TRUE) :
cannot compute exact p-value with ties
2:In wilcox.test.default(old, new, paired =TRUE) :
cannot compute exact p-value with zeroes
结果分析: p-value=0.06175,因此不能认为两法测谷丙转氨酶有差别。
例8-2 已知某地正常人尿氟含量的中位数为45.30μmol/L 。今在该地某厂随机抽取12名工人,测得尿氟含量见表8-2第(1)栏。问该厂工人的尿氟含量是否高于当地正常人的尿氟含量?
# 源代码:
data82 <-c(44.21,45.30,46.39,49.47,51.05,53.16,53.26,54.37,57.16,67.37,71.05,87.37)
wilcox.test(data82-45.30)# 运行结果:
>wilcox.test(data82,mu=45.30)
Wilcoxon signed rank test with continuity correction
data: data82
V =65, p-value =0.005099
alternative hypothesis:true location isnotequal to 45.3
Warning message:
In wilcox.test.default(data82, mu =45.3) :
cannot compute exact p-value with zeroes
结果分析:p值小于0.01,可以视为该工厂的工人尿氟含量高于正常人。
例8-3 对10例肺癌病人和12例矽肺0期工人用X光片测量肺门横径右侧距RD值(cm),结果见表8-5。问肺癌病人的RD值是否高于矽肺0期工人的RD值?
# 源代码:
lc <-c(2.78,3.23,4.2,4.87,5.12,6.21,7.18,8.05,8.56,9.6)
si <-c(3.23,3.5,4.04,4.15,4.28,4.34,4.47,4.64,4.75,4.82,4.95,5.1)
wilcox.test(lc,si)#运行结果:
>wilcox.test(lc,si,alternative="greater",exact=F)
Wilcoxon rank sum test with continuity correction
data: lc andsi
W =86.5, p-value =0.04318
alternative hypothesis:true location shift isgreater than 0
在R中,wilcox.test()函数可以用来做Wilcoxon秩和检验,也可以用于做Mann-Whitney U检验。当参数为单个样本,或者是两个样本相减,或者是两个参数,paired=F时,是Wilcoxon秩和检验。当paired = FALSE(独立样本)时,就是Mann-Whitney U检验,在这个题目中,用的就是Mann-Whitney U检验,虽然结果中W=86.5与书中的T=141.5不一样,但本质上是一样的,换算如下:W1=141.5-10*(10+1)/2=86.5;W2=111.5-12*(12+1)/2=33.5
例8-4 名吸烟工人和40名不吸烟工人的碳氧血红蛋白HbCO(%)含量见表8-6。问吸烟工人的HbCO(%)含量是否高于不吸烟工人的HbCO(%)含量?
答:本题书中用的是wilcox检验,其实用Ridit分析更合适一些,下面分别用这两种方法进行检验:
# 1. Wilcox.test(或Kruskal检验)
smoke <- c(1,8,16,10,4)
no.smoke <-c(2,23,11,4,0)
rank.c <- c(1:5)
group1 <- rep(rank.c,smoke)
group2 <- rep(rank.c,no.smoke)
data84 <- c(group1,group2)
group.f <-factor(c(rep(1,length(group1)),rep(2,length(group2))))
wilcox.test(data84~group.f)
# 或者进行kruskal.test检验
kruskal.test(data84~group.f)
# 运行结果:
> wilcox.test(data84~group.f)
Wilcoxon rank sum test with continuity correction
data: data84 by group.f
W = 1137, p-value = 0.0002181
alternative hypothesis: true location shift is not equal to 0
Warning message:
In wilcox.test.default(x = c(1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, :
cannot compute exact p-value with ties
> # 或者进行kruskal.test检验
> kruskal.test(data84~group.f)
Kruskal-Wallis rank sum test
data: data84 by group.f
Kruskal-Wallis chi-squared = 13.707, df = 1, p-value = 0.0002137
# 2. 此题用Ridit检验更合适一些。
data84 <- matrix(c(1,8,16,10,4,2,23,11,4,0),nrow=5)
dimnames(data84) <- list(cone=c("very low","low","median","a bit high","high"),
worker=c("smoker","no-smoker"))
library(Ridit)
ridit(data84,2)
chisq.test(data84)
# 运行结果:
> library(Ridit)
> ridit(data84,2)
Ridit Analysis:
Group Label Mean Ridit
----- ----- ----------
1 smoker 0.6159
2 no-smoker 0.387
Reference: Total of all groups
chi-squared = 13.707, df = 1, p-value = 0.0002137
> chisq.test(data84)
Pearsons Chi-squared test
data: data84
X-squared = 15.0785, df = 4, p-value = 0.004541
三、Kurskal-Wallis检验
Kurskal-Wallis检验是Wilcoxon方法(其实是Mann-Whitney检验)用于多于两个样本的时候的升级版。当对两个样本进行比较的时候,Kurskal-Wallis检验与Mann-Whitney检验是等价的。
例8-5 用三种药物杀灭钉螺,每批用200只活钉螺,用药后清点每批钉螺的死亡数、再计算死亡率(%),结果见表8-9。问三种药物杀灭钉螺的效果有无差别?
# 源代码:
drug <-rep(c("甲药","乙药","丙药"),each=5)
data <-c(32.5,35.5,40.5,46,49,16,20.5,22.5,29,36,6.5,9.0,12.5,18,24)
data85 <-data.frame(drug,data)
kruskal.test(data85$data~data85$drug)# 运行结果:
>kruskal.test(data85$data~data85$drug)
Kruskal-Wallis rank sum test
data: data85$data by data85$drug
Kruskal-Wallis chi-squared =9.74, df =2, p-value =0.007673
结果分析:Kruskal-Wallis的统计量是H值。
例8-6 比较小白鼠接种三种不同菌型伤寒杆菌9D、11C和DSC1后存活日数,结果见表8-10。问小白鼠接种三种不同菌型伤寒杆菌的存活日数有无差别?
# 源代码:
mice <-as.factor(c(rep("9D",10),rep("11C",9),rep("DSC1",11)))
data86 <-c(2,2,2,3,4,4,4,5,7,7,5,5,6,6,6,7,8,10,12,3,5,6,6,6,7,7,9,10,11,11)
data86 <-data.frame(mice,data)
kruskal.test(data86$data~data86$mice)# 运行结果:
>kruskal.test(data86$data~data86$mice)
Kruskal-Wallis rank sum test
data: data86$data by data86$mice
Kruskal-Wallis chi-squared =9.9405, df =2, p-value =0.006941
例8-7 四种疾病患者痰液内嗜酸性白细胞的检查结果见表8-11。问四种疾病患者痰液内的嗜酸性白细胞有无差别?注:这道例题与《医学统计学及SAS应用》(上海交通大学)的9.11类似
白细胞 | 支气管扩张 | 肺水肿 | 肺癌 | 病毒性呼吸道感染 |
- | 0 | 3 | 5 | 3 |
+ | 2 | 5 | 7 | 5 |
++ | 9 | 5 | 3 | 3 |
+++ | 6 | 2 | 2 | 0 |
# 源代码:
x1 <-c(0,2,9,6)
x2 <-c(3,5,5,2)
x3 <-c(5,7,3,2)
x4 <-c(3,5,3,0)
freq <-function(x){
count <-c()
for(i in1:4){
count1 <-c(rep(i,x[i]))
count <-append(count,count1)
}
return(count)}
data87 <-c(freq(x1),freq(x2),freq(x3),freq(x4))
group<-c(rep(1,sum(x1)),rep(2,sum(x2)),rep(3,sum(x3)),rep(4,sum(x4)))
kruskal.test(data87~group)# 运行结果:
>kruskal.test(data87~group)
Kruskal-Wallis rank sum test
data: data87 by group
Kruskal-Wallis chi-squared =15.5058, df =3, p-value =0.001432
例8-8 对例8-6资料(表8-10)作三个样本间的两两比较,Nemenyi检验。
# 源代码:
# 需要安装下列包:
# install.packages("pgirmess")
# install.packages("coin")
# install.packages("multcomp")
# library(pgirmess)
# library(coin)
# library(multcomp) # 安装并加载要用到的包
mice <-as.factor(c(rep("9D",10),rep("11C",9),rep("DSC1",11)))
data <-c(2,2,2,3,4,4,4,5,7,7,5,5,6,6,6,7,8,10,12,3,5,6,6,6,7,7,9,10,11,11)
data88 <-data.frame(mice,data)
kruskal.test(data~mice,data=data88)
kruskalmc(data~mice, data=data88, probs=0.05) # 使用kruskalmc函数做两两比较,但此方法不能给出具体的值
# 下面构建函数计算具体的p值
mult <-oneway_test(data ~mice, data =data88,
ytrafo =function(data) trafo(data, numeric_trafo =rank),
xtrafo =function(data) trafo(data, factor_trafo =function(x)
model.matrix(~x -1) %*%t(contrMat(table(x), "Tukey"))),
teststat ="max", distribution =approximate(B =90000))
pvalue(mult, method ="single-step") #计算具体的p值,这段代码是复制丁香园的[1]# 运行结果:
>kruskal.test(data~mice,data=data88)
Kruskal-Wallis rank sum test
data: data by mice
Kruskal-Wallis chi-squared =9.9405, df =2, p-value =0.006941
>kruskalmc(data~mice, data=data88, probs=0.05) # 使用kruskalmc函数做两两比较,但此方法不能给出具体的值
Multiple comparison test after Kruskal-Wallis
p.value:0.05
Comparisons
obs.dif critical.dif difference
11C-9D 10.3777778 9.683378 TRUE
11C-DSC1 0.4949495 9.472590 FALSE
9D-DSC1 10.8727273 9.208410 TRUE
>pvalue(mult, method ="single-step") #计算具体的p值
9D-11C 0.01870000
DSC1 -11C0.95110000
DSC1 -9D 0.01088889
例8-9 名受试对象在相同实验条件下分别接受4种不同频率声音的刺激,他们的反应率(%)资料见表8-12。问4种频率声音刺激的反应率是否有差别?
# 源代码:
freA <-c(8.40,11.60,9.40,9.80,8.30,8.60,8.90,7.80)
freB <-c(9.60,12.70,9.10,8.70,8.00,9.80,9.00,8.20)
freC <-c(9.80,11.80,10.40,9.90,8.60,9.60,10.60,8.50)
freD <-c(11.70,12.00,9.80,12.00,8.60,10.60,11.40,10.80)
matrix89 <-matrix(c(freA,freB,freC,freD),nrow=8,dimnames=list(no=1:8,c("A","B","C","D")))
friedman.test(matrix89)运行结果:
>friedman.test(matrix89)
Friedman rank sumtest
data: matrix89
Friedman chi-squared =15.1519, df =3, p-value=0.001691
例8-10 对例8-9资料作四个样本间的两两比较。
还没有查到关于Friedman检验的两两比较方法,以后补充。