R/multtest.TcGSA.R

In borishejblum/TcGSA: Time-Course Gene Set Analysis

#'Computing the P-value of the Likelihood Ratios Applying a Multiple Testing
#'Correction
#'
#'This function computes the p-value of the likelihood ratios and apply a
#'multiple testing correction.
#'
#'
#'@param tcgsa 
#'a TcGSA object.
#'
#'@param threshold 
#'
#'the threshold at which the FDR or the FWER should be
#'controlled.
#'
#'@param myproc 
#'a vector of character strings containing the names of the
#'multiple testing procedures for which adjusted p-values are to be computed.
#'This vector should include any of the following: "\code{Bonferroni}",
#'"\code{Holm}", "\code{Hochberg}", "\code{SidakSS}", "\code{SidakSD}",
#'"\code{BH}", "\code{BY}", "\code{ABH}", "\code{TSBH}" or "\code{none}".  
#'"\code{none}" indicates no adjustment for multiple testing.  See
#'\code{\link[multtest:mt.rawp2adjp]{mt.rawp2adjp}} for details.  Default is
#'"\code{BY}", the Benjamini & Yekutieli (2001) step-up FDR-controlling
#'procedure (general dependency structures).  In order to control the FWER (in
#'case of an analysis that is more a hypothesis confirmation than an
#'exploration of the expression data), we recommend to use "\code{Holm}", the
#'Holm (1979) step-down adjusted p-values for strong control of the FWER.
#'
#'@param exact 
#'logical flag indicating whether the raw p-values should be computed from the 
#'exact asymptotic mixture of chi-square, or simulated (longer and not better).
#'Default is \code{TRUE} and should be preferred.
#'
#'@param nbsimu_pval 
#'the number of observations under the null distribution to
#'be generated in order to compute the p-values. Default is \code{1e+06}.
#'
#'@return \code{multtest.TcGSA} returns an dataframe with 5 variables.  The
#'rows correspond to the gene sets under scrutiny.  The 1st column is the
#'likelihood ratios \code{LR}, the 2nd column is the convergence status of the
#'model under the null hypothesis \code{CVG_H0}, the 3rd column is the
#'convergence status of the model under the alternative hypothesis
#'\code{CVG_H1}, the 4th column is the raw p-value of the mixed likelihood
#'ratio test \code{raw_pval}, the 5th column is the adjusted p-value of the
#'mixed likelihood ratio test \code{adj_pval}.
#'
#'@author Boris P. Hejblum
#'
#'@seealso \code{\link{TcGSA.LR}}, 
#'\code{\link[multtest:mt.rawp2adjp]{mt.rawp2adjp}},
#'\code{\link{signifLRT.TcGSA}}
#'
#'@importFrom multtest mt.rawp2adjp
#'
#'@importFrom stats rchisq
#'
#'@export multtest.TcGSA
#'
#'@examples
#'
#'if(interactive()){
#'data(data_simu_TcGSA)
#'
#'tcgsa_sim_1grp <- TcGSA.LR(expr=expr_1grp, gmt=gmt_sim, design=design, 
#'                           subject_name="Patient_ID", time_name="TimePoint",
#'                           time_func="linear", crossedRandom=FALSE)
#'                           
#'mtt <- multtest.TcGSA(tcgsa_sim_1grp, threshold = 0.05, 
#'                      myproc = "BY", nbsimu_pval = 1000)
#'mtt
#'}
#'
multtest.TcGSA <- function(tcgsa, threshold=0.05, myproc="BY", 
						   exact=TRUE, nbsimu_pval = 1e6){
	
	emp <- tcgsa[["fit"]]
	emp$raw_pval <- NULL
	
	func <- tcgsa[["time_func"]]
	group.var <- tcgsa[["group.var"]]
	separateSubjects <- tcgsa[["separateSubjects"]]
	time_DF <- tcgsa[["time_DF"]]
	
	if(!exact){
		if(is.null(group.var)){
			if(!separateSubjects){
				if(func=="linear"){
					theodist <- c(stats::rchisq(nbsimu_pval/2,df=1), stats::rchisq(nbsimu_pval/2,df=2))
				}else if(func=="cubic"){
					theodist <- rchisqmix(nbsimu_pval,3,3)
				}else{
					theodist <-rchisqmix(nbsimu_pval,time_DF,time_DF)
				}
			}else{
				if(func=="linear"){
					theodist <- c(stats::rchisq(nbsimu_pval/2,df=0), stats::rchisq(nbsimu_pval/2,df=1) )
				}else if(func=="cubic"){
					theodist <- rchisqmix(nbsimu_pval,0,3)
				}else{
					theodist <-rchisqmix(nbsimu_pval,0,time_DF)
				}
			}
		}else{
			nbgp <- length(levels(group.var))
			if(func=="linear"){
				theodist <- rchisqmix(nbsimu_pval, 1*(nbgp-1), 0)
			}else if(func=="cubic"){
				theodist <- rchisqmix(nbsimu_pval, 3*(nbgp-1), 0)
			}else{
				theodist <-rchisqmix(nbsimu_pval, time_DF*(nbgp-1), 0)
			}
		}
		
		emp$raw_pval <- unlist(lapply(emp$LR, FUN=pval_simu, theo_dist=theodist))
		
		
	}else{
		if(is.null(group.var)){
			if(!separateSubjects){
				if(func=="linear"){
					emp$raw_pval <- sapply(emp$LR, pchisqmix, s=1, q=1, lower.tail=FALSE)
				}else if(func=="cubic"){
					emp$raw_pval <- sapply(emp$LR, pchisqmix, s=3, q=3, lower.tail=FALSE)
				}else{
					emp$raw_pval <- sapply(emp$LR, pchisqmix, s=time_DF, q=time_DF, lower.tail=FALSE)
				}
			}else{
				if(func=="linear"){
					emp$raw_pval <- sapply(emp$LR, pchisqmix, s=0, q=1, lower.tail=FALSE)
				}else if(func=="cubic"){
					emp$raw_pval <- sapply(emp$LR, pchisqmix, s=0, q=3, lower.tail=FALSE)
				}else{
					emp$raw_pval <- sapply(emp$LR, pchisqmix, s=0, q=time_DF, lower.tail=FALSE)
				}
			}
		}else{
			nbgp <- length(levels(group.var))
			if(func=="linear"){
				emp$raw_pval <- sapply(emp$LR, pchisqmix, s=1*(nbgp-1), q=0, lower.tail=FALSE)
			}else if(func=="cubic"){
				emp$raw_pval <- sapply(emp$LR, pchisqmix, s=3*(nbgp-1), q=0, lower.tail=FALSE)
			}else{
				emp$raw_pval <- sapply(emp$LR, pchisqmix, s=time_DF*(nbgp-1), q=0, lower.tail=FALSE)
			}
		}
	}
	
	if(myproc=="none" | length(emp$raw_pval)==1){
		emp$adj_pval <- emp$raw_pval
	}
	else{
		adj_pval <- mt.rawp2adjp(emp$raw_pval, proc=c(myproc), alpha=threshold)
		emp$adj_pval <- adj_pval$adjp[order(adj_pval$index),2]
	}
	return(emp)
}