R/test.R

In multtest: Resampling-based multiple hypothesis testing

#main user-level function for multiple hypothesis testing

MTP<-function(X,W=NULL,Y=NULL,Z=NULL,Z.incl=NULL,Z.test=NULL,na.rm=TRUE,test="t.twosamp.unequalvar",robust=FALSE,standardize=TRUE,alternative="two.sided",psi0=0,typeone="fwer",k=0,q=0.1,fdr.method="conservative",alpha=0.05,smooth.null=FALSE,nulldist="boot.cs",B=1000,ic.quant.trans=FALSE,MVN.method="mvrnorm",penalty=1e-6,method="ss.maxT",get.cr=FALSE,get.cutoff=FALSE,get.adjp=TRUE,keep.nulldist=TRUE,keep.rawdist=FALSE,seed=NULL,cluster=1,type=NULL,dispatch=NULL,marg.null=NULL,marg.par=NULL,keep.margpar=TRUE,ncp=NULL,perm.mat=NULL,keep.index=FALSE,keep.label=FALSE){
  ##sanity checks / formatting
  #X
  if(missing(X)) stop("Argument X is missing")
  if(inherits(X,"eSet")){ 
    if(is.character(Y)) Y<-pData(X)[,Y]
    if(is.character(Z)){
      if(Z%in%Y){
        Z<-Z[!(Z%in%Y)]
	warning(paste("Outcome Y=",Y,"should not be included in the covariates Z=",Z,". Removing Y from Z",sep=""))
	}
      Z<-pData(X)[,Z]
    }
    X<-exprs(X)
  }
  X<-as.matrix(X)
  dx<-dim(X)
  if(length(dx)==0) stop("dim(X) must have positive length")
  p<-dx[1]
  n<-dx[2]
  #W
  if(!is.null(W)){
    W[W<=0]<-NA
    if(is.vector(W) & length(W)==n) W <- matrix(rep(W,p),nrow=p,ncol=n,byrow=TRUE)
    if(is.vector(W) & length(W)==p) W <- matrix(rep(W,n),nrow=p,ncol=n)
    if(test%in%c("f","f.block","f.twoway","t.cor","z.cor")){
      warning("Weights can not be used with F-tests or tests of correlation parameters, arg W is being ignored.")
      W<-NULL
    }
  }
  #Y
  if(!is.null(Y)){
    if(is.Surv(Y)){
      if(test!="coxph.YvsXZ") stop(paste("Test ",test," does not work with a survival object Y",sep=""))
    }
    else{
      Y<-as.matrix(Y)
      if(ncol(Y)!=1) stop("Argument Y must be a vector")
    }
    if(nrow(Y)!=n) stop("Outcome Y has length ",nrow(Y),", not equal to n=",n)
  }
  if(test=="t.pair") n <- dx[2]/2
  #Z
  if(!is.null(Z)){
    Z<-as.matrix(Z)
    if(nrow(Z)!=n) stop("Covariates in Z have length ",nrow(Z),", not equal to n=",n,"\n")
    #Z.incl tells which columns of Z to include in model
    if(is.null(Z.incl)) Z.incl<-(1:ncol(Z))
    if(length(Z.incl)>ncol(Z)) stop("Number of columns in Z.incl ",length(Z.incl)," exceeds ncol(Z)=",ncol(Z))
    if(is.logical(Z.incl)) Z.incl<-(1:ncol(Z))[Z.incl]
    if(is.character(Z.incl) & length(Z.incl)!=sum(Z.incl%in%colnames(Z))) stop(paste("Z.incl=",Z.incl," names columns not in Z",sep=""))
    Za<-Z[,Z.incl]
    #Z.test tells which column of Z to test for an association
    if(test=="lm.XvsZ"){
      if(is.null(Z.test)){
        warning(paste("Z.test not specified, testing for association with variable in first column of Z:",colnames(Z)[1],sep=""))
	Z.test<-1
      }
      if(is.logical(Z.test)) Z.test<-(1:ncol(Z))[Z.test]
      if(is.character(Z.test) & !(Z.test%in%colnames(Z))) stop(paste("Z.test=",Z.test," names a column not in Z",sep=""))
      if(is.numeric(Z.test) & !(Z.test%in%(1:ncol(Z)))) stop("Value of Z.test must be >0 and <",ncol(Z))
      if(Z.test%in%Z.incl){
        Z.incl<-Z.incl[!(Z.incl%in%Z.test)]
	Za<-Z[,Z.incl]
      }
      Za<-cbind(Z[,Z.test],Za)
    }
    Z<-Za
    rm(Za)
  }
  #test
  TESTS<-c("t.onesamp","t.twosamp.equalvar","t.twosamp.unequalvar","t.pair","f","f.block","f.twoway","lm.XvsZ","lm.YvsXZ","coxph.YvsXZ","t.cor","z.cor")
  test<-TESTS[pmatch(test,TESTS)]
  if(is.na(test)) stop(paste("Invalid test, try one of ",TESTS,sep=""))
  #robust + see below with choice of nulldist
  if(test=="coxph.YvsXZ" & robust==TRUE)
    warning("No robust version of coxph.YvsXZ, proceding with usual version")
  #temp until fix
  if((test=="t.onesamp" | test=="t.pair") & robust==TRUE)
    stop("Robust test statistics currently not available for one-sample or two-sample paired test statistics.")
  #alternative
  ALTS<-c("two.sided","less","greater")
  alternative<-ALTS[pmatch(alternative,ALTS)]
  if(is.na(alternative)) stop(paste("Invalid alternative, try one of ",ALTS,sep=""))
  #null values
  if(length(psi0)>1) stop(paste("In current implementation, all hypotheses must have the same null value. Number of null values: ",length(psi0),">1",sep=""))
  #Error rate
  ERROR<-c("fwer","gfwer","tppfp","fdr")
  typeone<-ERROR[pmatch(typeone,ERROR)]
  if(is.na(typeone)) stop(paste("Invalid typeone, try one of ",ERROR,sep=""))
  if(any(alpha<0) | any(alpha>1)) stop("Nominal level alpha must be between 0 and 1")
  nalpha<-length(alpha)
  reject<-
    if(nalpha) array(dim=c(p,nalpha),dimnames=list(rownames(X),paste("alpha=",alpha,sep="")))
    if(test=="z.cor" | test=="t.cor") matrix(nrow=0,ncol=0) # deprecated for correlations, rownames now represent p choose 2 edges - too weird and clunky in current state for output.
    else matrix(nrow=0,ncol=0)
  if(typeone=="gfwer"){
    if(get.cr==TRUE) warning("Confidence regions not currently implemented for gFWER")
    if(get.cutoff==TRUE) warning("Cut-offs not currently implemented for gFWER")
    get.cr<-get.cutoff<-FALSE
    if(k<0) stop("Number of false positives can not be negative")
    if(k>=p) stop(paste("Number of false positives must be less than number of tests=",p,sep=""))
    if(length(k)>1){
      k<-k[1]
      warning("can only compute gfwer(k) adjp for one value of k at a time (using first value), try fwer2gfwer() function for multiple k")
    }
  }
  if(typeone=="tppfp"){
    if(get.cr==TRUE) warning("Confidence regions not currently implemented for TPPFP")
    if(get.cutoff==TRUE) warning("Cut-offs not currently implemented for TPPFP")
    get.cr<-get.cutoff<-FALSE
    if(q<0) stop("Proportion of false positives, q, can not be negative")
    if(q>1) stop("Proportion of false positives, q, must be less than 1")
    if(length(q)>1){
      q<-q[1]
      warning("Can only compute tppfp adjp for one value of q at a time (using first value), try fwer2tppfp() function for multiple q")
    }
  }
  if(typeone=="fdr"){
    if(!nalpha) stop("Must specify a nominal level alpha for control of FDR")
    if(get.cr==TRUE) warning("Confidence regions not currently implemented for FDR")
    if(get.cutoff==TRUE) warning("Cut-offs not currently implemented for FDR")
    get.cr<-get.cutoff<-FALSE
  }		
  #null distribution
  NULLS<-c("boot","boot.cs","boot.ctr","boot.qt","ic","perm")
  nulldist<-NULLS[pmatch(nulldist,NULLS)]
  if(is.na(nulldist)) stop(paste("Invalid nulldist, try one of ",NULLS,sep=""))
  if(nulldist=="boot"){
    nulldist <- "boot.cs"
    warning("nulldist='boot' is deprecated and now corresponds to 'boot.cs'. Proceeding with default center and scaled null distribution.")
  }
  if(nulldist!="perm" & test=="f.block") stop("f.block test only available with permutation null distribution. Try test=f.twoway")
  if((nulldist=="perm" | nulldist=="ic") & keep.rawdist==TRUE) stop("Test statistics distribution estimation using keep.rawdist=TRUE is only available with a bootstrap-based null distribution")
  if(nulldist=="boot.qt" & robust==TRUE) stop("Quantile transform method requires parametric marginal nulldist.  Set robust=FALSE")
  if(nulldist=="boot.qt" & standardize==FALSE) stop("Quantile transform method requires standardized test statistics.  Set standardize=TRUE")
  if(nulldist=="ic" & robust==TRUE) stop("Influence curve null distributions available only for (parametric) t-statistics.  Set robust=FALSE")
  if(nulldist=="ic" & standardize==FALSE) stop("Influence curve null distributions available only for (standardized) t-statistics.  Set standardize=TRUE")
  if(nulldist=="ic" & (test=="f" | test=="f.twoway" | test=="f.block" | test=="coxph.YvsXZ")) stop("Influence curve null distributions available only for tests of mean, regression and correlation parameters. Cox PH also not yet implemented.")
  if(nulldist!="ic" & (test=="t.cor" | test=="z.cor")) stop("Tests of correlation parameters currently only implemented for influence curve null distributions")
  if((test!="t.cor" & test!="z.cor") & keep.index) warning("Matrix of indices only returned for tests of correlation parameters")
  ### specifically for sampling null test statistics with IC nulldist
  MVNS <- c("mvrnorm","Cholesky")
  MVN.method <- MVNS[pmatch(MVN.method,MVNS)]
  if(is.na(MVN.method)) stop("Invalid sampling method for IC-based MVN null test statistics.  Try either 'mvrnorm' or 'Cholesky'")
  #methods
  METHODS<-c("ss.maxT","ss.minP","sd.maxT","sd.minP")
  method<-METHODS[pmatch(method,METHODS)]
  if(is.na(method)) stop(paste("Invalid method, try one of ",METHODS,sep=""))
  #estimate and conf.reg
  ftest<-FALSE
  if(test=="f" | test=="f.block"){
    ftest<-TRUE
    if(get.cr) stop("Confidence intervals not available for F tests, try get.cr=FALSE")
    if(!is.null(W)) warning("Weighted F tests not yet implemented, proceding with unweighted version")
  }

  
  #permutation null distribution - self contained in this if statement
  if(nulldist=="perm"){
    if(method=="ss.minP" | method=="ss.maxT") stop("Only step-down procedures are currently available with permutation nulldist")
    if(smooth.null) warning("Kernal density p-values not available with permutation nulldist")
    if(get.cr) warning("Confidence regions not available with permutation nulldist")
    if(get.cutoff) warning("Cut-offs not available with permutation nulldist")
    #if(keep.nulldist) warning("keep.nulldist not available with permutation nulldist")
    ptest<-switch(test,
                  t.onesamp=stop("One sample t-test not available with permutation nulldist"),
                  t.twosamp.equalvar=ifelse(robust,"wilcoxon","t.equalvar"),
                  t.twosamp.unequalvar="t",
                  t.pair="pairt",
                  f="f",
                  f.block="blockf",
                  f.twoway=stop("f.twoway not available with permutation nulldist"),
                  lm.XvsZ=stop("lm.XvsZ not available with permutation nulldist"),
                  lm.YvsXZ=stop("lm.YvsXZ not available with permutation nulldist"),
                  coxph.YvsXZ=stop("coxph.YvsXZ not available with permutation nulldist"),
                  t.cor=stop("t.cor not available with permutation nulldist"),
                  z.cor=stop("z.cor not available with permutation nulldist")
                  )
    pside<-switch(alternative,two.sided="abs",less="lower",greater="upper")
    pnonpara<-
      if(robust)"y"
      else "n"
    if(any(is.na(Y))){
      bad<-is.na(Y)
      Y<-Y[!bad]
      X<-X[,!bad]
      warning("No NAs allowed in Y, these observations have been removed.")
    }
    presult<-switch(method,
                    sd.maxT=mt.maxT(X,classlabel=Y,test=ptest,side=pside,B=B,nonpara=pnonpara),
                    sd.minP=mt.minP(X,classlabel=Y,test=ptest,side=pside,B=B,nonpara=pnonpara)
                    )
    if(typeone=="fwer" & nalpha){
      for(a in 1:nalpha) reject[,a]<-(presult$adjp<=alpha[a])
    }
    if(typeone=="gfwer"){
      presult$adjp<-fwer2gfwer(presult$adjp,k)
      if(nalpha){
        for(a in 1:nalpha) reject[,a]<-(presult$adjp<=alpha[a])
      }
      if(!get.adjp)
        presult$adjp<-vector("numeric",0)
    }
    if(typeone=="tppfp"){
      presult$adjp<-fwer2tppfp(presult$adjp,q)
      if(nalpha){
        for(a in 1:nalpha) reject[,a]<-(presult$adjp<=alpha[a])
      }
      if(!get.adjp)
        presult$adjp<-vector("numeric",0)
    }
    if(typeone=="fdr"){
      temp<-fwer2fdr(presult$adjp,fdr.method,alpha)
      reject<-temp$reject
      if(!get.adjp) presult$adjp<-vector("numeric",0)
      else presult$adjp<-temp$adjp
      rm(temp)
    }			
    #output results
    orig<-order(presult$index)
    if(keep.label) label <- as.numeric(Y)
    else label <- vector("numeric",0)
    out<-new("MTP",statistic=presult$teststat[orig],estimate=vector("numeric",0),sampsize=n,rawp=presult$rawp[orig],adjp=presult$adjp[orig],conf.reg=array(dim=c(0,0,0)),cutoff=matrix(nrow=0,ncol=0),reject=as.matrix(reject[orig,]),rawdist=matrix(nrow=0,ncol=0),nulldist=matrix(nrow=0,ncol=0),nulldist.type="perm",marg.null=vector("character",0),marg.par=matrix(nrow=0,ncol=0),label=label,index=matrix(nrow=0,ncol=0),call=match.call(),seed=vector("integer",0))
  }
  
  else{ # This should apply to all other MTP calls using the bootstrap and IC nulldists.
    if(nulldist=="boot.qt"){ # get parameter vals for quantile transform.
      # Get parameter values for the quantile transformed nulldist
      if(!is.null(marg.par)){
        if(is.matrix(marg.par)) marg.par <- marg.par
        if(is.vector(marg.par)) marg.par <- matrix(rep(marg.par,p),nrow=p,ncol=length(marg.par),byrow=TRUE)
        }
      if(is.null(ncp)) ncp = 0
      if(!is.null(perm.mat)){ 
        if(dim(X)[1]!=dim(perm.mat)[1]) stop("perm.mat must same number of rows as X.")
        }
    
      nstats <- c("t.twosamp.unequalvar","z.cor","lm.XvsZ","lm.YvsXZ","coxph.lmYvsXZ")
      tstats <- c("t.onesamp","t.twosamp.equalvar","t.pair","t.cor")
      fstats <- c("f","f.block","f.twoway")
      
      # If default , set values of marg.null to pass on.
      if(is.null(marg.null)){
	  if(any(nstats == test)) marg.null="normal"
	  if(any(tstats == test)) marg.null="t"
	  if(any(fstats == test)) marg.null="f"
        }
      else{ # Check to see that user-supplied entries make sense.  
        MARGS <- c("normal","t","f","perm")
        marg.null <- MARGS[pmatch(marg.null,MARGS)]
        if(is.na(marg.null)) stop("Invalid marginal null distribution. Try one of: normal, t, f, or perm")
        if(any(tstats==test) & marg.null == "f") stop("Choice of test stat and marginal nulldist do not match")
        if(any(fstats==test) & (marg.null == "normal" | marg.null=="t")) stop("Choice of test stat and marginal nulldist do not match")
        if(marg.null=="perm" & is.null(perm.mat)) stop("Must supply a matrix of permutation test statistics if marg.null='perm'")
        if(marg.null=="f" & ncp < 0) stop("Cannot have negative noncentrality parameter with F distribution.")
      }
    
      # If default (=NULL), set values of marg.par. Return as m by 1 or 2 matrix.
      if(is.null(marg.par)){
		marg.par <- switch(test,
                          t.onesamp = n-1,
                          t.twosamp.equalvar = n-2,
                          t.twosamp.unequalvar = c(0,1),
                          t.pair = floor(n/2-1),
                          f = c(length(is.finite(unique(Y)))-1,dim(X)[2]- length(is.finite(unique(Y))) ),
                          f.twoway = {
                            c(length(is.finite(unique(Y)))-1, dim(X)[2]-(length(is.finite(unique(Y)))*length(gregexpr('12', paste(Y, collapse=""))[[1]]))-2)
                            },
                          lm.XvsZ = c(0,1),
                          lm.YvsXZ = c(0,1),
                          coxph.YvsXZ = c(0,1),
                          t.cor = n-2,
                          z.cor = c(0,1)
                          )
      marg.par <- matrix(rep(marg.par,dim(X)[1]),nrow=dim(X)[1],ncol=length(marg.par),byrow=TRUE)
              }
     else{ # Check that user-supplied values of marg.par make sense (marg.par != NULL)
       if((marg.null=="t" | marg.null=="f") & any(marg.par[,1]==0)) stop("Cannot have zero df with t or F distributions. Check marg.par settings")
       if(marg.null=="t" & dim(marg.par)[2]>1) stop("Too many parameters for t distribution.  marg.par should have length 1.")
       if((marg.null=="f" | marg.null=="normal") & dim(marg.par)[2]!=2) stop("Incorrect number of parameters defining marginal null distribution.  marg.par should have length 2.")
     }
    }

    ##making a closure for the particular test
    theta0<-0
    tau0<-1
    stat.closure<-switch(test,
                         t.onesamp=meanX(psi0,na.rm,standardize,alternative,robust),
                         t.twosamp.equalvar=diffmeanX(Y,psi0,var.equal=TRUE,na.rm,standardize,alternative,robust),
                         t.twosamp.unequalvar=diffmeanX(Y,psi0,var.equal=FALSE,na.rm,standardize,alternative,robust),
                         t.pair={
                           uY<-sort(unique(Y))
                           if(length(uY)!=2) stop("Must have two class labels for this test")
                           if(trunc(ncol(X)/2)!=ncol(X)/2) stop("Must have an even number of samples for this test")
                           X<-X[,Y==uY[2]]-X[,Y==uY[1]]
                           Y<-NULL
                           n<-dim(X)[2]
                           meanX(psi0,na.rm,standardize,alternative,robust)
                         },
                         f={
                           theta0<-1
                           tau0<-2/(length(unique(Y))-1)
                           FX(Y,na.rm,robust)
                         },
                         f.twoway={
                           theta0<-1
                           tau0 <- 2/((length(unique(Y))*length(gregexpr('12', paste(Y, collapse=""))[[1]]))-1)
                           twowayFX(Y,na.rm,robust)
                         },
                         lm.XvsZ=lmX(Z,n,psi0,na.rm,standardize,alternative,robust),
                         lm.YvsXZ=lmY(Y,Z,n,psi0,na.rm,standardize,alternative,robust),
                         coxph.YvsXZ=coxY(Y,Z,psi0,na.rm,standardize,alternative),
                         t.cor=NULL,
                         z.cor=NULL)
    ##computing observed test statistics
    if(test=="t.cor" | test=="z.cor") obs<-corr.Tn(X,test=test,alternative=alternative,use="pairwise")
    else obs<-get.Tn(X,stat.closure,W)
    ##or computing influence curves
    if(nulldist=="ic"){
      rawdistn <- matrix(nrow=0,ncol=0)
      nulldistn<-switch(test,
                        t.onesamp=corr.null(X,W,Y,Z,test="t.onesamp",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        t.pair=corr.null(X,W,Y,Z,test="t.pair",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        t.twosamp.equalvar=corr.null(X,W,Y,Z,test="t.twosamp.equalvar",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        t.twosamp.unequalvar=corr.null(X,W,Y,Z,test="t.twosamp.unequalvar",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        lm.XvsZ=corr.null(X,W,Y,Z,test="lm.XvsZ",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        lm.YvsXZ=corr.null(X,W,Y,Z,test="lm.YvsXZ",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        t.cor=corr.null(X,W,Y,Z,test="t.cor",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat),
                        z.cor=corr.null(X,W,Y,Z,test="z.cor",alternative,use="pairwise",B,MVN.method,penalty,ic.quant.trans,marg.null,marg.par,perm.mat)
                        )
    }

    ## Cluster Checking
    if ((!is.numeric(cluster))&(!inherits(cluster,c("MPIcluster", "PVMcluster", "SOCKcluster"))))
       stop("Cluster argument must be integer or cluster object")
    ## Create cluster if cluster > 1 and load required packages on nodes
    if(is.numeric(cluster)){
      if(cluster>1){
    ## Check installation of packages
      have_snow <- qRequire("snow")
      if(!have_snow) stop("The package snow is required to use a cluster. Either snow is not installed or it is not in the standard library location.")
      if (is.null(type))
         stop("Must specify type argument to use a cluster. Alternatively, provide a cluster object as the argument to cluster.")
      if (type=="SOCK")
         stop("Create desired cluster and specify cluster object as the argument to cluster directly.")
      if ((type!="PVM")&(type!="MPI"))
         stop("Type must be MPI or PVM")
      else if (type=="MPI"){
         have_rmpi <- qRequire("Rmpi")
         if(!have_rmpi) stop("The package Rmpi is required for the specified type. Either Rmpi is not installed or it is not in the standard library location.")
      }
      else if (type=="PVM"){
         have_rpvm <- qRequire("rpvm")
         if(!have_rpvm) stop("The package rpvm is required for the specified type. Either rpvm is not installed or it is not in the standard library location.")
      }
      cluster <- makeCluster(cluster, type)
      clusterEvalQ(cluster, {library(Biobase); library(multtest)})
      if (is.null(dispatch)) dispatch=0.05
      }
    }
    else if(inherits(cluster,c("MPIcluster", "PVMcluster", "SOCKcluster"))){
      clusterEvalQ(cluster, {library(Biobase); library(multtest)})
      if (is.null(dispatch)) dispatch=0.05
    }

    ##computing the nonparametric bootstrap (null) distribution
    if(nulldist=="boot.cs" | nulldist=="boot.ctr" | nulldist=="boot.qt"){
      nulldistn<-boot.null(X,Y,stat.closure,W,B,test,nulldist,theta0,tau0,marg.null,marg.par,ncp,perm.mat,alternative,seed,cluster,dispatch,keep.nulldist,keep.rawdist)
     if(inherits(cluster,c("MPIcluster", "PVMcluster", "SOCKcluster")))  stopCluster(cluster)
    rawdistn <- nulldistn$rawboot
    nulldistn <- nulldistn$muboot
    }

    
    ##performing multiple testing
    #rawp values
    rawp<-apply((obs[1,]/obs[2,])<=nulldistn,1,mean)
    if(smooth.null & (min(rawp,na.rm=TRUE)==0)){
      zeros<-(rawp==0)
      if(sum(zeros)==1){
        den<-density(nulldistn[zeros,],to=max(obs[1,zeros]/obs[2,zeros],nulldist[zeros,],na.rm=TRUE),na.rm=TRUE)
	rawp[zeros]<-sum(den$y[den$x>=(obs[1,zeros]/obs[2,zeros])])/sum(den$y)
      }
      else{
        den<-apply(nulldistn[zeros,],1,density,to=max(obs[1,zeros]/obs[2,zeros],nulldistn[zeros,],na.rm=TRUE),na.rm=TRUE)
	newp<-NULL
	stats<-obs[1,zeros]/obs[2,zeros]
	for(i in 1:length(den)){
          newp[i]<-sum(den[[i]]$y[den[[i]]$x>=stats[i]])/sum(den[[i]]$y)
	}
        rawp[zeros]<-newp		
      }
      rawp[rawp<0]<-0
    }
    #c, cr, adjp values
    pind<-ifelse(typeone!="fwer",TRUE,get.adjp)
    if(method=="ss.maxT") out<-ss.maxT(nulldistn,obs,alternative,get.cutoff,get.cr,pind,alpha)
    if(method=="ss.minP") out<-ss.minP(nulldistn,obs,rawp,alternative,get.cutoff,get.cr,pind,alpha)
    if(method=="sd.maxT") out<-sd.maxT(nulldistn,obs,alternative,get.cutoff,get.cr,pind,alpha)
    if(method=="sd.minP") out<-sd.minP(nulldistn,obs,rawp,alternative,get.cutoff,get.cr,pind,alpha)
    if(typeone=="fwer" & nalpha & (test!="t.cor" & test !="z.cor")){
      for(a in 1:nalpha) reject[,a]<-(out$adjp<=alpha[a])
    }
    #augmentation procedures
    if(typeone=="gfwer"){
      out$adjp<-as.numeric(fwer2gfwer(out$adjp,k))
      out$c<-matrix(nrow=0,ncol=0)
      out$cr<-array(dim=c(0,0,0))
      if(nalpha){
        for(a in 1:nalpha) reject[,a]<-(out$adjp<=alpha[a])
      }
      if(!get.adjp) out$adjp<-vector("numeric",0)
    }
    if(typeone=="tppfp"){
      out$adjp<-as.numeric(fwer2tppfp(out$adjp,q))
      out$c<-matrix(nrow=0,ncol=0)
      out$cr<-array(dim=c(0,0,0))
      if(nalpha){
        for(a in 1:nalpha) reject[,a]<-(out$adjp<=alpha[a])
      }
      if(!get.adjp) out$adjp<-vector("numeric",0)
    }
    if(typeone=="fdr"){
      out$c<-matrix(nrow=0,ncol=0)
      out$cr<-array(dim=c(0,0,0))
      temp<-fwer2fdr(out$adjp,fdr.method,alpha)
      reject<-temp$reject
      if(!get.adjp) out$adjp<-vector("numeric",0)
      else out$adjp<-temp$adjp
      rm(temp)
    }
    #output results
    if(!keep.nulldist) nulldistn<-matrix(nrow=0,ncol=0)
    if(!keep.rawdist) rawdist<-matrix(nrow=0,ncol=0)
    if(nulldist!="boot.qt"){  
      marg.null <- vector("character")
      marg.par <- matrix(nrow=0,ncol=0)
    }
    if(!keep.label) label <- vector("numeric",0)
    if(!keep.index) index <- matrix(nrow=0,ncol=0)
    if(test!="z.cor" & test !="t.cor") index <- matrix(nrow=0,ncol=0)
    if(keep.index & (test!="z.cor" | test !="t.cor")){
      index <- t(combn(p,2))
      colnames(index) <- c("Var1","Var2")
    }
    names(out$adjp)<-names(rawp)
    estimates <- obs[3,]*obs[1,]
    if(ftest) estimates <- vector("numeric",0)
    if(test=="t.onesamp" | test=="t.pair") estimates <- obs[3,]*obs[1,]/sqrt(n)
    out<-new("MTP",statistic=(obs[3,]*obs[1,]/obs[2,]),
      estimate=estimates,
      sampsize=n,rawp=rawp,adjp=out$adjp,conf.reg=out$cr,cutoff=out$c,reject=reject,
      rawdist=rawdistn,nulldist=nulldistn,nulldist.type=nulldist,
      marg.null=marg.null,marg.par=marg.par,label=label,index=index,
      call=match.call(),seed=as.integer(seed))
  }
  return(out)
}

#funtions to compute cutoffs and adjusted pvals

ss.maxT<-function(null,obs,alternative,get.cutoff,get.cr,get.adjp,alpha=0.05){
  p<-dim(null)[1]
  B<-dim(null)[2]
  nalpha<-length(alpha)
  mT<-apply(null,2,max)
  getc<-matrix(nrow=0,ncol=0)
  getcr<-array(dim=c(0,0,0))
  getp<-vector(mode="numeric")
  if(get.cutoff | get.cr){
    getc<-array(dim=c(p,nalpha),dimnames=list(dimnames(null)[[1]],paste("alpha=",alpha,sep="")))
    if(get.cr) getcr<-array(dim=c(p,2,nalpha),dimnames=list(dimnames(null)[[1]],c("LB","UB"),paste("alpha=",alpha,sep="")))
    for(a in 1:nalpha){
      getc[,a]<-rep(quantile(mT,pr=(1-alpha[a])),p)
      if(get.cr) getcr[,,a]<-cbind(ifelse(rep(alternative=="less",p),rep(-Inf,p),obs[3,]*obs[1,]-getc[,a]*obs[2,]),ifelse(rep(alternative=="greater",p),rep(Inf,p),obs[3,]*obs[1,]+getc[,a]*obs[2,]))
    }
  }
  if(get.adjp) getp<-apply((obs[1,]/obs[2,])<=matrix(mT,nrow=p,ncol=B,byrow=TRUE),1,mean)
  if(!get.cutoff) getc<-matrix(nrow=0,ncol=0)
  list(c=getc,cr=getcr,adjp=getp)
}

ss.minP<-function(null,obs,rawp,alternative,get.cutoff,get.cr,get.adjp,alpha=0.05){
  p<-dim(null)[1]
  B<-dim(null)[2]
  nalpha<-length(alpha)
  getc<-matrix(nrow=0,ncol=0)
  getcr<-array(dim=c(0,0,0))
  getp<-vector(mode="numeric")
  R<-apply(null,1,rank)
  if(get.cutoff | get.cr){
    getc<-array(dim=c(p,nalpha),dimnames=list(dimnames(null)[[1]],paste("alpha=",alpha,sep="")))
    if(get.cr) getcr<-array(dim=c(p,2,nalpha),dimnames=list(dimnames(null)[[1]],c("LB","UB"),paste("alpha=",alpha,sep="")))
    for(a in 1:nalpha){
      q<-quantile(apply(R,1,max),1-alpha[a])
      for(j in 1:p){
        getc[j,a]<-min(c(null[j,R[,j]>=q],max(null[j,])))
      }
      if(get.cr) getcr[,,a]<-cbind(ifelse(rep(alternative=="less",p),rep(-Inf,p),obs[3,]*obs[1,]-getc[,a]*obs[2,]),ifelse(rep(alternative=="greater",p),rep(Inf,p),obs[3,]*obs[1,]+getc[,a]*obs[2,]))
    }
  }
  if(get.adjp){
    R<-matrix(apply((B+1-R)/B,1,min),nrow=p,ncol=B,byrow=TRUE)
    getp<-apply(rawp>=R,1,mean)
  }
  if(!get.cutoff) getc<-matrix(nrow=0,ncol=0)
  list(c=getc,cr=getcr,adjp=getp)
}

sd.maxT<-function(null,obs,alternative,get.cutoff,get.cr,get.adjp,alpha=0.05){
  p<-dim(null)[1]
  B<-dim(null)[2]
  nalpha<-length(alpha)
  ord<-rev(order(obs[1,]/obs[2,]))
  mT<-null[ord[p],]
  getc<-matrix(nrow=0,ncol=0)
  getcr<-array(dim=c(0,0,0))
  getp<-vector(mode="numeric")
  if(get.cutoff | get.cr){
    getc<-array(dim=c(p,nalpha),dimnames=list(dimnames(null)[[1]],paste("alpha=",alpha,sep="")))
    for(a in 1:nalpha) getc[ord[p],a]<-quantile(mT,pr=1-alpha[a])
  }
  if(get.adjp) getp[ord[p]]<-mean((obs[1,]/obs[2,])[ord[p]]<=mT)
  for(j in (p-1):1){
    mT<-pmax(mT,null[ord[j],])
    if(get.adjp) getp[ord[j]]<-mean((obs[1,ord[j]]/obs[2,ord[j]])<=mT)
    if(get.cutoff | get.cr){
      for(a in 1:nalpha) getc[ord[j],a]<-quantile(mT,pr=(1-alpha[a]))
    }
  }
  c.ind<-rep(TRUE,nalpha)
  for(j in 2:p){
    if(get.adjp) getp[ord[j]]<-max(getp[ord[j]],getp[ord[j-1]])
    if(get.cutoff | get.cr){
      for(a in 1:nalpha){
        if(c.ind[a]){
          if((obs[1,]/obs[2,])[ord[j-1]]<=getc[ord[j-1],a]){
            getc[ord[j:p],a]<-Inf
            c.ind[a]<-FALSE
          }
	}
      }
    }
  }
  if(get.cr){
    getcr<-array(dim=c(p,2,nalpha),dimnames=list(dimnames(null)[[1]],c("LB","UB"),paste("alpha=",alpha,sep="")))
    for(a in 1:nalpha){
      getcr[,,a]<-cbind(ifelse(rep(alternative=="less",p),rep(-Inf,p),obs[3,]*obs[1,]-getc[,a]*obs[2,]),ifelse(rep(alternative=="greater",p),rep(Inf,p),obs[3,]*obs[1,]+getc[,a]*obs[2,]))
    }
  }
  if(!get.cutoff) getc<-matrix(nrow=0,ncol=0)
  list(c=getc,cr=getcr,adjp=getp)
}

sd.minP<-function(null,obs,rawp,alternative,get.cutoff,get.cr,get.adjp,alpha=0.05){
  p<-dim(null)[1]
  B<-dim(null)[2]
  nalpha<-length(alpha)
  ord<-order(rawp)
  R<-apply(null,1,rank) #B x p
  mR<-R[,ord[p]]
  getc<-matrix(nrow=0,ncol=0)
  getcr<-array(dim=c(0,0,0))
  getp<-vector(mode="numeric")
  if(get.cutoff | get.cr){
    getc<-array(dim=c(p,nalpha),dimnames=list(dimnames(null)[[1]],paste("alpha=",alpha,sep="")))
    for(a in 1:nalpha){
      q<-quantile(mR,pr=1-alpha[a])
      getc[ord[p],a]<-min(c(null[ord[p],R[,ord[p]]>=q],max(null[ord[p],])))
    }
  }
  if(get.adjp){
    mP<-(B+1-mR)/B
    getp[ord[p]]<-mean(rawp[ord[p]]>=mP)
  }
  for(j in (p-1):1){
    mR<-pmax(mR,R[,ord[j]])
    if(get.adjp){
      mP<-(B+1-mR)/B
      getp[ord[j]]<-mean(rawp[ord[j]]>=mP)
    }
    if(get.cutoff | get.cr){
      for(a in 1:nalpha){
        q<-quantile(mR,pr=1-alpha[a])
	getc[ord[j],a]<-min(c(null[ord[j],R[,ord[j]]>=q],max(null[ord[j],])))
      }
    }
  }
  c.ind<-rep(TRUE,nalpha)
  for(j in 2:p){
    if(get.adjp) getp[ord[j]]<-max(getp[ord[j]],getp[ord[j-1]])
    if(get.cutoff | get.cr){
      for(a in 1:nalpha){
        if(c.ind[a]){
          if((obs[1,]/obs[2,])[ord[j-1]]<=getc[ord[j-1],a]){
            getc[ord[j:p],a]<-Inf
            c.ind[a]<-FALSE
          }
	}
      }
    }
  }
  if(get.cr){
    getcr<-array(dim=c(p,2,nalpha),dimnames=list(dimnames(null)[[1]],c("LB","UB"),paste("alpha=",alpha,sep="")))
    for(a in 1:nalpha){
      getcr[,,a]<-cbind(ifelse(rep(alternative=="less",p),rep(-Inf,p),obs[3,]*obs[1,]-getc[,a]*obs[2,]),ifelse(rep(alternative=="greater",p),rep(Inf,p),obs[3,]*obs[1,]+getc[,a]*obs[2,]))
    }
  }
  if(!get.cutoff) getc<-matrix(nrow=0,ncol=0)
  list(c=getc,cr=getcr,adjp=getp)
}

#functions to convert FWER adjp to AMTP (gFWER, TPPFP) adjp:
fwer2gfwer<-function(adjp,k=0){
  ord<-order(adjp)
  m<-length(adjp)
  if(any(k>=m)) stop(paste("number of rejections k=",k," must be less than number of hypotheses=",m,sep=""))
  newp<-NULL
  for(j in k) newp<-rbind(newp,c(rep(0,j),adjp[ord[1:(m-j)]]))
  rownames(newp)<-k
  colnames(newp)<-ord
  newp<-matrix(newp[,order(ord)],ncol=m,byrow=FALSE)
  return(t(newp))
}

fwer2tppfp<-function(adjp,q=0.05){
  ord<-order(adjp)
  m<-length(adjp)
  newp<-NULL
  if(any(q>1)|any(q<0)) stop(paste("proportion of false positives q=",q," must be in [0,1]",sep=""))
  for(l in q) newp<-rbind(newp,adjp[ord][ceiling((1:m)*(1-l))])
  rownames(newp)<-q
  colnames(newp)<-names(ord)
  newp<-matrix(newp[,order(ord)],ncol=m,byrow=FALSE)
  return(t(newp))
}

#function to compute rejection indicator for FDR methods
fwer2fdr<-function(adjp,method="both",alpha=0.05){
  get.cons<-function(adjp,alpha,ord,M,nalpha){
    newp<-NULL
    for(m in 1:M){
      #try ceiling/floor
      k<-if(m%%2) 0:((m-1)/2) else 0:(m/2)
      f<-2*adjp[ord][m-k]
      u<-2*(k+1)/m
      l<-2*k/m
      if(sum(f<=u)){
        ind<-min(which(f<=u))
	newp[ord[m]]<-
	if(f[ind]>=l[ind]) f[ind]
	else l[ind]
      }
      else newp[ord[m]]<-1
    }
    newp[newp>1]<-1
    a<-alpha/2
    rejections<-matrix(nrow=M,ncol=nalpha)
    for(i in 1:nalpha) rejections[,i]<-(fwer2tppfp(adjp,a[i])<=a[i])
    return(list(reject=rejections,adjp=newp))
  }
  get.restr<-function(adjp,alpha,ord,M,nalpha){
    newp<-NULL
    ginv<-function(x) 1-(1-x)^2
    for(m in 1:M){
      k<-m:1
      f<-adjp[ord][k]
      u<-1-(k-1)/m
      l<-1-k/m
      if(sum(f<=u)){
        ind<-min(which(f<=u))
	newp[ord[m]]<-
          if(f[ind]>=l[ind]) ginv(f[ind])
          else ginv(l[ind])
      }
      else
	newp[ord[m]]<-1
    }
    newp[newp>1]<-1
    a<-1-sqrt(1-alpha)
    rejections<-matrix(nrow=M,ncol=nalpha)
    for(i in 1:nalpha) rejections[,i]<-(fwer2tppfp(adjp,a[i])<=a[i])
    return(list(reject=rejections,adjp=newp))
  }
  ord<-order(adjp)
  nalpha<-length(alpha)
  M<-length(adjp)
  if(method=="both"){
    rejections<-array(dim=c(M,nalpha,2),dimnames=list(NULL,paste("alpha=",alpha,sep=""),c("conservative","restricted")))
    newp<-matrix(nrow=M,ncol=2,dimnames=list(NULL,c("conservative","restricted")))
    temp<-get.cons(adjp,alpha,ord,M,nalpha)
    rejections[,,"conservative"]<-temp$reject
    newp[,"conservative"]<-temp$adjp
    temp<-get.restr(adjp,alpha,ord,M,nalpha)
    rejections[,,"restricted"]<-temp$reject
    newp[,"restricted"]<-temp$adjp
    rm(temp)
  }
  else{
    rejections<-matrix(nrow=M,ncol=nalpha,dimnames=list(NULL,paste("alpha=",alpha,sep="")))
    newp<-NULL
    if(method=="conservative") temp<-get.cons(adjp,alpha,ord,M,nalpha)
    else temp<-get.restr(adjp,alpha,ord,M,nalpha)
    rejections<-temp$reject
    newp<-temp$adjp
    rm(temp)
  }
  return(list(reject=rejections,adjp=newp))
}