# TODO: Add comment
# 
# Author: Korsching 10.9.2009, 22.07.2014, based on R.Voss
###############################################################################



adjust.p <- function(rawp, proc = c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY"))
{
	# adjusted p-values for multiple testing procedures
	# input: vector with raw p values , output: vector with adjusted p-values
	# from pvalues.adjust.ssc and mt.raw2p.adjp in multtest.ssc from ArrayAnalyzer

	# Sandrine Dudoit, http://www.stat.berkeley.edu/~sandrine, 
	# Yongchao Ge, gyc@stat.berkeley.edu.
	#  http://www.maths.lth.se/help/R/.R/library/multtest/html/mt.rawp2adjp.html 
	#  R.Voss  rev. 23.10.03  10:20
	#  DATE: 11-03-04   R.Voss
	#  from mt.raw2p.adjp in multtest.ssc from ArrayAnalyzer
	# see also plot.pvalues and pvalues.reject and adjust.pvalues
	# sort=T returns the adjusted p-values in increasing order
	# Adjusted p-values are computed for simple FWER and FDR controlling procedures 
	# based on a vector of raw (unadjusted) p-values.
	# DESCRIPTION
	# This function computes adjusted p-values for simple multiple testing procedures 
	# from a vector of raw (unadjusted) p-values. The procedures include the Bonferroni, 
	# Holm (1979), Hochberg (1988), 
	# and Sidak procedures for strong control of the family-wise Type I error rate (FWER), 
	# and the Benjamini & Hochberg (1995) and 
	# Benjamini & Yekutieli (2001) procedures for (strong) control of the false discovery rate (FDR). 
	# USAGE
	# mt.rawp2adjp(rawp, proc=c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY"))
	# ARGUMENTS
	# rawp A vector of raw (unadjusted) p-values for each hypothesis under consideration.
	# These could be nominal p-values, for example, from t-tables, or permutation p-values 
	# as given in mt.maxT and mt.minP. If the mt.maxT or mt.minP functions are used, raw p-values
	# should be given in the original data order, rawp[order(index)]. 
	# proc 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:
	# "Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY".  
	# VALUE:
	# adjp: A matrix of adjusted p-values, with rows corresponding to hypotheses
	#  and columns to multiple testing procedures. 
	#  sort == T: Hypotheses are sorted in increasing order of their raw (unadjusted) p-values  
	# DETAILS
	# Family Wise ERrors  FWER
	# Bonferroni: Bonferroni single-step adjusted p-values for strong control of the FWER
	# Holm:  Holm (1979) step-down adjusted p-values for strong control of the FWER. 
	# Hochberg: Hochberg (1988) step-up adjusted p-values for strong control of the FWER
	#     (for raw (unadjusted) p-values satisfying the Simes inequality). 
	# SidakSS: Sidak single-step adjusted p-values for strong control of the FWER
	#     (for positive orthant dependent test statistics). 
	# SidakSD: Sidak step-down adjusted p-values for strong control of the FWER
	#     (for positive orthant dependent test statistics). 
	# False discovery rates FDR
	# BH:   adjusted p-values for the Benjamini & Hochberg (1995) step-up FDR 
	#     controlling procedure (independent and positive regression dependent test statistics). 
	# BY:   adjusted p-values for the Benjamini & Yekutieli (2001) step-up FDR 
	#     controlling procedure (general dependency structures). 
	# REFERENCES
	# Y. Benjamini and Y. Hochberg (1995). Controlling the false discovery rate: a practical 
	# and powerful approach to multiple testing. J. R. Statist. Soc. B. Vol. 57: 289-300.
	# Y. Benjamini and D. Yekutieli (2001). The control of the false discovery rate in multiple
	# hypothesis testing under dependency. Annals of Statistics. Accepted.
	# S. Dudoit, J. P. Shaffer, and J. C. Boldrick (Submitted). Multiple hypothesis testing in 
	# microarray experiments.
	# Y. Ge, S. Dudoit, and T. P. Speed. Resampling-based multiple testing for microarray data 
	# hypothesis, Technical Report #633 of UCB Stat. http://www.stat.berkeley.edu/~gyc
	# Y. Hochberg (1988). A sharper Bonferroni procedure for multiple tests of significance, 
	# Biometrika. Vol. 75: 800-802.
	# S. Holm (1979). A simple sequentially rejective multiple test procedure. Scand.
	# J. Statist.. Vol. 6: 65-70.
	
	# exclude NA and NaN values
	rawp.na <- is.na(rawp)
	rawp.na.len <- length(rawp.na)
	if(sum(is.na(rawp))>0){ rawp <- rawp[!rawp.na] }
	
	# ini
	m <- length(rawp)
	index <- order(rawp)
	spval <- rawp[index]
	
	# methods
	if(is.element("Bonferroni", proc)) {
		cat("\nAdjusting p-values: Bonferroni single-step adjusted p-values for strong control of the FWER \n")
		tmp <- m * spval
		tmp[tmp > 1] <- 1
	}
	if(is.element("Holm", proc)) {
		cat("\nAdjusting p-values: Holm (1979)step-down adjusted p-values for strong control of the FWER \n")
		tmp <- (m:1) * spval
		tmp[tmp > 1] <- 1
		tmp <- cummax(tmp)
	}
	if(is.element("Hochberg", proc)) {
		cat("\nAdjusting p-values: Hochberg (1988)step-up adjusted p-values for strong control of the FWER \n")
		tmp <- rev((m:1) * spval)
		tmp[tmp > 1] <- 1
		tmp <- rev(cummin(tmp))
	}
	if(is.element("SidakSS", proc)) {
		cat("\nAdjusting p-values: Sidak single-step adjusted p-values for strong control of the FWER \n")
		tmp <- 1 - (1 - spval)^m
	}
	if(is.element("SidakSD", proc)) {
		cat("\nAdjusting p-values: Sidak step-down adjusted p-values for strong control of the FWER \n")
		tmp <- 1 - (1 - spval)^(m:1)
		tmp <- cummax(tmp)
	}
	if(is.element("BH", proc)) {
		cat("\nAdjusting p-values: Benjamini & Hochberg (1995) step-up FDR \n")
		tmp <- rev((m/(1:m)) * spval)
		tmp[tmp > 1] <- 1
		tmp <- rev(cummin(tmp))		# cummin() : NA sensitive !
	}
	if(is.element("BY", proc)) {
		cat("\nAdjusting p-values: Benjamini & Yekutieli (2001) step-up FDR \n")
		a <- sum(1/(1:m))
		tmp <- rev((m * a)/(1:m) * spval)
		tmp[tmp > 1] <- 1
		tmp <- rev(cummin(tmp))
	}
	
	# restore initial order
	tmp <- tmp[order(index)]
	
	# reinclude NA values   [suboptimal]
	tmp2 <- vector(mode="numeric",length=rawp.na.len)
	tmp2[rawp.na] <- NA
	counter <- 1
	for(i in 1:rawp.na.len){
		if(!is.na(tmp2[i])){
			tmp2[i] <- tmp[counter]
			counter <- counter +1
		}
	}
	
    return(tmp2)
}