#-------------------------------------------------------------------------------------------------------
#                                                                                                       
#               R-PROGRAM - APPLICATION                                                                 
#                                                                                                       
#-------------------------------------------------------------------------------------------------------
#
#
#  This is an R-programm for a two-stage test for experiments with a large number of hypotheses.        
#  First and second stage p-values have to be imported, then the programm results the test decision     
#  if a hypothesis is right or wrong for a given false discovery rate alpha.                            
#
#    functions:                                                                                         
#    ----------                                                                                         
#    pi0.est........estimated number of null hypotheses                                                 
#    g..............function to compute gamma as a function of gamma2                                   
#    int............integrates over the function g     
#    fdr............computes the estimated false discovery rate                                         
#    fdrminus.......computes the estimated false discovery rate minus the parameter alpha               
#    gamma2.........the intervall [0,1] is searched for a value which sets the function fdrminus zero.  
#                   the solution is gamma2                                                              
#    testdecision...vector resulting the test decision for each hypothesis.                             
#                   TRUE: The alternative hypothesis is true;                                           
#
#    Parameters:                                                                                        
#                                                                                                       
#    p1.............input vector of first stage p-values                                                
#    p..............input vector of the p-values from the pooled sample                                 
#                   (NA for hypotheses not carried to the 2nd stage)                                    
#    gamma1.........first stage selection boundary                                                      
#    n1.............first stage sample size                                                             
#    n2.............second stage sample size                                                            
#    lambda.........parameter to estimate the number of null hypotheses                                 
#    alpha..........targeted false discovery rate                                                       
#_______________________________________________________________________________________________________



pi0.est<-function(p1,lambda) min(sum(p1>lambda)/((1-lambda)*length(p1)),1)

g<-function(z,n1,n2,g2)         
    (1-pnorm((qnorm(1-g2)-sqrt(n1/(n1+n2))*z)/sqrt(n2/(n1+n2))))*dnorm(z)

int<-function(gamma1,n1,n2,g2)
    integrate(g,lower=qnorm(1-gamma1),upper=Inf,n1=n1,n2=n2,g2=g2) 

fdr<-function(p,p1,gamma1,n1,n2,g2,lambda)
    pi0.est(p1,lambda)*length(p1)*int(gamma1,n1,n2,g2)$value/(max(sum(p<g2,na.rm=TRUE),1))
    
fdrminus<-function(p,p1,gamma1,n1,n2,g2,lambda,alpha) 
    fdr(p,p1,gamma1,n1,n2,g2,lambda)-alpha

gamma2<-function(p,p1,gamma1,n1,n2,lambda,alpha) 
    uniroot(fdrminus,c(0.00000001,0.01),tol = 1e-200,p=p,p1=p1,gamma1=gamma1,n1=n1,n2=n2,lambda=lambda,alpha=alpha)

testdecision<-function(p,p1,gamma1,n1,n2,lambda,alpha)  
    (p1<gamma1) & (p<gamma2(p,p1,gamma1,n1,n2,lambda,alpha)$root)




#-------------------------------------------------------------------------------------------------------
#                   Example                                                                             
#-------------------------------------------------------------------------------------------------------


# The mean value of normal distributions is simulated with 4950 null hypotheses (mean=0)                
# and 50 alternative hypotheses (mean=1)                                                                
z<-c(rep(0,4950),rep(1,50))
mean1<-rnorm(5000,0,1/sqrt(4))+z
# P-value of the first stage
p1<-1-pnorm(mean1*sqrt(4),0,1)

# Simulation of second stage data
mean2<-rnorm(5000,0,1/sqrt(12))+z
# P-value from the pooled sample for hypotheses carried to the second stage
# hypotheses not carried to the second stage are NA
p<-ifelse(p1<0.1,1-pnorm((mean1*4/sqrt(16)+mean2*12/sqrt(16)),0,1),NA)

testdecision(p,p1,gamma1=0.1,n1=4,n2=12,lambda=0.5,alpha=0.05)