# ANOVA
# Dennis E. Slice (c) 2012
# Null hypothesis data

# Clear the cobwebs
rm(list = ls())

# Function to read integer.
# Borrowed and modified from internet:
# http://stackoverflow.com/questions/5974967/what-is-the-correct-way-to-ask-for-user-input-in-an-r-program 
getInt <- function(s = "Enter a postive integer:") {
	n = -1
	while (n < 1) {
		n <- readline(s)
		n <- ifelse(grepl("\\D", n), -1, as.integer(n))
		if (is.na(n)) {
			break
		} # breaks when hit enter
	}
	# Function returns last result of last statement.
	# Here, just n...
	n
}

# We could put this in a loop and play around all day long.
# Some parameters - change for different plots
nGroups <- getInt("nGroups: ")
# Equal n = balanced design
# Unequal n = unbalanced, requires adj. of comp. and interpretation 
n = getInt("n (per sample): ")
nTotal = n * nGroups
# Distribution parameters - something other than 0,1
mu = 30
sigma = 5

# Group factor
g = rep(1:nGroups, n) # 1 2 3...1 2 3...
print(g)
g = factor(g)
print(g)
y = rnorm(nTotal, mean = mu, sd = sigma)
print(y)

# boxplot y by groups
# range: mean +-4sd = 99.994% of all observations
boxplot(	
	y ~ g,
	ylim = c(mu - 4 * sigma, mu + 4 * sigma), 
	main=paste("H0 - n=",n,", ",nGroups," groups",sep=""),
	xlab="GROUP",
	ylab="Y"
)

# Analysis of Variance
result <- aov(y~g)
print("Analysis of Variance")
print(result)
print("Summary")
print(summary(result))