# This program is written specifically for the teenage deviant
# behavior data as the covariate has to be discreet with values
# 4,5,...,21.

#Modification can be made for other kinds of discreet covariate

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

wayne<- read.table("data",header=T)
wayne[,5][wayne[,5]==9]<-8
wayne[,5][wayne[,5]==1]<-2
wayne[,6][wayne[,6]==1]<-2
sizes<-as.vector(table(wayne[,4],wayne[,5],wayne[,6]))
sizea<-table(wayne[,6])
sizeb<-table(wayne[,5])
sizec<-table(wayne[,4])
wayne2<-wayne[order(wayne[,6],wayne[,5],wayne[,4]),]
y<-wayne2[,2]
x<-wayne2[,3]

ja<-t(rep(1,5))/5
jb<-t(rep(1,7))/7
jc<-t(rep(1,2))/2
pa<-diag(rep(1,5))-t(matrix(sizea/n,5,5))
pb<-diag(rep(1,7))-t(matrix(sizeb/n,7,7))
pc<-diag(rep(1,2))-t(matrix(sizec/n,2,2))
#Factor A is parental education
#Factor B is grade
#Factor C is gender
ha.mat<-kronecker(pa,kronecker(jb,jc))
hb.mat<-kronecker(kronecker(ja,pb),jc)
hc.mat<-kronecker(ja,kronecker(jb,pc))
n <- sum(sizes)
k <- length(sizes)
d<-diag(1/(sizes-1))
m<-solve(diag(rep(1,k))+2*d)

group <- rep(1:k,sizes)


#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
#            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                                               #
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
#P-value for testing  that higher parental education then less marijuna use
 s<-sizea/n*(length(sizea):1)
 ta <- t(s)%*%ha.mat%*%int
 cvc <- t(s)%*%ha.mat%*%v.mat%*%t(ha.mat)%*%s
 testst<-sqrt(n)*ta/sqrt(cvc)
 1-pnorm(testst)
 testst
#P-value for testing  that higher grade then less marijuna use
 s<-sizeb/n*(length(sizeb):1)
 ta <- t(s)%*%hb.mat%*%int
 cvc <- t(s)%*%hb.mat%*%v.mat%*%t(hb.mat)%*%s
 testst<-sqrt(n)*ta/sqrt(cvc)
 1-pnorm(testst)
 testst

#P-value for testing that females have less marijuna use
 s<-sizec/n*(length(sizec):1)
 ta <- t(s)%*%hc.mat%*%int
 cvc <- t(s)%*%hc.mat%*%v.mat%*%t(hc.mat)%*%s
 testst<-sqrt(n)*ta/sqrt(cvc)
 1-pnorm(testst)
 testst