# This program is used to calibrate bandwidth for the teenage
# deviant behavior data as the covariate is assumed to take integer
# values between 4 and 21

# A one-way design is used

# For bandwidth calibration, change the constant in
#               an<-const./k^0.2
# in the following function
#
# When the covariate is continuous, an<-const./k^0.26

weight.f<-function(x){
        k <- length(x)
        an <- 13.9/k^0.2
        x.mat <- (matrix(4:21, 18, 18, byrow = T) - 4:21)/an
        x.mat <- replace(x.mat, abs(x.mat) > 1, 1)
        k.mat <- (15/16) * (1 - x.mat^2)^2
        ct <- count.f(x)
        wt <- k.mat * ct
        tot <- apply(wt, 2, sum)
        tot <- replace(tot, tot == 0, 1)
        w.mat <- k.mat/tot
        w.mat
}
  
cal1.f<-function(mat,count){
        apply(mat * count, 2, sum)/sum(count)
}

count.f<-function(x){
        tabulate(x, 21)[ - (1:3)]}


r1<-matrix(integer(15000),5000,3)
           # the row vector of r1 will contain 1 (the test
           # significant at alpha=0.10,0.05 or 0.01) and 0
           # (not significant at the levels.
sizes<-rep(15,8) # Specify the sample sizes to be used per group
n<- sum(sizes)
k <- length(sizes)
ctr.mat<-diag(rep(1,8))-matrix(rep(1,64),8,8)/8
           # omnibus test will be used for calibration
set.seed(368)

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


for(i in 1:5000)
{
index<-sample(1:nrow(wayne),n)
x<-(wayne[,3])[index]
y<-(wayne[,2])[index]

#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
#            Calculate the test statistic                         #
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#

 midrank <- rank(y)  #The midranks of the responses
 
 
 w.arr<-lapply(split(x,group),weight.f)
 ct<-count.f(x)
 wm.arr<-lapply(w.arr,cal1.f,count=ct)
 wm.v<-unlist(wm.arr,use.names=F)
 ct1<-unlist(lapply(split(x,group),count.f),use.names=F)
 wm.v<-rep(wm.v,ct1)
 midrank<-midrank[order(group,x)]
 vmr.v<-(1/n)*wm.v*midrank
 
 int <- rowsum(vmr.v,group)

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

 ww.mat<-matrix(unlist(w.arr,use.names=F),18,18*k)
 ww.mat<-apply(ww.mat,1,rep,times=ct1)
 ww.mat<-ww.mat*midrank/n
 wt2.mat<-rowsum(ww.mat,group)
 wt3.mat<-t(wt2.mat)
 wt.mat<-apply(wt3.mat,2,rep,times=ct)
 v1.mat<-(n-1)/n*var(wt.mat)

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

 wt1.mat<-ww.mat*midrank/n
 wt1.mat<-rowsum(wt1.mat,group)
 wt.mat<-wt1.mat-(wt2.mat)^2

 vu.mat<-rep(as.vector(t(wt.mat)),times=ct1)*wm.v^2
 vu.mat<-rowsum(vu.mat,group)
 v2.mat<-diag(vu.mat)*n
 

#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
#calculate the estimate of the covariance matrixc                       #
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#

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

#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
#calculate the statisticc                                               #
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
 s<-sizes/n*(1:k)
 ta <- t(s)%*%ctr.mat%*%int
 cvc.mat <- t(s)%*%ctr.mat%*%v.mat%*%t(ctr.mat)%*%s
 testst<-sqrt(n)*ta/sqrt(cvc.mat)
 qt<-qnorm(c(0.9,0.95,0.99))
 r1[i,]<-(testst>qt)
}

apply(r1,2,table)/5000