samr <- function(data, params, resp.type=c("twoclass","onesample","paired","survival","categorical","quantitative","timearea"), s0=NULL, s0.perc=NULL, nperm=50, permy=NULL){

#SAM method. copyright june 2000: Goss, Tibshirani and Chu.
# coded by r tibshirani; 
# y is response measure: 1,2 for two twoclass groups,
# y=1,1,1 for onesample problem,
# -1, 1,2-,2 for paired groups, with -1 paired with 1 etc
# or survival time, or categorical for unordered groups 1,2,3,..
# quantitative for continuous ordered y

# new resp.type- timearea
#
# params is a list containing user-defined parameters; at this point just has
#   min.foldchange (NULL means no fold change criterion, o/w should be> 1)
#
# s0 is the exchangeability factor; you can specify
# s0 as an actual value
# s0.perc, the percentile of sd values to use for s0
#  or if both s0 and s0.perc are null (the default), then s0 is automatically estimated

# permy is values of paired y, for within-style paired comparson (Naras- ignore)

# returns
#  evo= expected order statistics (length p=# of genes)
#  tt=numer/(sd+s0)  test statistics on original data (and the ingredients)
#  ttstar0= p by nperm matrix of test statistics on permted data
#  ttstar= ttstar0 with columns sorted (largest values in row 1)
#  also returns permuted values: foldchange.star,  ystar, sdstar, statusstar (for survival data)
        

  this.call=match.call()  

  resp.type=match.arg(resp.type)

  x=data$x
  y=data$y
  
  if(resp.type=="survival"){status=data$status}

  n <- nrow(x)
  ny <- length(y)

  sd <- NULL;numer <- NULL

  ystar.keep <- matrix(0,ncol=nperm,nrow=length(y))
  sdstar.keep <- matrix(0,ncol=nperm,nrow=nrow(x))

  status.keep <- NULL
  if(resp.type=="survival"){status.keep <- ystar.keep}

# estimate s0 if necessary
  
  if(is.null(s0)){
     if(resp.type=="twoclass"){junk <- ttest.func(x,y);numer <- junk$numer;sd <- junk$sd}
     if(resp.type=="onesample"){junk <- onesample.ttest.func(x,y);numer <- junk$numer;sd <- junk$sd}
     if(resp.type=="paired"){junk <- paired.ttest.func(x,y,diffsum.reg=diffsum.reg);
      numer <- junk$numer$sd <- junk$sd}
     if(resp.type=="survival"){junk <- cox.func(x,y,status);numer <- junk$numer;sd <- junk$sd}
     if(resp.type=="categorical"){junk <- fisher.func(x,y);numer <- junk$numer;sd <- junk$sd}
     if(resp.type=="quantitative"){junk <- quantitative.func(x,y);numer <- junk$numer;sd <- junk$sd}
     if(resp.type=="timearea"){junk <- timearea.func(x,y);numer <- junk$numer;sd <- junk$sd}
     
     if(!is.null(s0.perc)){s0 <- quantile(junk$sd,s0.perc)}
     if(is.null(s0.perc)){s0=est.s0(junk$tt,junk$sd)$s0.hat}
  }
 
# compute test statistics on original data

  if(resp.type=="twoclass"){tt <- ttest.func(x,y,s0=s0)$tt}
  if(resp.type=="onesample"){tt <- onesample.ttest.func(x,y,s0=s0)$tt}
  if(resp.type=="paired"){tt <- paired.ttest.func(x,y,s0=s0)$tt}
  if(resp.type=="survival"){tt <- cox.func(x,y,status,s0=s0)$tt}
  if(resp.type=="categorical"){tt <- fisher.func(x,y,s0=s0)$tt}
  if(resp.type=="quantitative"){tt <- quantitative.func(x,y,s0=s0)$tt}
  if(resp.type=="timearea"){tt <- timearea.func(x,y,s0=s0)$tt}

 
# compute test statistics on permuted  data
  
  ttstar <- matrix(0,nrow=nrow(x),ncol=nperm)
  foldchange.star=NULL
  if(params$min.foldchange>0){foldchange.star <- matrix(0,nrow=nrow(x),ncol=nperm)}
  for(b in 1:nperm){
    cat(c("perm=",b),fill=T)
    xstar <- x


    if(resp.type=="onesample"){
         ystar <- sample(c(-1,1),replace=T,size=ny)
         ttstar[,b] <- onesample.ttest.func(xstar,ystar,s0=s0)$tt
    }



    if(resp.type=="paired"){
      if(is.null(permy)){permy <- y}
      istar <- sample(c(-1,1),size=ny/2,replace=T)
      ystar <- rep(0,length(permy))
      ystar[permy<0] <- permy[permy<0]*istar
      for(j in 1:(ny/2)){ystar[permy==j] <-  -ystar[permy==-j]}
      ttstar[,b] <- paired.ttest.func(xstar,ystar,s0=s0)$tt
      ystar.keep[,b] <- ystar
      if(params$min.foldchange>0){foldchange.star[,b]=foldchange.paired(xstar,ystar,data$logged2)}

    }

    if(resp.type=="twoclass"){
      ystar <- sample(y,size=ny)
      junk <- ttest.func(xstar,ystar,s0=s0)
      ttstar[,b] <- junk$tt
      ystar.keep[,b] <- ystar
      sdstar.keep[,b] <- junk$sd
      if(params$min.foldchange>0){foldchange.star[,b]=foldchange.twoclass(xstar,ystar,data$logged2)}

    }

    if(resp.type=="survival"){
      o <- sample(1:length(y),size=ny)
      ttstar[,b] <- cox.func(xstar,y[o],status=status[o],s0=s0)$tt
      ystar.keep[,b] <- y[o]
      status.keep[,b] <- status[o]
    }
    if(resp.type=="categorical"){
      o <- sample(1:length(y),size=length(y))
      junk <- fisher.func(xstar,y[o],s0=s0)
      ttstar[,b] <- junk$tt
      ystar.keep[,b] <- y[o]
      sdstar.keep[,b] <- junk$sd
    }
    if(resp.type=="quantitative"){
      o <- sample(1:length(y),size=length(y))
      junk <- quantitative.func(xstar,y[o],s0=s0)
      ttstar[,b] <- junk$tt
      ystar.keep[,b] <- y[o]
      sdstar.keep[,b] <- junk$sd
    }
    if(resp.type=="timearea"){
      o <- sample(1:length(y),size=length(y))
      junk <- timearea.func(xstar,y[o],s0=s0)
      ttstar[,b] <- junk$tt
      ystar.keep[,b] <- y[o]
      sdstar.keep[,b] <- junk$sd
    } 

  }
  
# sort columns of statistics from permuted samples, and compute expected order statistics
  
  ttstar0 <- ttstar
  for(j in 1:ncol(ttstar)) {
    ttstar[, j] <- -1 * sort(-1 * ttstar[, j])
  }
  for(i in 1:nrow(ttstar)) {
    ttstar[i,  ] <- sort(ttstar[i,  ])
  }
  evo <- apply(ttstar, 1, mean)
  evo <- evo[length(evo):1]
  ystar <- ystar.keep;
  statusstar <- status.keep
  sdstar <- sdstar.keep

# estimation of pi0= prop of null genes
  
  qq<-quantile(ttstar,c(.25,.75))
  pi0 <- sum(tt>qq[1] & tt< qq[2])/(.5*length(tt))

 # compute fold changes, when applicable
  
  foldchange=NULL
  if(params$min.foldchange>0){
    if(resp.type=="twoclass"){ foldchange=foldchange.twoclass(data$x,data$y,data$logged2)}
    if(resp.type=="paired"){ foldchange=foldchange.paired(data$x,data$y,data$logged2)}
  }



  return(list(evo=evo,tt=tt,ttstar=ttstar,ttstar0=ttstar0,n=n,pi0=pi0,foldchange=foldchange, 
              foldchange.star=foldchange.star,  s0=s0, numer=numer,sd=sd, ystar=ystar,
              statusstar=statusstar, sdstar=sdstar, resp.type=resp.type, call=this.call))
}