# This function performs randomization tests for ANOVAs of any design.

# Inputs are your ANOVA model and the number of randomizations you want to do, s.
	rand.anova <- function(my.formula, my.data, resp, s) {

# Do the standard ANOVA, calculate number of tests to be done.
	emp.aov <- aov(my.formula,my.data)
	emp.aov.summary <- summary(emp.aov)

	ntab <- NULL
	ntab <- length(emp.aov.summary)
	ntest <- NULL
	total.tests <- NULL
	if (ntab==1)
		{ntest <- dim(emp.aov.summary[[1]])[1] - 1} else {
		for (i in 1:length(emp.aov.summary)) {
			if (dim(emp.aov.summary[[i]][[1]])[2] < 5)
				{ntest[i] <- 0} else {
				ntest[i] <- dim(emp.aov.summary[[i]][[1]])[1] - 1
			}}
	}
	total.tests <- sum(ntest)

# Extract the empirical F-stats into a matrix where each row corresponds to an F table and each column to an F stat.
	fstats <- NULL
	if (ntab==1)
		{fstats <- emp.aov.summary[[1]][1:ntest,4]} else {
		fstats <- rep(NA,ntab*max(ntest))
		dim(fstats) <- c(ntab,max(ntest))
		for (i in 1:ntab) {
			if (ntest[i]==0) {next}
			if (ntest[i]==1) {
				fstats[i,1] <- emp.aov.summary[[i]][[1]][1,4]} else {
				for (j in 1:ntest[i]) {
					fstats[i,j] <- emp.aov.summary[[i]][[1]][j,4]
				}
			}
		}
	}

# Make a matrix that will fit the null distribution of F-values for each test into columns.
	Fnull <- rep(NA,length(fstats)*s)
	fstats <- as.matrix(fstats)
	dim(Fnull) <- c(dim(fstats)[2],dim(fstats)[1],s)

# Randomize the response variable and do the ANOVA on teh randomized variable s times.
# Take the F-values from each randomized ANOVA and place them into a row of Fnull.
	for (i in 1:s) {
		my.data$perm <- sample(resp,length(resp),replace=FALSE)
		perm.formula <- update(my.formula,perm ~ .)
		perm.aov.rep <- summary(aov(perm.formula,my.data))

		if (ntab==1) {
			for (k in 1:ntest) {
				Fnull[1,k,i] <- perm.aov.rep[[1]][k,4]}} else {
			for (j in 1:ntab) {
				if (ntest[j]==0) {next}
				if (ntest[j]==1) {
					Fnull[1,j,i] <- perm.aov.rep[[j]][[1]][1,4]} else {
				for (k in 1:ntest[j]) {
					Fnull[k,j,i] <- perm.aov.rep[[j]][[1]][k,4]
				}
				}
			}
		}
	}

# Make a matrix to house the p-values from the randomization tests and calculate a one-sided p-value 
# for each effect.
	pval <- rep(NA,ntab*max(ntest+1))
	dim(pval) <- c(ntab,max(ntest+1))
	if (ntab==1) {
		for (j in 1:ntest) {
			pval[1,j] <- sum(Fnull[1,j,] > fstats[j,1]) / s}
	} else {
		for (i in 1:ntab) {
			if (ntest[i]==0) {next}
			if (ntest[i]==1) {
				pval[i,1] <- sum(Fnull[1,i,] > fstats[i,1]) / s}
			else {
				for (k in 1:ntest[i]) {
					pval[i,k] <- sum(Fnull[k,i,] > fstats[i,k]) / s}
				}}
	}

# Compile results and print.
	results <- emp.aov.summary
	if (ntab==1) {
		results[[1]] <- cbind(results[[1]],rep(NA,ntest+1))
		colnames(results[[1]])[6] <- "Prand"
		results[[1]][1:(ntest+1),6] <- t(pval)
	} else {
		for (i in 1:ntab) {
			if (ntest[i]==0) {next}
			if (ntest[i]==1) {
				results[[i]][[1]] <- cbind(results[[i]][[1]],rep(NA,2))
				colnames(results[[i]][[1]])[6] <- "Prand"
				results[[i]][[1]][1:(ntest[i]+1),6] <- pval[i,1:(ntest[i]+1)]
			} else {
				results[[i]][[1]] <- cbind(results[[i]][[1]],rep(NA,ntest[i]+1))
				colnames(results[[i]][[1]])[6] <- "Prand"
				results[[i]][[1]][1:(ntest[i]+1),6] <- pval[i,1:(ntest[i]+1)]
		}
	}}
	results
}