非参数检验之KW单因素方差分析和JT检验

knitr::opts_chunk$set(tidy = TRUE, 
                      warning = FALSE,
                      message = FALSE)

setwd("C:/Users/213yi/Desktop/非参数统计/4-12")
library(showtext) #载入库
library(ggplot2)
library(MASS)
library(dunn.test)
library(ggthemes)
library(ggplot2)
library(rJava)
library(xlsx)
library(plyr)
library(epade)
library(MASS)
library(RColorBrewer)
library(datarium)

药物治疗效果的检验

KW单因素方差分析(H检验)

drug=c(80,203,236,252,284,368,457,393,133,180,100,160,156,295,320,448,465,481,279,194,214,272,330,386,475)
gr.drug=factor(rep(1:4,c(8,4,7,6)),labels = c("A","B","C","D"))
kruskal.test(drug,gr.drug)
#结果显著

编写Dun函数

#Dunn_test(drug, k = 4, nj = c(8, 4, 7, 6))
Dunn_test<-function(data,k,nj,alpha=0.05){
  #输出的是谁和谁比较，djj是多少（相除得到的），Z1-α的值，
  #第i个处理是从第几个到第几个，去索引rank
  n<-length(data)
  r<-rank(data)
  sy<-c()
  label = c("A","B","C","D","E","F","J","H","I")
  
  for (j in 1:k){
    sy[j]<-cumsum(nj)[j]-nj[j]+1
  }
  #sy=1  9 13 20
  rrankmean=c()
  rsum=0
  j=0
  for (i in 1:n){
    if ((i %in% sy)&(i!=1)){
      j=j+1
      rsum=rsum/nj[j]
      rrankmean<-c(rrankmean,rsum)
      rsum=r[i]
    }else {
      rsum=rsum+r[i]
    }
  }
  rrankmean<-c(rrankmean,rsum/nj[k])
  
  resolution=as.data.frame(matrix(ncol = 3))
  colnames(resolution)<-c("比较","dij","Z(1-α*/2)")
  
  
  for (i in 1:(k-1)){
    for (j in (i+1):k){
      con<-c()
      con[1]<-paste(label[i],"vs",label[j])
      con[2]<-abs(rrankmean[i]-rrankmean[j])/sqrt(n*(n+1)/12*(1/nj[i]+1/nj[j]))
      con[3]<-qnorm(1-2*alpha/(k-1)/k/2)#Bonferroni修正
      resolution<-rbind(resolution,con)
    }
  }
  resolution=resolution[-1,]
  return(resolution)
}

对药效进行检验修正

Dunn_test(drug, k = 4, nj = c(8, 4, 7, 6))

• 由于加了绝对值，其实是双边检验，α*要除以2

• 也可以这样理解：是因为绝对值的加入只能>0,因此只有0.5的总概率

• 可以看出，由于2.75618>2.63825

• 只有B vs C的样本拒绝了原假设，即B和C的样本分布不同

内置Dun函数的检验(修正 Bonferroni)

dunn.test(drug,gr.drug,method = "bonferroni")

• 通过内置函数的检验，发现显著的是B和C
• 计算的$d_{j,j^{、}}$=-2.756186和上述的绝对值相符
• 说明自行编写的Dunn函数正确

可视化boxplot

data<-data.frame(drug,gr.drug)
boxplot(drug~gr.drug,data)

• 可以看出，B和C确实是最离谱的分布，证明了Dunn检验的真实性

收入和学历的关系

参数方差分析

filee=read.table("employee.txt",header = TRUE)
plot(density(filee$salary),main="salary分布密度图")
boxplot(salary~educ,data=filee,main="edu对salary的影响箱线图")

• 从上述箱线图可以看出，高学历的人群收入无论是下分位点、还是上分位点或者均值，都相对较高
• 由于参数方差分析要求每个处理都是正太的样本，但是通过上图发现，并不是完全如此
• 下面，对其中部分进行检验

前提检验

#对受教育年限17年的进行正态性检验
dt=filee[which(filee$educ==17),2]
plot(dt)
shapiro.test(dt)
ks.test(dt,rnorm(10000,mean(dt),sd(dt)))

• 发现：样本量过少，定量的正态性检验结果做不出来

aov

aov_result <- aov(salary ~ educ, data = filee)
summary(aov_result)

• 这个检验是应用统计学学过的参数方差检验
• 用到的是SST(总平方和)、SSE(误差组内的平方和)、SSA(处理平方和)
• 在 0.01的水平下是显著的
• 方差分析表明收入和教育有关

进行多重比较

filee<-read.table("employee.txt",header = TRUE)
fi=filee
fi$educ[filee$educ <= 12] <- 'A'
fi$educ[filee$educ >= 13 & filee$educ <= 16] <- 'B'
fi$educ[filee$educ >= 17] <- 'C'
fi$educ <- factor(fi$educ, labels = c('A', 'B', 'C'))

#fi是聚类类后，filee是聚类前
#fi可以，filee不行
pairwise.t.test(fi$salary,fi$educ)

• 配对 t 检验结果表明低学历和中等学历以及低学历和高等学历之间的工资水平存在显著差异

非参数方差分析

filee<-read.table("employee.txt",header = TRUE)
fi=filee
fi$educ[filee$educ <= 12] <- 'A'
fi$educ[filee$educ >= 13 & filee$educ <= 16] <- 'B'
fi$educ[filee$educ >= 17] <- 'C'
fi$educ <- factor(fi$educ, labels = c('A', 'B', 'C'))

kruskal.test(fi$salary,fi$educ)

• 参数分布有要求：正态性，独立同方差的假设
• 而非参数KW单因素方差分析只用到了rank
• 只是假定分布连续、分布形状相似
• 更可靠，此时结果拒绝他们分布相同的原假设
• 接下来进行Dunn的检验

Dunn检验

fileeorder=filee[order(filee$educ),]
Dunn_test(fileeorder$salary,7,c(4,11,7,3,3,1,1))

• Dun的结果不理想

• 没有任意2个样本拒绝了分布均值不同的假设

• 所以，下面还是对样本的类别重新进行聚类

filee<-read.table("employee.txt",header = TRUE)
fi=filee
fi$educ[filee$educ <= 12] <- 'A'
fi$educ[filee$educ >= 13 & filee$educ <= 16] <- 'B'
fi$educ[filee$educ >= 17] <- 'C'
fi$educ <- factor(fi$educ, labels = c('A', 'B', 'C'))

Dunn_test(fi$salary,k=3,c(15,10,5))

• Dunn 检验结果也表明低学历和中等学历以及低学历和高等学历之间的工资水平存在显著差异

JT检验

JTnorm<-function(x,group){
  N<-length(x)
  index<-unique(group)
  k<-length(index)
  ns<-NULL
  for(i in 1:k){
    ns<-c(ns,sum(group==index[i]))
  }
  Diffval<-NULL
  for(i in 1:(k-1)){
    xi<-x[which(group==index[i])]
    for(j in(i+1):k){
      xj<-x[which(group==index[j])]
      Diffvali<-0
      for(l in 1:length(xj)) {
        Diffvali<-Diffvali+sum(xi-xj[l]<0)
      }
      Diffval<-c(Diffval,Diffvali)
    }
  }
  m.val<-(N^2-sum(ns^2))/4
  sd.val<-sqrt((N^2*(2*N+3)-sum(ns^2*(2*ns+3)))/72)
  zval<-(sum(Diffval)-m.val)/sd.val 
  pval<-pnorm(zval,0,1,lower.tail = FALSE)
  list(ns=ns,k=k,Diffval=Diffval,J=sum(Diffval),m.val=m.val,sd.val=sd.val,zval=zval,pval=pval)
}

验证例4.5

jump=c(125,136,116,101,105,109,122,114,131,120,119,127,128,142,128,134,135,131,140,129)
gr.jump=factor(rep(1:3,c(6,6,8)))
JTnorm(jump,gr.jump)

• 注意 为什么我的lower.tail=F失效，必须是FALSE
• JT结果很显著0.0008，即样本的位置呈现出上升或下降趋势

手写作业

定理 4.1 的证明（提示：根据期望、方差和协方差定义出发进行证明，可以参考定理 3.1 或者第一章 P22 页推论 1.3 的证明）

P143 页：习题 4.1，要求写出 H 检验过程(使用大样本卡方分布近似，并用 R 函数 kruskal.test()来验证求解过程是否正确)

drug=c(83,64,67,62,70,85,81,80,78,88,89,79,90,95)
gr.drug=factor(rep(1:3,c(5,4,5)),labels = c("A","B","C"))
kruskal.test(drug,gr.drug)
#结果显著