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