######################Adaptive Sum test with 2-directional searches to test for for association, possibly including interactions
######################   in the possible presence of nuisance parameters. 
###########References: 
#            Pan W, Basu S, Shen X (2011). Adaptive Tests for Detecting Gene-Gen
e and Gene-Environment Interactions.

######################Wei Pan, Division of Biostatistics, U of Minnesota
######################Email: weip@biostat imn.edu
#### V1, May 25, 2011


library("mvtnorm")

#Input:
#      Y: a vector of 0-1 response, nx1;
#      X: design matrix for covariates (e.g. gene-gene interactions) of
#           interest (i.e. to be tested), nxk;
#      Z: design matrix for nuisance covariates, nxk2;
#      B: # of simulations to calculate p-values
#Goal:
#      fit model: logit(Pr(Y=1)) = a0+ a1*X+a2*Z ,
#      test H0: a1=0.
#Output: p-values for aSum with the natural order (i.e. order given by the
#        columns of X), aSum with the order determined by Uj/SEj,
#        aSum with the order opposite to that determined by Uj/SEj,
#        and aSUm2d with 2 directional searches ordered by Uj/SEj;
#        In the notation given in Pan et al (2011), the OUTPUT is a vector
#        of p-values for aSum-Loc, aSum-Ord, aSum-ReverseOrd, 
#        aSum2-Ord (=min{P_aSum-Ord, P_aSum-ReverseOrd)

#Examples: suppose that one has two groups of variables in design matrices
#          X1 and X2 respectively, and Y is a binary response vector,
#          1) if one would like to test the marginal effects of X1, then
#             Tests5V1(Y, X=X1);
#          2) if one would like to test the effects of X1 after ajusting for X2,
#             Tests5V1(Y, X=X1, Z=X2);
#          3) if one would like to test the effects of X1 in the presence of 
#             possible interactions with X2, then
#             Tests5V1(Y, X=cbind(X1, X12), Z=X2), 
#             where X12 is the design matrix for interactions;
#          4) if one would like to test on the interaction between X1 and X2
#             after adjusting for the main effects of X1 and X2, then
#             Tests5V1(Y, X=X12, Z=cbind(X1,X2)), 
#             where X12 is the design matrix for interactions;

############################################
############################################
############################################

aSum2dV1<-function(Y, X, Z=NULL, B=200){

X<-as.matrix(X)

n<-length(Y)
k<-ncol(X)
k2<-ncol(Z)

#######construction of the score vector and its cov matrix:
if (is.null(Z)){
## NO nuisance parameters:
   Xg <- X
   Xbar<-apply(Xg, 2, mean)
   Xgb<-Xg
   for(i in 1:nrow(Xg))
      Xgb[i,]<-Xg[i,]-Xbar
   ##score vector:
   U<-t(Xg) %*% (Y-mean(Y))
   ##cov of the score stats:
   CovS<- mean(Y)*(1-mean(Y))*(t(Xgb) %*% Xgb)

} else {
## with nuisance parameters:
tdat1<-data.frame(trait=Y, Z)
fit1<-glm(trait~.,family="binomial",data=tdat1)
pis<-fitted.values(fit1)
Us<-matrix(0, nrow=n, ncol=k)
for(i in 1:k){
    tdat2<-data.frame(X1=X[,i], Z)
    fit2<-glm(X1~.,data=tdat2)
    X1mus<-fitted.values(fit2)
    Us[, i]<-(Y - pis)*(X[,i] - X1mus)
    }
U<-apply(Us, 2, sum)
CovS<-matrix(0, nrow=k, ncol=k)
for(i in 1:n)
  CovS<-CovS + Us[i,] %*% t(Us[i,])
}

K=length(U)

      svd.CovS<-svd(CovS)
      CovSsqrt<-svd.CovS$u %*% diag(sqrt(svd.CovS$d))

   #order U's components from large to small in Z-stats:
      Ords<-order(U/sqrt(diag(CovS)), decreasing=T)
      Uords<-U[Ords]
      CovSords<-CovS[Ords, Ords]

   Sums<-Sumords<-Sumords2<-rep(0, length(U))
   for (k in 1:length(U)){
     Sums[k]<-Sum(U[1:k], CovS[1:k,1:k])
     Sumords[k]<-Sum(Uords[1:k], CovSords[1:k,1:k])
     Sumords2[k]<-Sum(Uords[K:(K-k+1)], CovSords[K:(K-k+1), K:(K-k+1)])
     }
   Sum1<-min(Sums)
   Sumord1<-min(Sumords)
   Sumord2<-min(Sumords2)
 
   #######simulations:
   n<-length(Y)
   Sum0s<-Sum0ords<-Sum0ords2<-Sum0ordsBoth<-rep(0, B)
   for(b in 1:B){
      U00<-rnorm(K, 0, 1)
      U0<-CovSsqrt %*% U00

     #order U's components from large to small in Z-stats:
      Ords<-order(U0/sqrt(diag(CovS)), decreasing=T)
      U0ords<-U0[Ords]
      CovS0ords<-CovS[Ords, Ords]

     for (k in 1:length(U0)){
       Sums[k]<-Sum(U0[1:k], CovS[1:k,1:k])
       Sumords[k]<-Sum(U0ords[1:k], CovS0ords[1:k,1:k])
       Sumords2[k]<-Sum(U0ords[K:(K-k+1)], CovS0ords[K:(K-k+1), K:(K-k+1)])
       }
     Sum0s[b]<-min(Sums)
     Sum0ords[b]<-min(Sumords)
     Sum0ords2[b]<-min(Sumords2)
     Sum0ordsBoth[b]<-min(c(Sum0ords[b], Sum0ords2[b]))
     }

return(c(
         lessThan(Sum0s, Sum1)/B,
         lessThan(Sum0ords, Sumord1)/B,
         lessThan(Sum0ords2, Sumord2)/B,
         lessThan(Sum0ordsBoth, min(c(Sumord1, Sumord2)) )/B
        ))

}  

lessThan<-function(x, x0){
  if (x0==0) sum(x <= x0)
  else sum(x<x0)
  }

##########SumTest########################
Sum<-function(U, CovS){
     #it's possible U=0 and Cov(U)=0:
     if (abs(sum(U))<1e-10 || sum(CovS)<1e-10) pTsum<-1 else{
      a<-rep(1, length(U))
      Tsum<- sum(U)/(sqrt(as.numeric(t(a) %*% CovS %*% (a))))
      pTsum <- as.numeric( 1-pchisq(Tsum^2, 1) )
      }
pTsum
}