##########################################################################
#########multiple-sample Rank-Sum Test for clustered data #######################
##########################################################################
###Note: this function works for m groups and M clusters
###Note: the input is Cluster(from 1 to M), X(the value), 
### and grp(the Groups are from 1 to m) 
### Data must be arranged in the way illustrated in the sample below
### the output of the function is S.j, E(S.j), the chisq-statistic
### with degree of freedom of m-1, and p-values
### Written by: Xu, Ping
### Distributed "as is" without warrantee of any kind; use at your own risk!
##########################################################################


clus.rank.sum<-function(Cluster,X,grp) {

#####calculate quantity 2 (using the pooled estimate of F)
      n<-length(X)
      F.hat<-numeric(n)
      for (i in 1:n){
           F.hat[i]<-(sum(X<=X[i])+sum(X<X[i]))/(2*n)
                    }
#####calculate quantity 1 (using ECD-F for each cluster)
#### M is No. of clusters, n is No. of observations
     M<-length(unique(Cluster)) 
     n.i<-table(Cluster)
     F.prop<-numeric(n)
     for(ii in 1:n){
               F.j<-numeric(M)
               for (i in 1:M){
                    F.j[i]<-(sum(X[Cluster==i]<X[ii])+0.5*sum(X[Cluster==i]==X[ii]))/(n.i[i])
                             } 
               F.prop[ii]<-sum(F.j[-Cluster[ii]])
                   }   

###########calculate S(j)=E(W*|X,g=j), where m is the number of groups
m<-length(unique(grp))
a<-matrix(0,m,M)
b<-1+F.prop
for(j in 1:m){
    for (i in 1:M){
      gik.j<-ifelse(grp==j,1,0)
      a[j,i]<-sum((gik.j[Cluster==i]*b[Cluster==i])/(n.i[i]))
                  } 
              }   
c<-1/(M+1)
S.j<-c*(apply(a,1,sum))

#########Calculate E(S)=E(W*)
n.ij<-matrix(0,m,M)
for (i in 1:m){
n.ij[i,]<-table(Cluster[grp==i])           
}
d<-apply(n.ij,1,FUN=function(x){x/n.i})
E.S.j<-(1/2)*(apply(d,2,sum))

#######Calculate estimate of variance of S
W.hat<-matrix(0,m,M) 
a<-t(d)       #####first calculate W.hat for each cluster
for (i in 1:M){
     for (j in 1:m){
         gik.j<-ifelse(grp[Cluster==i]==j,1,0)
         b<-1/(n.i[i]*(M+1))
         c<-(gik.j)*(M-1)
         d<-sum(a[j,-i])
         W.hat[j,i]<-b*sum((c-d)*F.hat[Cluster==i])
              }
                  }
E.W<-matrix(0,m,M)
for (j in 1:m){
E.W[j,]<-(M/(2*(M+1)))*(a[j,]-sum(a[j,])/M)
              }                            ##second, calculate E(W)

##########calculate sample variance
dev.W<-W.hat-E.W
term.old<-matrix(0,m,m)
for (i in 1:M){
    term<-dev.W[,i]%*%t(dev.W[,i])
    term.old<-term+term.old
              }
V.hat<-(1/M)*term.old

######calculate the test statistic
library(MASS)
T<-(t(S.j-E.S.j)%*%ginv(V.hat)%*%(S.j-E.S.j))*(1/M)
p.value<-p.value<-pchisq(T, df=(m-1),lower.tail=F)

 #calculate the test statistic
list(S=S.j,E.S=E.S.j,Chisq.stat=T,df=m-1,p.value=p.value)
      }




################### test using the data in the table
#######creat a sample data set following table 1 in the paper

ID<-c(1,2,3,4,5,6,7,8,9)
Cluster<-c(1,1,2,2,2,2,3,3,3)
Member<-c(1,2,1,2,3,4,1,2,3)
X<-c(1,4,2,4,6,7,4,7,8)
grp<-c(1,2,1,1,2,2,2,1,2)
my.data<-data.frame(cbind(ID,Cluster,Member,X,grp))

#########obtain chisq-stat using function clus.rank.sum
clus.rank.sum(Cluster,X,grp)

> clus.rank.sum(Cluster,X,grp)
$S
[1] 0.578125 0.921875

$E.S
[1] 0.6666667 0.8333333

$Chisq.stat
         [,1]
[1,] 1.392647

$df
[1] 1

$p.value
          [,1]
[1,] 0.2379586