# TODO: Add comment
# 
# Author: E.Korsching 10.9.2009
###############################################################################



levene.bartlett <- function(x, method = "l", measure = "mean", alpha.critical = 0.95)
{
	# variance tests: h0 variance equal
	# x: dataframe, samples in cols, NAs will be removed, degrees of freedom df will be number of cols
	# method: "l": Levene (especially good for deviations from normal distributed data)
	#         "b": Bartlett (also good to detect deviations from normal distributed data)
	# measure (for Levene only): "mean", "median", "trim": 10% trimmed mean
	#      mean provided the best power for symmetric, moderate-tailed, distributions
	#      median performed best when the underlying data followed a X^2_4 (i.e., skewed) distribution
	#      trimmed mean good for Cauchy distribution (i.e., heavy-tailed)
	# erg: 0: var equal, 1: var not equal
	
	#def
	my.length <- function(x)
	{
		length(x[!is.na(x)])
	}
	prae <- function(x, z)
	{
		abs(x - z)
	}
	subterm.z1 <- function(x, z)
	{
		x[1] * (x[2] - z)^2
	}
	#ini
	k <- ncol(x)
	nr <- nrow(x)
	df <- k - 1
	#degrees of freedom
	N.samples <- apply(x, 2, my.length)
	N <- sum(N.samples)
	#
	if(method == "l") {
		# transform x in zz by
		if(measure == "mean") {
			zz <- apply(x, 2, mean, na.rm = T)
		}
		if(measure == "median") {
			zz <- apply(x, 2, median, na.rm = T)
		}
		if(measure == "trim") {
			zz <- apply(x, 2, mean, na.rm = T, trim = 0.1)
		}
		z <- matrix(0, nr, k)
		for(i in 1:k) {
			z[, i] <- prae(x[, i], zz[i])
		}
		z.mean <- apply(z, 2, mean, na.rm = T)
		# per group mean
		Z.mean <- sum(apply(z, 2, sum, na.rm = T), na.rm = T)/sum(
			N.samples, na.rm = T)
		#overall mean
		numerator.sum <- sum(apply(cbind(N.samples, z.mean), 1, 
			subterm.z1, Z.mean), na.rm = T)
		numerator <- (N - k) * numerator.sum
		denominator.sum <- 0
		for(i in 1:k) {
			denominator.sum <- denominator.sum + sum((z[, i] -
				z.mean[i])^2, na.rm = T)
		}
		denominator <- (k - 1) * denominator.sum
		erg <- numerator/denominator
		upper.critical.value.F <- qf(alpha.critical, df, N - k)
		#critical value with df1=k-1 and df2=N-k
		cat("\n h0: var homogenous")
		cat("\n W = ", round(erg, 2))
		cat("\n Upper critical value of F distribution (alpha = ",
			1 - alpha.critical, ") = ", round(
			upper.critical.value.F, 2))
		cat("\n h0 acceptance intervall ( 0 , ", alpha.critical, " )")
		if(erg > upper.critical.value.F) {
			cat("\n h0 is rejected   (W > critical value)\n")
			erg <- 1
		}
		else {
			cat("\n h0 accepted   (W <= critical value)\n")
			erg <- 0
		}
	}
	if(method == "b") {
		si2 <- apply(x, 2, var, na.method = "available")
		# var of groups
		sp2 <- 0
		for(i in 1:k) {
			sp2 <- sp2 + (((N.samples[i] - 1) * si2[i])/(N - k))
		}
		numerator.right <- 0
		for(i in 1:k) {
			numerator.right <- numerator.right + ((N.samples[
				i] - 1) * log(si2[i]))
		}
		numerator <- (N - k) * log(sp2) - numerator.right
		denominator.sum <- 0
		for(i in 1:k) {
			denominator.sum <- denominator.sum + (1/(N.samples[
				i] - 1))
		}
		denominator.right <- denominator.sum - (1/(N - k))
		denominator.mid <- 1/(3 * (k - 1))
		denominator <- 1 + denominator.mid * denominator.right
		erg <- numerator/denominator
		upper.critical.value.chi.sqr <- qchisq(alpha.critical, df)
		#critical value
		cat("\n h0: var homogenous")
		cat("\n T = ", round(erg, 2))
		cat("\n Upper critical value of chi square distribution (alpha = ",
			1 - alpha.critical, ") = ", round(
			upper.critical.value.chi.sqr, 2))
		cat("\n h0 acceptance intervall ( 0 , ", alpha.critical, " )")
		if(erg > upper.critical.value.chi.sqr) {
			cat("\n h0 is rejected   (T > critical value)\n")
			erg <- 1
		}
		else {
			cat("\n h0 accepted   (T <= critical value)\n")
			erg <- 0
		}
	}
	return(erg)
}


#levene.bartlett(x=cbind(rnorm(n=100,mean=1,sd=1),runif(n=100,min=0,max=1)), method = "l", measure = "mean", alpha.critical = 0.95)
#levene.bartlett(x=cbind(rnorm(n=100,mean=1,sd=1),rnorm(n=100,mean=0,sd=1)), method = "l", measure = "mean", alpha.critical = 0.95)