date()


#########################################################################
#Generating data                                                         
#The k samples are stored in one vector y=(y_11,...,y_1n1,...,y_k1,...,  
#y_knk)' and the covariates are stored in one vector x similarly.        
#The vector SIZES contains the sizes for each sample. k is the number of 
#samples.                                                              
#########################################################################
r1<-numeric(5000)
r2<-numeric(5000)
r3<-numeric(5000)
sizes<-c(16,24,12,20,18,15)
n <- sum(sizes)
k <- length(sizes)
ctr.mat1 <- matrix(c(1/6,1/6,1/6,-1/6,-1/6,-1/6),1,6)
ctr.mat2<-matrix(rep(c(1/3,-1/6,-1/6),4),2,6)
ctr.mat3<-kronecker(matrix(c(1,-1),1,2),cbind(rep(1,2),-diag(rep(1,2))))

set.seed(368)

group <- rep(1:k,sizes)
#group

allsizes <- rep(sizes,sizes)
#allsizes


an <- 9/allsizes^0.26
#an


for(i in 1:5000)
{
index<-sample(1:length(tbbm),n)
x<-htvel[index]
y<-tbbm[index]



#########################################################################
#Calculattion of the test statistic;
#########################################################################


 midrank <- rank(y)  #The midranks of the responses
 
 x.mat <- t(matrix(x,n,n))
 k.mat <- (x.mat-x)/an
 k.mat <- replace(k.mat, abs(k.mat)>1,1)
 k.mat <- (15/16)*(1-k.mat^2)^2
 #k.mat

 temp1 <- rowsum(k.mat,group)
 temp <- replace(temp1,temp1==0,1)

 tmp <- apply(temp,2,rep,times=sizes)
 w.mat <- k.mat/tmp
 #w.mat

 wr.mat<- (1/n)*w.mat*midrank

 int.pre <- apply(wr.mat,1,mean)
 int <- rowsum(int.pre,group)
# int

########################################################################
#calculate \hat \delta_{1,i}^2 and \hat \delta_{1,i_1,i_2},i=1,...,k
########################################################################

 tp1.mat <- rowsum(wr.mat,group)
 v1.mat<-(n-1)/n*var(t(tp1.mat))

#########################################################################
#calculate \hat \delta_{2,i}^2, i=1,...,k
#########################################################################

 wrr.mat <- wr.mat * midrank
 tp4.mat <- rowsum(wrr.mat,group)
 uij.vec <- (apply(w.mat,1,mean))^2
 tp5.mat <- t(tp4.mat)*uij.vec
 tp6.mat <- rowsum(tp5.mat,group)
 v21.vec <- diag(tp6.mat)
 wrsq.mat <- (rowsum(wr.mat,group))^2
 tp7.mat <- n*rowsum(t(wrsq.mat)*uij.vec,group)
 v22.vec <- diag(tp7.mat)
 v2.mat <- diag(v21.vec-v22.vec)
# v2.mat

#########################################################################
#calculate the estimate of the covariance matrix
#########################################################################

 v.mat <- v1.mat + v2.mat
# v.mat

#########################################################################
#calculate the statistic
#########################################################################

 ta1 <- ctr.mat1%*%int
 cvc.mat1 <- solve(ctr.mat1%*%v.mat%*%t(ctr.mat1))

 ta2 <- ctr.mat2%*%int
 cvc.mat2 <- solve(ctr.mat2%*%v.mat%*%t(ctr.mat2))

 ta3 <- ctr.mat3%*%int
 cvc.mat3 <- solve(ctr.mat3%*%v.mat%*%t(ctr.mat3))

 r1[i]<-n*t(ta1)%*%cvc.mat1%*%ta1
 r2[i]<-n*t(ta2)%*%cvc.mat2%*%ta2
 r3[i]<-n*t(ta3)%*%cvc.mat3%*%ta3
}
print("Treatment effect:")
table(r1>=qchisq(0.90,1))/5000
table(r1>=qchisq(0.95,1))/5000
table(r1>=qchisq(0.99,1))/5000
print("Age effect")
table(r2>=qchisq(0.90,2))/5000
table(r2>=qchisq(0.95,2))/5000
table(r2>=qchisq(0.99,2))/5000
print("Interaction")
table(r3>=qchisq(0.90,2))/5000
table(r3>=qchisq(0.95,2))/5000
table(r3>=qchisq(0.99,2))/5000