1.将K组数据混合排序,得到各观测值在混合样本中的秩,记第组样本的第个观测值订的秩为 Rii·对各组观测值的秩求和得Ri=n=i Ri, i=l,..,k.该组数据的平均秩为 R; = Ri/ni:所有数据的平均秩为 R=Z=1 Ri/N=(N+1)/2.对于固定的样本n1,,nk,共有M = N!/ II'=1 ni!种混合排序结果,在零假设下,每种结果出现的概率应该相等,概率值为1/M
1. 将 组数据混合排序, 得到各观测值在混合 样本中的秩. 记第 组样本的第 个观测值 的 秩为 . 对各组观测值的秩求和得 该组数据的平均秩为 所有数据的平均秩为 对于固定的样本 , 共有 种混合排序结果. 在零假设下, 每种结果出现的概 率应该相等, 概率值为
2.Kruskal-Wallis统计量kh1212R2- 3(N + 1)ZZCn;(R; - R)?=HN(N + 1)N(N+1)ni=1i=1在零假设Ho :Fi(α) =···= F(α) = F(α)下,按每种混合排序结果的概率为1/M,可计算Kruskal-Wallis统计量H ,在零假设的分布.(n1, n2, n3)书上对不同和水平α,给出了满足P(H ≥)= α. 的临界值c表
2. Kruskal-Wallis统计量 在零假设 下, 按 每种混合排序结果的概率为 , 可计算 Kruskal-Wallis 统计量 , 在零假设的分布. 书上对不同 和水平 , 给出了满 足 的临界值 c 表
3. 将观测值代入统计量H,得到相应的数值4.比较零分布和观测值相应的统计量H的值,5.做出是否拒绝零假设的决定.在N大时,并且对每个i,ni/N趋于某个非零数入i≠0统计量H在零假设下近似于x(%-1)分布另外在大样本时,统计量Zh=1 n N+1Ri-/(k-1)2F*Zh=1 Zn=i(Rij - R)2/(N - k)近似于F(k-1,N-k)分布
3. 将观测值代入统计量 , 得到相应的数值. 4. 比较零分布和观测值相应的统计量 的值. 5. 做出是否拒绝零假设的决定
123生活方式例中,3.7(3.5)7.3 (12)9.0 (14)一个月后3.7(3.5)5.2 (7)4.9 (6)减少的重量3.0(2)5.3 (8)7.1 (11)(5)(单位500k)3.95.7 (9)8.7 (13)及秩2.7(1)6.5 (10)15秩和Ri46443秩平均R119.2Ho: Q1= 02= U3;H1:至少有一个等式不成立这里N=14.算出H大约等于9.4114用精确检验方法得到p-值:0.001347858查表找到P(H≥8.52)=0.0048而P(H≥9.51)=0.001031(用大样本近似,p-值0.00895)
例中, 至少有一个等式不成立 用精确检验方法得到p-值:0.001347858 查表 ( 用大样本近似,p-值0.00895)
KW.test=function(m1=5,m2=5,m3=4,Hvalue=9.4114)# this programis for m1=5, m2=5, and m3 can be any integerm<-m1+m2+m3;Jh5=function(m)(a<-rep(0,5)for (i in 1:(m-4)(for (j in (i+1):(m-3)for (k in (j+1):(m-2))(for (l in (k+1):(m-1))[for (f in (I+1):m)[a<-rbind(a,c(ij,k,I,f))))a[2:nrow(a),J)JTid1<-Jh5(m1+m2+m3);n1<-nrow(JTid1);JTid2<-Jh5(m2+m3);n2<-nrow(JTid2);nn<-n1*n2;const<-1:m;y<-0for (i in 1:n1)[for (j in 1:n2)[temp1<-c(JTid1[i,]);temp2<-(const[-temp1])[c(JTid2i,])];temp3<-const[-c(temp1,temp2)];y<-c(y,12/(m*(m+1)*((sum(temp1)^2/m1+(sum(temp2)^2/m2+(sum(temp3))^2/m3)-3*(m+1))))y<-y[2:(nn+1)];pvalue<-(sum(y>=Hvalue))/nn;y<-sort(y);aaa<-aa<-y[1];tempc<-1for (i in 2:nn)(if (y[i]-aa)>10^{-12)[aaa<-c(aaa,y[i]);aa<-y[i];tempc<-c(tempc,1-(i-1)/nn)))out<-cbind(aaa,tempc)list(c("(m1,m2,m3)"=c(m1,m2,m3),"H"=Hvalue,"pval"=pvalue),out))
KW.test=function(m1=5,m2=5,m3=4,Hvalue=9.4114){ # this program is for m1=5, m2=5, and m3 can be any integer m<-m1+m2+m3;Jh5=function(m){a<-rep(0,5) for (i in 1:(m-4)){for (j in (i+1):(m-3)){ for (k in (j+1):(m-2)){for (l in (k+1):(m-1)) {for (f in (l+1):m){a<-rbind(a,c(i,j,k,l,f))}}}}} a[2:nrow(a),]} JTid1<-Jh5(m1+m2+m3);n1<-nrow(JTid1);JTid2<-Jh5(m2+m3); n2<-nrow(JTid2);nn<-n1*n2;const<-1:m;y<-0 for (i in 1:n1){for (j in 1:n2){temp1<-c(JTid1[i,]); temp2<-(const[-temp1])[c(JTid2[j,])]; temp3<-const[-c(temp1,temp2)]; y<-c(y,12/(m*(m+1))*((sum(temp1))^2/m1+(sum(temp2))^2/m2+ (sum(temp3))^2/m3)-3*(m+1))}} y<-y[2:(nn+1)];pvalue<-(sum(y>=Hvalue))/nn; y<-sort(y);aaa<-aa<-y[1];tempc<-1 for (i in 2:nn){if ((y[i]-aa)>10^{-12}) {aaa<-c(aaa,y[i]);aa<-y[i];tempc<-c(tempc,1-(i-1)/nn)}} out<-cbind(aaa,tempc) list(c("(m1,m2,m3)"=c(m1,m2,m3),"H"=Hvalue,"pval"=pvalue),out)}