# TODO: Add comment
# 
# Author: E.Korsching Mar 14, 2012
###############################################################################


mtt <- function(
		# Procedure to calculate differences between multiple groups of signal values
		# see the results
		x,						# signal data, data.frame
		categories,				# group structure for the given columns e.g. c("a","b","c","a",...)
		permutations=1,			# permutation test of results, 1: no testing
		reg=NULL,					# rules: regulation wanted between groups (categories), all other pairwise comparisons: no regulation
		method="lsd",			# mean analysis after one way ANOVA - "t-test", "bonferroni", "lsd", "sidak", "tukey", "scheffe"
		proc="BH",					# multiple testing correction - "noadjust", "Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY"
		var.equal=T,				####automate?
		palpha=0.05,				# alpha error threshold
		plot.it=T,					# plot graphs
		gin,						# genom information, annotation
		label.cols=c("Gene.Symbol","label"),			# annotation columns of gin (one or more)
		mf.row=c(2,2),				# page design for plotting
		cex=0.5,					# font size
		plottype="interval",		# LSD plot types ("fold.change","bar","interval")
		result="short"			# return: "full": full results-permutations=1, "short": short results-permutations=1
#		seed=100,					# seed  - (see permutation procedure)
#		diagnostic=T,				# diagnostics on (T) - (see permutation procedure)
#		maxt=F						# compare t-values (T) or compare p values (F) - (see permutation procedure)
)
{
	# ini
	color.set(color="cgh", example=F, colorcode=F)	#set color space
	
	# check parameters
	method <- match.arg(method,c("t-test","bonferroni","lsd","sidak","tukey","scheffe"))	#lsd info: protected (<alpha!) Fisher's LSD (Least Significant Difference) test
	proc <- match.arg(proc,c("noadjust","Bonferroni", "Holm", "Hochberg", "SidakSS", "SidakSD", "BH", "BY"))	#BH,BY: Benjamini-Hochberg/Yekutieli
	
	nr <- nrow(x)
	nc <- ncol(x)
	if(nc < 4){ cat("\n 1: x must exhibit at least 2 groups with at least 2 columns each"); break }
	if(missing(gin)){ cat("\n 1b: gin data must be provided - row names like in x"); break }
	
	rnames <- row.names(x)
	group.label <- unique(categories)[order(unique(categories))]			# save the category names as group labels    ??? uc used 2 times ?????
	k <- length(group.label)										# number of groups
	if(k < 2){ cat("\n 2: x must exhibit at least 2 groups"); break }
	comparisons <- k*(k-1)/2						# number of comparisons having k groups		
	comparison <- t(outer(group.label,group.label,paste,sep="."))[lower.tri(matrix(0,k,k))]		# compute all unique combinations of categories
	
	#order categories  and  order x according to groups(categories)
	categ.order.vec <- order(categories)
	categories <- categories[categ.order.vec]
	x <- x[,categ.order.vec]
	
	#calculate group sizes
	categ.ordered.n <- as.vector(table(categories))		# as.vector(table(c("b","b","a","b","a","b","b")))
	if(sum(categ.ordered.n < 2)>0){ cat("\n 3: each group must exihibit at least 2 columns"); break }
	
	categories <- as.factor(categories)				# category col labels to factor

#	xname <- deparse(substitute(x))				# extract data.frame name
#	cat("\nANOVA on: ",xname)		# cause a problem by complex data: runtime 1h instead of 18s
	
	cat("\nANOVA : group labels:", group.label )
	cat("\nordered categories:",as.character(categories) )
	
	
	####### ANOVA (one-way ANOVA) #######
	
	# assumptions: independence, normality (the distributions of the residuals are normal), variance equal, always one-tailed
	#  H0: mean1 = mean2 = ... = meanN (all),  HA: mean1 != mean2 != ... != meanN (at least one)
	# If H0 is true: all groups are from the same population -  the groups with their different sample means simply represent points out of the same sampling distribution
	# If H0 is true: "between group variance" will be equal to the "within group variance"
	# "between group variance" (or Mean Square due to Treatments or MSTR) is an estimate of the variance of the population if H0 is true
	# "within group variance" (or Mean Square Error or MSE)
	# If H0 is true: "between group variances" = (or close to) "within group variances"
	# test statistic: F = MSTR / MSE	[F statistic]
	# H0: means are equal: accept H0: F-calculated < F-critical or P-value > alpha
	
	anova <- data.frame(t( apply(x, 1, anova.ek, categories=categories, categ.ordered.n=categ.ordered.n, group.label=group.label) ))
	# return values: fp.value, f.value, between.var(msb), within.var(mse), between.df1, within.df2, between.ssb, within.sse, within.means(vec! check length)
	
	# export some data from anova structure
	within.df2 <- anova[1,"within.df2"]
	within.means <- data.frame(anova[ ,9:ncol(anova)])
	dimnames(within.means)[[2]] <- group.label		#col names must be restored
	
	
	####### Test of significant ANOVA results on mean differences #######
	
	# pairwise multiple testing : critical point
	if(method=="t-test"){ crit.point <- qt(1-palpha/2, within.df2) }							# t-test
	if(method=="lsd"){ crit.point <- qt(1-palpha/2, within.df2) } 								# LSD
	if(method=="bonferroni"){ crit.point <- qt(1 - palpha/(2*p), within.df2) }					# Bonferroni - less conservative / more powerful, also for t tests, number of paired tests are low
	if(method=="sidak"){ crit.point <- qmvt(1 - palpha, p, df.residual, rho = 0) }				# Sidak
	if(method=="tukey"){ crit.point <- (1/sqrt(2)) * qtukey(1 - palpha, k, within.df2) }		# Tukey - sampla size equal (unequal: Tukey-Kramer), less conservative, if many/all pairs will be tested
	if(method=="scheffe"){ crit.point <- sqrt((k-1)* qf(1 - palpha, k-1, within.df2)) }			# Scheffe - when sample sizes are unequal, when most conservative test is desired
	
	#R also include: Tukey HSD test as a post-hoc test: TukeyHSD(your.aov) 
	
	# info test
	cat("\n\nTest of significant ANOVA results on mean differences")
	cat("\ngroups:", k,", comparisons:", comparisons,", method:", method,", (within) degrees of freedom:", within.df2, ", critical point:", round(crit.point,2) )
	
	# fold change
	fold.change <- NULL
	name.fc <- NULL
	cat("\n\npairs 1: ")				# fc: no dot.product procedure available ####really?
	for(i in 1:k){
		for(j in i:k){
			if(group.label[i] != group.label[j]){
				name.fc <- c(name.fc,paste(group.label[i],group.label[j],sep="."))				# names vector
				fold <- within.means[,group.label[i]]/within.means[,group.label[j]]		# fold change between groups
				fold <- ifelse(fold < 1,-1/fold,fold)
				fold.change <- cbind(fold.change,fold)					# fold.change within groups
			}
		}
	}
	cat(name.fc)					#show work (for comparison purpose with "comparison <- t(outer..." variable)
	fold.change <- data.frame(fold.change)
	dimnames(fold.change) <- list(dimnames(within.means)[[1]],name.fc)
	cat("\npairs 2:",comparison )
	
	ni <- as.vector(table(match(categories,group.label)))		# as.vector(table(match(c("a","a","a","b","a","b","b"),c("a","b")))) => {4, 3}
	nmat <- matrix(rep(ni,dim(x)[1]),nrow=dim(x)[1],byrow=T)	# matrix as long as x with identical ni vectors
	dimnames(nmat) <- list(row.names(x), group.label)
	
	Lmat <- contrast.matrix(categories, comparison)
	
	pair.diff <- as.matrix(within.means) %*% Lmat			#working - but uncommon - extra transpose not nessecary!
	pair.diff <- data.frame(pair.diff)
	dimnames(pair.diff) <- list(row.names(x), comparison)
	
	if(method=="t-test"){
		stderr  <- NULL
		if(var.equal){	#standard t test
			df0 <- as.matrix(nmat) %*% abs(Lmat) -2			# n1+n2-2   -- variance equal
			stderr <- sqrt((as.matrix(within.vars*(nmat-1)) %*% abs(Lmat))/df0 * as.matrix(1/nmat) %*% abs(Lmat))
		} else {		#WELCH -  variance unequal: by loosing minimal power -use this method as standard method, gives also some ci intervals
			var.n <- within.vars/nmat							# vars/ni
			stderr <- sqrt(as.matrix(var.n) %*% abs(Lmat))	# Welch method: sqrt( var.x/nx + var.y/ny )
			# df0 <- (varx/nx+vary/ny)^2/((varx/nx)^2/(nx-1)+(vary/ny)^2/(ny-1))
			df0 <-  (as.matrix(var.n) %*% abs(Lmat))^2/as.matrix((var.n)^2/(nmat-1)) %*% abs(Lmat)
		}
		
		statistic <- pair.diff/stderr
		
		pval <- 2*(1- pt(q=as.matrix(abs(statistic)), df=df0) )
		pvalues <- matrix( pval, nrow=dim(pair.diff)[1], ncol=dim(pair.diff)[2], byrow=F, dimnames=list(row.names(x), comparison) )
		cp <- qt((1- palpha/2), df0)										# critical point for pairwise t-test (non-simultaneous)
		crit.point <- matrix( cp, nrow=dim(pair.diff)[1], ncol=dim(pair.diff)[2], byrow=F, dimnames=list(row.names(x), comparison) )
	}else{	
		stderr <- sqrt(anova[ ,"within.var"]*as.matrix(1/nmat) %*% abs(Lmat))		# mse*(1/nx+1/ny) standard error for all other methods exept t-test						
	}
	
	ci <- stderr*crit.point				# confidence limit (two sided)
	
	if(method=="lsd" | method=="bonferroni"| method=="sidak"){
		statistic <- pair.diff/stderr
		pval <- round( 2*(1- pt(q=as.matrix(abs(statistic)), df=within.df2) ) ,5)
		if(method=="bonferroni" | method=="sidak"){ pval <- 1 - (1 - pval)^p }
		pvalues <- matrix( pval, nrow=dim(pair.diff)[1], ncol=dim(pair.diff)[2], byrow=F, dimnames=list(row.names(x), comparison) )	# lsd method Millard/Krause p.342
	}
	
	# for tukey, scheffe no pvalues can be computed 
	#  estimate: artificial pvalue of 0.0001  if the ci interval does not include zero
	
	ci <- data.frame(ci)
	names(ci) <- comparison
	lower <- data.frame(pair.diff-ci)
	upper <- data.frame(pair.diff+ci)
	
	if(method=="tukey" | method=="scheffe"){
		pvalues <- lower	# initialize variable
		pvalues <- ifelse( (lower < 0 & upper < 0) | (lower > 0 & upper > 0), 0.0001, 1 )	# fill variable
		pvalues <- data.frame(pvalues)
	}
	
	statistic <- pair.diff/stderr
	
	
	####### multiple testing correction of p value #######
	if(proc!="noadjust"){
		cat("\nAdjusting p-values, method:",proc)
		pvalues <- apply(pvalues, 2, adjust.p, proc=proc)
	}
	
#	## if permutations >1 take pvalues from this test for further selections e.g. in rnr.R
#	if(permutations > 1){
#		#### pboot <- call perm sub function here
#		pvalues <- pboot/permutations
#		pvalues[pvalues == 0] <- 1 		# nothing found, no sig
#	}
	
	# add row.names etc.
	pvalues <- data.frame(round(pvalues,5), row.names=rnames, check.rows=F)
	
	
	####### filter the pvalues with reg/noreg criteria (direct p or corrected p or bootstrap p) #######
	psig <- rnr(
			x=pvalues,
			reg=reg,
			sorting=T,
			palpha=palpha,
			return.logic=T		#return filter vector with T/F entries per row
	)		
	nsig <- sum(psig)	
	cat("\nSignificant results: ",nsig," of:",nr)
	cat("\n by applying mean test method '",method,"', multiple testing correction '",proc,"'")
	cat("\n and filtering by reg/noreg criteria (based on ANOVA p and given alpha threshold)")
	
	# create significants
	if(nsig > 0){
		rn.sig <- rnames[psig]
		pvalue.sig <- round(pvalues[psig,],5)		# select sig pvalues concerning mean test, eventually multiple testing correction and rules 
		anova.sig <- anova[psig,]
		lower.sig <- lower[psig,]
		upper.sig <- upper[psig,]
		pvalue.boot.sig <- round(pvalues[psig,],5)
		within.means.sig <- within.means[psig,]
		pair.diff.sig <- pair.diff[psig,]
		fc.sig <- fold.change[psig,]
		se.sig <- stderr[psig]
	}else{
		cat("\nNo result - no pvalues are saved")
	}
	
	
	####### plot #######
	if(plot.it & nsig>0){
		par(mfrow=mf.row)
		nrs <- sum(psig)
		
		for(i in 1:nrs){
			figure <- i
			main <- paste(rn.sig[i], substring(gin[rn.sig[i],"Gene.Symbol"],1,15), sep="  ")

			if(nr==1){		#one/many line(s) data types
				stat <- data.frame(t(rbind(pair.diff=pair.diff.sig[i,],lower=lower.sig[i,],upper=upper.sig[i,],pvalue=pvalue.sig,fc=fc.sig[i,],se=se.sig[i])))
			}else{
#				stat <- data.frame(t(rbind(pair.diff.sig=pair.diff.sig[i,],lower=lower.sig[i,],upper=upper.sig[i,],pvalue=pvalue.sig[i,],fc=fc.sig[i,],se=se.sig[i],p.perm=pvalue.boot.sig[i,])))
				stat <- data.frame(t(rbind(pair.diff=pair.diff.sig[i,],lower=lower.sig[i,],upper=upper.sig[i,],pvalue=pvalue.sig[i,],fc=fc.sig[i,],se=se.sig[i])))
			}
			
			stat$SIG <- factor(cut(stat$pvalue, breaks = c(0, 0.0001, 0.001, 0.01, 0.05, 100), labels = c("****", "***", "**", "*", ""), include.lowest = T))
			
			if(method=="tukey" | method=="scheffe"){ 
				stat$SIG <- ifelse((stat$lower < 0 & stat$upper < 0) |(stat$lower > 0 & stat$upper > 0),"****","")
				stat$pvalue <- rep("",p)
			}
			
			if(method=="t-test"){ method.str <- paste("(two-sided t-test,var.eq.:",var.equal,")") }
			if(method=="lsd"){ method.str <- "(Fisher's Protected LSD)" }
			if(method=="bonferroni"){ method.str <- "(Bonferroni method)" }
			if(method=="tukey"){ method.str <- "(Tukey method)" }
			if(method=="sidak"){ method.str <- "(Sidak method)" }
			if(method=="scheffe"){ method.str <- "(Scheffe method)" }
			
			if(method=="lsd"){
				sub <- paste("Fig.",figure,". 95 % non-simultaneous confidence intervals and",sep="")
			}else{
				sub <- paste("Fig.",figure,". 95 % simultaneous confidence intervals and",sep="")
			}
			
			if(plottype=="bar"){ sub <- paste("Fig.",figure,". Difference and pooled standard errors of group means  and",sep="") }
			if(plottype=="fold.change"){ sub <- paste("Fig.",figure,". Fold change values of group means and t-test",sep="") }
			if(plottype=="fold.change"){
				sub <- paste(sub,"\nsig. indicators for ",comparisons," specified comparisons ",method.str,sep="")
			}else{
				sub <- paste(sub,"\np values for ",comparisons," specified comparisons ",method.str,sep="")
			}
			
#			if(permutations > 1){
#				T1.err <- round(sum(stat$p.perm < palpha)/length(stat$p.perm),4)
#				sub <- paste(sub,"\nRight p-values from permutation test: ",permutations,"samples, incl.",ngenes," genes,TI-error=",T1.err)		
#				if(maxt)sub <- paste(sub,", maxT ") else sub <- paste(sub,", minP ")
#				if(proc != "noadjust")sub <- paste(sub,", method:",proc)
#				cor.p.perm <- round(cor(stat$pvalue,stat$p.perm),2)
#				sub <- paste(sub,"\nCorrelation between ",method.str," and pvalues from permutation test= ",cor.p.perm)		
#			}
			
			sub <- paste(sub,"\np alpha: ",palpha,", MSE: ",round(anova.sig$within.var[i],5),", crit.point: ",round(crit.point,3),", ANOVA p < ",round(anova.sig$fp.value[i],5),sep="")
			if(proc!="noadjust"){ sub <- paste(sub,",corrected by:",proc) }
			
#			stat[is.infinite(stat)] <- 0		#does not work on list/data.frame
			
			if(plot.it){
				plot.lsd(
						x=stat,
						main=main,
						sub=sub,
						figure=i,
						sig=T,
						palpha=palpha,
						plottype=plottype,
						cex=cex
				)
			}
		}
	}
	
	#write
	if(method=="t-test"){ method.str <- paste("(two-sided t-test,var.eq.:",var.equal,")") }
	if(method=="lsd"){ method.str <- "(Fisher's Protected LSD)" }
	if(method=="bonferroni"){ method.str <- "(Bonferroni method)" }
	if(method=="tukey"){ method.str <- "(Tukey method)" }
	if(method=="sidak"){ method.str <- "(Sidak method)" }
	if(method=="scheffe"){ method.str <- "(Scheffe method)" }
	
	cat("\n\np-values for ",comparisons," comparisons ",method.str," alpha: ",palpha)
	cat("\nType I error, probability of significant regulation if H0 is true")
	cat("\nprobability of comparisons indicated as significant  with p < 0.05\n")
	
#	if(permutations >1){ FDR <- median(pboot)/m }			# false discovery rate, see Sorin Draghici p.228
	
	T1.error <- apply( pvalues, 2, function(x){ round(sum(x < 0.05)/length(x),4) } )
	print(T1.error)
	
#	if(permutations > 1){
#		cat("\nOverall false discovery rate (FDR)")
#		cat("\nfor",m,"genes and ",permutations," permutations =",round(FDR,6))
#	}
	
	#end - clean up
	palette("default")
	cat("\n >>>> end >>>> \n")
	if(result=="short"){
		if(nsig > 0){
			return(list(
						p.Significant=pvalue.sig,
						mean.sig=within.means.sig,
						difference.sig=pair.diff.sig,
						fold.change.sig=fc.sig,
						GIN=gin[psig, label.cols, drop=F]
					) )
		}else{
			return(list(
						p.all=pvalues,
						mean.all=within.means,
						difference.all=pair.diff,
						fold.change.all=fold.change
					) )
		}
	}
	if(result=="full"){
		if(nsig > 0){
			return(list(
						p.all=pvalues,
						p.Significant=pvalue.sig,
						mean.sig=within.means.sig,
						difference.sig=pair.diff.sig,
						fold.change.sig=fc.sig,
						GIN=gin[psig, label.cols, drop=F]
					) )
		}else{
			return(list(
						p.all=pvalues,
						mean.all=within.means,
						difference.all=pair.diff,
						fold.change.all=fold.change
					) )
		}
	}
}