drug=c(80、203、236、252、284、368、457、393、133、130、160、156、295、295、465、465、481、279、194、214、273、133、160、160、160、150、160、160、160、160、160、156、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、236、2366、236、236、236、236、236、266、266、266、266、266、26、265、26、265、26、26、26、26、26、26、265、265、265、2656、265656、265、265、265、26、26、26、26、26、26、26、265、265、265、265、26、265、26、26、26、26、26、26、26、26、265、26、26、265、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、26、2 gr.drug=factor(rep(1:4,c(8，4，7，6)，labels = c("A","B","C","D")) kruskal.test(drug,gr.drug) #结果显着 

#Dunn_test(drug, k = 4, nj = c(8, 4, 7, 6)) Dunn_test<-function(data,k,nj,alpha=0.05){   #输出是谁相比，djj多少(相除得到)，Z1-α的值，   #第一个处理是从第几个到第几个，去索引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样本分布不同

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

• 通过内置函数的检查，发现B和C
• 计算的 d j , j 、 d_{j,j^、} dj,j、=-2.756186与上述绝对值一致
• 说明自己写的Dunn函数正确

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_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的检验

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 检验结果也表明低学历和中等学历以及低学历和高等学历之间的工资水平存在显著差异

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)
}

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，即样本的位置呈现出上升或下降趋势

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)
#结果显著