#================== 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 distributino based on the parametric bootstrap
 


aSPUboot <- 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)
   
  if (is.null(Z)){
  ## NO nuisance parameters:
   Xg <- X
   U<-t(Xg) %*% (Y-mean(Y))
  } else {
  tdat1<-data.frame(trait=Y, Z)
  fit1<-glm(trait~.,family="binomial",data=tdat1)
  pis<-fitted.values(fit1)
  Us<-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])
  }
  U<-t(XUs) %*% (Y - pis)
 }

  ##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")
  ## perm
  pPerm0 = rep(NA,length(pow))
  T0s = numeric(B)
  s <- sample(1:10^5,1)
  Y0 <- numeric(n)

  for (j in 1:length(pow)){
    set.seed(s) # to ensure the same samples are drawn for each pow
      for (b in 1:B){
            if (is.null(Z)) {
            Y0 <- sample(Y, length(Y))
            ##  Null score vector:
            U0<-t(Xg) %*% (Y0-mean(Y0))
            } else {
            ## with nuisance parameters:
            for(i in 1:n) Y0[i] <- sample(c(1,0), 1, prob=c(pis[i], 1-pis[i]) )
            tdat0<-data.frame(trait=Y0, Z)
            fit0<-glm(trait~.,family="binomial",data=tdat0)
            pis0<-fitted.values(fit0)
            U0<-t(XUs) %*% (Y0 - pis0)
            }

        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)
}