[R语言统计]秩转换的非参数检验

来源：互联网发布：matlab 矩阵转置编辑：程序博客网时间：2024/05/01 04:56

非参数检验(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 and new
V = 11.5, p-value = 0.06175
alternative hypothesis: true location shift is not equal 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 is not equal 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 and si
W = 86.5, p-value = 0.04318
alternative hypothesis: true location shift is greater 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 39名吸烟工人和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类似

白细胞

支气管扩张

肺水肿

肺癌

病毒性呼吸道感染

+++

# 源代码：
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 in 1: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 - 11C 0.95110000
DSC1 - 9D  0.01088889

例8-9 8名受试对象在相同实验条件下分别接受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 sum test

data: matrix89
Friedman chi-squared = 15.1519, df = 3, p-value = 0.001691

例8-10 对例8-9资料作四个样本间的两两比较。

还没有查到关于Friedman检验的两两比较方法，以后补充。

参考资料：

关于Nemenyi检验的方法是参照丁香园的，原贴地址如下，帖子的作者有把写过《医学统计学及SAS应用》之R语言实现这本电子书，下载地址如下。

Nemenyi检验方法：http://www.dxy.cn/bbs/thread/26928446#26928446

《医学统计学及SAS应用》之R语言实现：http://d.dxy.cn/detail/6143147

0 0