#R code to calculate the correction factor as described in the paper
#Kim and Schliekleman (In Press) Prioritizing Hypothesis Tests in High Throughput Data. Bioinformatics

#This file contains two functions. FWERfunc() returns the quantity [1-FWER(lambda)]-[1-alpha0] 
# as a function of the correction factor lambda. This is used in the second function CorrectionFactor to solve for the 
# value of the correction factor. 




FWERfunc <- function(lambda,s,alpha0) 
{#Returns the mathematical function whose root is the correction factor for optimized filtering. 

#input:
#   lambda=corretion factor
#		s=the vector of test set sizes
#		alpha0=nominal FWER

#output:   
#   [1-FWER(lambda)]-[1-alpha0] 
#   FWER(lambda) is the FWER when the correction factor is lambda. The goal is to have [1-FWER(lambda)]=[1-alpha0]. 
#   The value of lambda that makes [1-FWER(lambda)]=[1-alpha0] is the correct value. 

  k=length(s)
  t=(1-(alpha0/(lambda*s[1])))^s[1]
  for(j in 1:(k-1))
   {
     t=t*(1-(alpha0/(lambda*s[j+1])))^(s[j+1]-s[j])
   }
  #t=1-FWER. This is equation (2) in the paper. We solve FWER=alpha or 1-FWER=1-alpha

  return(t-(1-alpha0))
}



CorrectionFactor<-function(set,alpha,tol0)
{#function to calculate the correction factor for optimized filtering. 

#input: 
#   set=vector s=[s1,s2,s3...] specifying the tested sets. That is, the first tested set is 1 to s1, the second tested set is 1 to s2, etc. 
#	  alpha=the desired FWER for set after the weight is applied. 
#   tol0: the tolerance factor for the uniroot function. 
  
#output:  
#   opt.lambda = the correction factor  

#This function requires an additional function FWERfunc that returns [1-FWER(lambda)]-[1-alpha0]. This is used by the function
# uniroot() in order to solve for the correction factor lambda. 
  

if(length(set)==1) return(1)
opt.lambda=uniroot(FWERfunc,lower=0.01,upper=length(set),s=set,alpha0=alpha,tol=tol0)$root 
return(opt.lambda)

}



#Usage:
#Calculate the correction factor for the case with target FWER=0.05, test set=(20 t0 1000 in incrments of 20), tolerance for uniroot=0.001. 
CorrectionFactor(seq(from=20,to=1000,by=20),0.05,0.001)