#**********************************************************************
# File name: analf.q                                                   *
# For: This program is only useful for data analysis of Wayne Osgood  *
#      data set since it specifically uses the coding of the          *
#      covariate which takes values from 4-18. It can be modifeid     *
#      for any kind of discrete covariates.                           *
# Usage: In the current directory, enter Splus. Then store the        *
#        response variable and covariate in vectors "x" and "y".      *
#        They should be in a one-way layout regardless of the original*
#        design and only one index is used to keep track. They should *
#        also be sorted according to the categories. Then store the   *
#        corresponding sample sizes into the vector "sizes". Ex. for a*
#        2X3 design with sample sizes:                                *
#              A\B   1   2   3                                        *
#               1    10 15  20                                        *
#               2    25 30  35                                        *
#        "y" is arranged like (y_11(1),...,y_11(10),...,y_13(1),...,  *
#        y_13(20),y_21(1),...,y_21(25),...,y_23(1),...,y_23(35)) and  *
#        so is "x". "sizes" is a vector like (10,15,20,25,30,35).     *
#        Finally, define the matrix "h.mat" which is a square matrix  *
#        matrix of the contrast and "h" is the diagonal element.      *
#        Modify the matrix "m" in the program. Now, exit Splus and    *
#        run "Splus BATCH anal.q anal". Using "more  anal" in Unix,   *
#        you can see the P_value using finite-sample correction.      *
# Programmer: Yunling Du                                              *
# Date: 4/16/98                                                       *
#**********************************************************************

weight.f<-function(x){
        k <- length(x)
        an <- 14.8/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)]
}



n <- sum(sizes)
k <- length(sizes)
d<-diag(1/(sizes-1))
m<-solve(diag(rep(1,16))+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)*rep(c(rep(0,3),rep(1,14),0),k)
 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<-sizes/n*v2.mat

#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
#calculate the statisticc                                               #
#-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-#
 
 qh<-t(int)%*%h.mat%*%int
 f1<-h^2*((sum(diag(v2.mat)))^2-2*sum(diag(d%*%m%*%v2.mat%*%v2.mat)))/ sum(
     diag(h.mat%*%v2.mat%*%h.mat%*%v2.mat))
 f2<-((sum(diag(solve(d)%*%v.mat)))^2-2*sum(diag(solve(d)%*%m%*%v.mat%*%
     v.mat)))/sum(diag(solve(d)%*%v.mat%*%v.mat))
 testst<- n/(h*sum(diag(v2.mat)))*qh
 1-pf(testst,f1,f2)