#================== Sum of powered score (SPU) test
# Y: disease lables; =0 for controls, =1 for cases;
# X: genotype data; each row for a subject, and each column for an SNP;
#    the value of each element is the # of the copies for an allele;
# Z: covariates (if any)
# B: number of permutation
#==========================================================================
# Output is the p-values of the SPS test in the order of pow and aSPU
#========================================================================== 
# The null distribution based on the permutation of residuals


aSPUperm <- function(Y, X, Z=NULL,  pow=c(1:8, Inf), B=1000){
	
    n <- length(Y)
    if (is.null(X) && length(X)>0) X=as.matrix(X, ncol=1)
    k <- ncol(X)
	
    #### Score vector:   
   if (is.null(Z)){
   ## NO nuisance parameters:
    XUs<-Xg <- X
    r<-Y-mean(Y)
    U<-as.vector(t(Xg) %*% r)
   } else {
	   tdat1<-data.frame(trait=Y, Z)
	   fit1<-glm(trait~.,family="binomial",data=tdat1)
	   pis<-fitted.values(fit1)
	   XUs<-matrix(0, nrow=n, ncol=k)
	   Xmus = X
	   for(i in 1:k){
	       tdat2<-data.frame(X1=X[,i], Z)
	       fit2<-glm(X1~.,data=tdat2)
	       Xmus[,i]<-fitted.values(fit2)
	       XUs[, i]<-(X[,i] - Xmus[,i])
	       }
	   r<-Y - pis
	   U<-t(XUs) %*% r
  }
  
  ##observed statistics
  Ts=rep(NA,length(pow))
  for (j in 1:length(pow)){
  if (pow[j]<Inf) Ts[j] = sum(U^pow[j]) else Ts[j] = max(abs(U))
  }
  ## cat("statistic calculated","\n")
  
  ## residual permutation
  pPerm0 = rep(NA,length(pow))
  T0s = numeric(B)
  s <- sample(1:10^5,1)
    
  for (j in 1:length(pow)){
    set.seed(s) # to ensure the same samples are drawn for each pow
    for (b in 1:B){
        r0 <- sample(r, length(r))
        U0 <- as.vector(t(XUs) %*% r0)
        if (pow[j] < Inf){ T0s[b] = round(sum( U0^pow[j]), digits = 8) }
        if (pow[j] == Inf) {T0s[b] = round( max(abs(U0)), digits = 8) }
	}

	
   pPerm0[j] = round((sum(abs(Ts[j])<=abs(T0s[1:(B-1)]))+1)/(B), digits=8)
   P0s = (B-rank(abs(T0s))+1)/(B)
   if (j==1) minp0=P0s else minp0[which(minp0>P0s)]=P0s[which(minp0>P0s)]
   # cat("j=",j,"\n")
   }

   cat("P0s caculated","\n")
   Paspu<-(sum(minp0<=min(pPerm0))+1)/(B+1)
   pvs <- c(pPerm0, Paspu)
return(pvs) 
}