## Time-stamp: <2006-07-26 19:37:41 zaykind> [written by Dmitri Zaykin]
#
# Script parameters:
#
# fname: data file name; 1st column affection status: 0 or 1; then SNPS
# recoded as 1,2,3 for AA,Aa,aa
#
# k1,k2,k3: numbers of principal components to use (e.g. 1, 2, 3,)
#
# nsim: number of permutations for case-control LD p-values
#

# args: DataFile, num PC (k1-k3), NumOfSimulations
DeltaTest <- function(fname, k1, k2, k3, nsim)
{
   c <- ncol(read.table(fname))
   library("MASS")
   # read data
   p <- matrix(scan(fname, quiet=TRUE), ncol=c, byrow=T)
   options(warn=-1)
   co <- p[p[,1] != 1,]
   ca <- p[p[,1] != 0,]
   nco <- nrow(co)
   nca <- nrow(ca)

   cat("\nNumber of cases and controls: ", nca, ",", nco, "\n")

   # get LD matrices
   ca2 <- cor(ca[,2:c], use="pairwise")
   co2 <- cor(co[,2:c], use="pairwise")
   k <- c-1

   # norm-based statistic
   s0 <- Tr(crossprod(ca2 - co2)) / (4*(k-1)*k)
   ei1 <- eigen(ca2)
   ei2 <- eigen(co2)
   
   t0 <- Sangle(t(ei1$vectors[,1:k1]), t(ei2$vectors[,1:k1]))
   v0 <- Sangle(t(ei1$vectors[,1:k2]), t(ei2$vectors[,1:k2]))
   u0 <- Sangle(t(ei1$vectors[,1:k3]), t(ei2$vectors[,1:k3]))
   # Obtain statistics empirical distribution
   cat("# PCAs:", k1, k2, k3, "; Stats: ", s0, t0, u0, v0, "\n")
   cnt1 <- cnt2 <- cnt3 <- cnt4 <- 0
   for(i in 1:nsim) {
       # randomly allocate cases and controls
       idx1 <- sample(1:nrow(p), nca, replace=FALSE)
       idx2 <- c(1:nrow(p))[-idx1]
       ca <- p[idx1, ]
       co <- p[idx2, ]
       # recompute statistics
       ca2 <- cor(ca[,2:c], use="pairwise")
       co2 <- cor(co[,2:c], use="pairwise")
       si <- Tr(crossprod(ca2 - co2)) / (4*(k-1)*k)
       ei1 <- tryCatch(eigen(ca2), error = function(e) NA)
       ei2 <- tryCatch(eigen(co2), error = function(e) NA)
       ti <- tryCatch(Sangle(t(ei1$vectors[,1:k1]), t(ei2$vectors[,1:k1])), error = function(e) NA)
       vi <- tryCatch(Sangle(t(ei1$vectors[,1:k2]), t(ei2$vectors[,1:k2])), error = function(e) NA)
       ui <- tryCatch(Sangle(t(ei1$vectors[,1:k3]), t(ei2$vectors[,1:k3])), error = function(e) NA)

       if(identical(NA,si) 
       || identical(NA,ti) 
       || identical(NA,vi) 
       || identical(NA,ui)) i <- i-1
       else {
          if(si >= s0) cnt1 <- cnt1+1
          if(ti <= t0) cnt2 <- cnt2+1
	  if(vi <= v0) cnt3 <- cnt3+1
          if(ui <= u0) cnt4 <- cnt4+1
       }
   }

   cat("LD norm p-value = ", x1 <- cnt1/nsim, "\n")
   cat("LD subspace1 p-value = ", x2 <- cnt2/nsim, "\n")
   cat("LD subspace2 p-value = ", x3 <- cnt3/nsim, "\n")
   cat("LD subspace3 p-value = ", x4 <- cnt4/nsim, "\n")
   x <- list("pv1"=x1, "pv2"=x2, "pv3"=x3, "pv4"=x4)
   return(x)
}

Tr <- function(x) sum(diag(x))

Sangle <- function(u, v) {
  ifelse(identical(u, NA) || identical(v, NA),
    NA, Tr(u %*% t(v) %*% v %*% t(u)))
}