# TODO: Add comment
# 
# Author: E.Korsching 10.9.2009, adopted from R.Voss
###############################################################################

adjust.p <- function(rawp, proc = c("Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY"), sort = F)
{
	# Adjusted p-values for simple multiple testing procedures
	# Output: One vector with adjusted p-values
	# use pvalues.adjust to to get a matrix with a set of all (in proc) corrections for one raw-p-vector
	# from pvalues.adjust.ssc and mt.raw2p.adjp in multtest.ssc from ArrayAnalyzer
	# see also plot.pvalues and pvalues.reject
	# sort=T returns the adjusted p-values in increasing order
	# speed up the loops use cummax and cumin (from adjust.p.ssc)

	# AUTHOR's
	# 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.

	m <- length(rawp)
	index <- order(rawp)
	spval <- rawp[index]

	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)
		#  tmp[1] <- min(na.rm = T, c(m * spval[1], 1))
		#  for(i in 2:m)
		#   tmp[i] <- max(na.rm = T, c(tmp[i - 1], min(na.rm = T, c((m - i + 1) * spval[i], 1))))
	}
	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))
		#  tmp <- spval
		#  for(i in (m - 1):1)
		#   tmp[i] <- min(na.rm = T, c(tmp[i + 1], (m - i + 1) * spval[i], 1))
	}
	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)
		#  tmp <- spval
		#  tmp[1] <- 1 - (1 - spval[1])^m
		#  for(i in 2:m)
		#   tmp[i] <- max(na.rm = T, c(tmp[i - 1], 1 - (1 - spval[i])^(m - i + 1)))
	}
	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))
		#  tmp <- spval
		#  for(i in (m - 1):1)
		#   tmp[i] <- min(na.rm = T, c(tmp[i + 1], (m/i) * spval[i], 1))
	}
	if(is.element("BY", proc)) {
		# Benjamini & Yekutieli (2001)        
		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))
		#  tmp <- spval
		#  a <- sum(1/(1:m))
		#  tmp[m] <- min(na.rm = T, c(a * spval[m], 1))
		#  #noting we need to set tmp[m]
		#  for(i in (m - 1):1)
		#   tmp[i] <- min(na.rm = T, c(tmp[i + 1], ((m * a)/i) * spval[i], 1))
	}

    if(sort){return(tmp)}else{return(tmp[order(index)])}
}