## Data will be a three column format; the first two columns will be the paired responses
## and the third column will be the cluster number

data <- read.table("C:\\Somnath\\Research\\signed rank\\alzheimer.txt",header=F)

Xij <- data[,1]-data[,2]
ni <- as.vector(table(data[,3]))
g <- length(ni)
n <- sum(ni)

 


cni <- cumsum(ni)
cni <- c(0,cni)

Fi <- function(x,i) { Xi <- Xij[(cni[i]+1):(cni[i+1])];
                      (sum(abs(Xi)<=x)+sum(abs(Xi)<x))/(2*ni[i])}
Ftot <- function(x) { st <- 0;
                      for (i in 1:g) st <- st + Fi(x,i);
                      return(st)}
Fcom <- function(x) { st <- 0;
                      for (i in 1:g) st <- st + Fi(x,i)*ni[i];
                      return(st/n)}

# SIGNED RANK TEST STATISTIC
TS <- VTS <- 0
for (i in 1:g) {

Xi <- Xij[(cni[i]+1):(cni[i+1])]
first <- (sum(Xi>0)-sum(Xi<0))/length(Xi)
second <- 0
third <- 0
for (x in Xi) { second <- second + sign(x)*(Ftot(abs(x))-Fi(abs(x),i));
                third <- third + sign(x)*Fcom(abs(x))}

TS <- TS + first+second/length(Xi)
VTS <- VTS + (first+ (g-1)*third/length(Xi))^2
                }

Z <- TS/sqrt(VTS)

unlist(list("Test Statistic"=Z, "P -value"= 2*(1-pnorm(abs(Z))) ))