但是，一旦数据不满足一定的分布条件，也就是不满足方差分析的条件时，我们就需要更换为非参数的方法，这里我们测试了一个非参数检验的方法示例，先执行Kruskal-Wallis检验比较整体差异，再执行Behrens-Fisher的非参数多重比较查看两两差异。
但是Behrens-Fisher 检验得到的是两组间的 p 值，对于显著性“abc”的标注，通过 p 值再结合中位数水平需要手动判断，大体方法和前面讲的类似。
但是，npmc包的最后更新版本（1.0.7）无法正确运行在R 3.0以上的版本。HK Zhang 在Rstudio中做了调试和编译，对npmc.R的原代码做了一些细微的更改以支持最新的mvtnorm包。npmc包最后编译运行成功。（注：R>4.0 目前也不行了，3.8可行）
最新原程序共享在GitHub:
https://github.com/HK-Zhang/Rice/tree/master/npmc
也可以通过下面链接直接下载这个包
http://pan.baidu.com/s/1nuqdXcX
还是用21.1的测试数据。
library(reshape2)
library(npmc)
library(ggplot2)
data <- read.table("test_data.txt",header=T,sep="\t")
我们先来测试一组数据，选取group1和value1的数据来测试看下。
#Kruskal-Wallis 检验，整体差异
dat <- data[c("group1", "value1")]
names(dat) <- c('class', 'var')
dat$class <- factor(dat$class)
fit <- kruskal.test(var~class, dat)
p_value <- fit$p.value
#获取 Kruskal-Wallis 检验 p 值
stat_kruskal <- NULL
if (p_value < 0.001) {
sig <- '***' } else if(p_value >= 0.001 & p_value < 0.01) {
sig <- '**' } else if(p_value >= 0.01 & p_value < 0.05) {
sig <- '*' } else {
sig <- '' }
stat_kruskal <- rbind(stat_kruskal, c("group1","value1",p_value, sig))
#Behrens-Fisher 检验（npmc 包），非参数多重比较
BF_test <- npmc(dat, df = 2, alpha = 0.05)
BF <- BF_test$test$BF
BF结果如下图所示，可以看出是group1内每个小组和小组之间value1的统计检验结果，pvalue在p.value.2s中，只不过每个小组用的1,2,3对应的index顺序。
info <- BF_test$info
info里面包含index和class的信息。
BF$group <- paste("group1", "value1", sep = '/')
所以，我们现在把BF中的index换成相应的组名。
for (i in 1:nrow(BF)){
x <- strsplit(as.character(BF$cmp[i]),"-")[[1]]
BF[i,'class1'] <-rownames(info)[which(info$group.index == x[1])]
BF[i,'class2'] <-rownames(info)[which(info$group.index == x[2])]
}
stat_BF <- NULL
stat_BF <- rbind(stat_BF, BF[c(12:14, 2:11)])
#这样子，我们就获得了每个组之间的统计检验显著性结果。
但是，如果我们想要标记abc，就需要自己手动来标记了，首先需要排序中位数，然后按21.1的方法来手动标记abc。
#下面，我们就计算中位数，然后按从高到低来排序。
dat_new <- cbind(aggregate(dat$var, by = list(dat$class), FUN = median),
aggregate(dat$var, by = list(dat$class), FUN = sd))
dat_new$group <- "group1"
dat_new$var <- "value1"
dat_new <- dat_new[c(5, 6, 1,2, 4)]
names(dat_new) <- c('group', 'var', 'class', 'median',"sd")
dat_new <- dat_new[order(dat_new$median, decreasing = TRUE), ]
class <- as.vector(dat_new$class)
下面，我们就来手动标记abc，一般做法如下（前面也讲过）：
（1）首先根据均值大小，将各组由高往低排排序，均值最高的组标注为“a”；
（2）将均值最高的组与第二高的组相比，若差异显著，则第二组标注为“b”；若不显著，继续比较其与均值第三高的组的差异；
（3）若均值最高的组与第二高的组不显著，均值第二高的组与第三高的组显著，则第二高的组就同样标注为“a”，第三高的组标注为“b”；若均值最高的组与第二高的组不显著、均值第二高的组与第三高的组不显著，但均值最高的组与第三高的组显著，则第二高的组就标注为“ab”，第三高的组标注为“b”；
（4）然后以标注为“b”的组的均值为标准，以此类推，继续循环往后比较，直到最小均值的组被标记，且比较完毕为止。
m = 1; n = 2; l = 1
dat_new[m,'sig'] <- letters[l] #最高的赋值为a
while(n<length(class)){
for(n in (m+1):length(class)) {
#从stat_BF中获取class[m]和class[n]的pvalue，但是首先需要判断一下顺序，index小的放前面
tmp1=info$group.index[which(info$class.level == class[m])]
tmp2=info$group.index[which(info$class.level == class[n])]
if(tmp1<tmp2){
index1=class[m]
index2=class[n]
}else{
index1=class[n]
index2=class[m]
}
tmp_p= stat_BF[intersect(which(stat_BF$class1==index1),which(stat_BF$class2==index2)),]$p.value.2s
#如果p小于显著性值，则letters依次望下排。
if(tmp_p<0.05){
dat_new[n,'sig'] <- letters[l+1]
if (n - m != 1) {
for (x in (m+1):(n-1)) {
tmp1=info$group.index[which(info$class.level == class[x])]
tmp2=info$group.index[which(info$class.level == class[n])]
if(tmp1<tmp2){
index1=class[x]
index2=class[n]
}else{
index1=class[n]
index2=class[x]
}
tmp_p= stat_BF[intersect(which(stat_BF$class1==index1),which(stat_BF$class2==index2)),]$p.value.2s
if(tmp_p<0.05){
dat_new[x,'sig'] <- letters[l]
}else{
dat_new[x,'sig'] <- paste(letters[l], letters[l+1], sep = '')
}
}
}
l = l + 1
m = n
break
}
}
dat_new[(m:n),'sig'] <- letters[l]
}
ggplot(data = dat_new, aes(x = class, y = median)) +
geom_col(aes(fill = class), color = NA, show.legend = FALSE) +
geom_errorbar(aes(ymin = median - sd, ymax = median + sd), width = 0.2) +
geom_text(aes(label = sig, y = median + sd + 0.3*(median+sd)),size=5) +
labs(x = '', y = '')+
theme(text = element_text(size = 20))
下面，像21.1一样，我们写个双层循环，把group1和group2分别和value进行显著性检验。
groups<-c('group1', 'group2')
values<-c('value1', 'value2', 'value3', 'value4')
stat_kruskal <- NULL
stat_BF <- NULL
data_list <- NULL
for(i in groups)
{
for (j in values)
{
cat("i:", i,"\n")
cat("j:", j,"\n")
#Kruskal-Wallis 检验，整体差异
dat <- data[c(i,j)]
names(dat) <- c('class', 'var')
dat$class <- factor(dat$class)
fit <- kruskal.test(var~class, dat)
p_value <- fit$p.value
#获取 Kruskal-Wallis 检验 p 值
if (p_value < 0.001) {
sig <- '***' } else if(p_value >= 0.001 & p_value < 0.01) {
sig <- '**' } else if(p_value >= 0.01 & p_value < 0.05) {
sig <- '*' } else {
sig <- '' }
stat_kruskal <- rbind(stat_kruskal, c(i,j,p_value, sig))
#Behrens-Fisher 检验（npmc 包），非参数多重比较
BF_test <- npmc(dat, df = 2, alpha = 0.05)
BF <- BF_test$test$BF
info <- BF_test$info
BF$group <- paste(i, j, sep = '/')
for (q in 1:nrow(BF)){
x <- strsplit(as.character(BF$cmp[q]),"-")[[1]]
BF[q,'class1'] <-rownames(info)[which(info$group.index == x[1])]
BF[q,'class2'] <-rownames(info)[which(info$group.index == x[2])]
}
stat_BF <- rbind(stat_BF, BF[c(12:14, 2:11)])
tmp_BF<-BF[c(12:14, 2:11)]
dat_new <- cbind(aggregate(dat$var, by = list(dat$class), FUN = median),
aggregate(dat$var, by = list(dat$class), FUN = sd))
dat_new$group <- i
dat_new$var <- j
dat_new <- dat_new[c(5, 6, 1,2, 4)]
names(dat_new) <- c('group', 'var', 'class', 'median',"sd")
dat_new <- dat_new[order(dat_new$median, decreasing = TRUE), ]
class <- as.vector(dat_new$class)
m = 1; n = 2; l = 1
dat_new[m,'sig'] <- letters[l]
while(n<length(class)){
for(n in (m+1):length(class)) {
tmp1=info$group.index[which(info$class.level == class[m])]
tmp2=info$group.index[which(info$class.level == class[n])]
if(tmp1<tmp2){
index1=class[m]
index2=class[n]
}else{
index1=class[n]
index2=class[m]
}
tmp_p= tmp_BF[intersect(which(tmp_BF$class1==index1),which(tmp_BF$class2==index2)),]$p.value.2s
if(tmp_p<0.05){
dat_new[n,'sig'] <- letters[l+1]
if (n - m != 1) {
for (x in (m+1):(n-1)) {
tmp1=info$group.index[which(info$class.level == class[x])]
tmp2=info$group.index[which(info$class.level == class[n])]
if(tmp1<tmp2){
index1=class[x]
index2=class[n]
}else{
index1=class[n]
index2=class[x]
}
tmp_p= tmp_BF[intersect(which(tmp_BF$class1==index1),which(tmp_BF$class2==index2)),]$p.value.2s
if(tmp_p<0.05){
dat_new[x,'sig'] <- letters[l]
}else{
dat_new[x,'sig'] <- paste(letters[l], letters[l+1], sep = '')
}
}
}
l = l + 1
m = n
break
}
}
dat_new[(m:n),'sig'] <- letters[l]
}
data_list <- rbind(data_list, dat_new)
}
}
整体的统计检验结果：
每对里面每个类别的Behrens-Fisher非参数多重比较结果：
最后包含abc的文件如下图所示：
我们按分面画出每对的统计检验结果。
ggplot(data = data_list, aes(x = class, y = median)) +
geom_col(aes(fill = class), color = NA, show.legend = FALSE) +
geom_errorbar(aes(ymin = median - sd, ymax = median + sd), width = 0.2) +
facet_grid(var~group, scale = 'free' , space = 'free_x') +
geom_text(aes(label = sig, y = median + sd + 0.3*(median+sd)),size=5) +
labs(x = '', y = '')+
theme(text = element_text(size = 20))
