#FILENAME: 20100810_TEQ_for_Ebayes.R
#Andrew Winterman


# This is a Two One Sided T-test (TOST) of equivalence. 
#Meaning it tests for evidence of equivalence rather than 
#absence of evidence of difference, and does so using a 
#pair of t - tests, one testing null:
#	H_0 :      D <= - eps
#versus alternative
#	H_A : 		D > -eps

#The other testing:
#	H_0 :      D >= eps
#versus alternative
#	H_A : 		D < eps

#D is typically a difference or contrast statistics, 
#often the difference between the means of two groups

#We reject the null that abs(D) > eps if and only if both 
#tests reject their nulls. 
#Hence we reject if the confidence interval for D is contained in (-eps : eps).
#Reference: Stefan Wellek's book:
#Testing Statistical Hypothoses of Equivalence and Noninferiority
#Second Edition
#year 2010
#publisher  CRC Press


tEQ <- function(Fit.bayes, eps )
{
#	Fit.bayes is the output of EBayes, eps defines the size 
#	of the equivalence region. If the confidence interval for D is contained 
# within (-eps,eps) we call the two groups equivalent. 

#	 eps must either be a single number 
#	or a vector equal to the length equal to the number of rows 
# in Fit.bayes$coefficients.
# If you specify more than one contrast, this funtion will calculate 
#	p-values for each of them.
	eps <- abs(eps)
	d <- as.matrix(Fit.bayes$coeff) 
#		d holds the contrast coefficients 
	SE <- as.matrix(sqrt(Fit.bayes$s2.post)*Fit.bayes$stdev.unscaled)
# 		SE holds the ebayes modified standard errors.  
#		s2.post is the residual standard error from lmFit modified by ebayes
#		stdev.unscaled is the unscaled standard deviations 
#		computed from the unscaled covariance matrix.	 
	dgf <- Fit.bayes$df.resid
#		Is eps specified per feature or singly?
	if( length(eps) != 1 & length(eps) != nrow(d) )
	{
		stop("eps must either have length one or nrow(d)")
	}			
	p <- matrix(NA, nrow = nrow(d), ncol = ncol(d))
#		p will hold p-values for the test of equivalence 
	for ( i in 1:ncol(d) ) 
	{		  #Moving through each contrast (colums of d) 
		SEc <- SE[,i]
		dc <- d[,i]
#			dc and SEc are the ith column of d 
#			and SE respectively
		j <-  (abs(dc) < eps) & !is.na(SEc) & (dgf != 0)
#			contrast coefficient inside equivalence region, not
#			missing, and with nonzero degrees of freedom
		h <-  (abs(dc) >= eps)  & !is.na(SEc) & (dgf != 0)
#			contrast coefficient outside equivalence region, not 
#			missing, and with nonzero degrees of freedom
		
		####
		
		p[h,i] <- 1 			

#			This if() else statement checks the length of
#			eps, and then computes the distance (in 
#			standard error units) between the per feature 
#			measurements  of difference, dc (in the i-th
#			contrast of course), and the error bars. This
#			forms a confidence interval, whose level
#			provides our p-value. The level is easy to
#			calculate with the Cumulative Distribution
#			Function
if (length(eps)>1)
			{
			k <- ( eps[j] - abs( dc[j] )) / SEc[j]
			} 
		
#				Computing the maximum possible interval
#				which still remains inside the error bars
else
			{
			k <- ( eps - abs( dc[j] )) / SEc[j]
			}
# 				Computing the maximum possible interval 
#				which still remains inside the region of equivalence.
			p[j,i] <- pt(k, dgf[j], lower.tail = FALSE) 
#				Finding the quantile for the given 
#				interval, with the cumulative distribition #				function, Pr(X <= x)
		}
colnames(p) <- paste(colnames(Fit.bayes$contrasts),".tEQ", sep = "")
p
}