# Lecture 25: 2-way anova
# Dennis E. Slice (c) 2013

# Absterge the environment
rm(list=ls())

# Get the pause utility
source("pause.R")

# Load the Sokal and Rohlf rat data
X = read.table("25_anova_2way_SR_Box11_1.RData")
print(X)

# Make names more meaningful
levels(X$fat) = c("fresh","rancid")
levels(X$sex) = c("female","male")

# Test fat
print("aov(data=X,intake~fat)")
print(summary(aov(data=X,intake~fat)))
print("tapply(X$intake,X$fat,mean)")
print(tapply(X$intake,X$fat,mean))

# Test sex
print("aov(data=X,intake~sex)")
print(summary(aov(data=X,intake~sex)))
print("tapply(X$intake,X$sex,mean)")
print(tapply(X$intake,X$sex,mean))

# Test fat and sex
print("aov(data=X,intake~sex+fat)")
print(summary(aov(data=X,intake~sex+fat)))
print("tapply(X$intake,X$fat:X$sex,mean)")
print(tapply(X$intake,X$fat:X$sex,mean))

# Test fat and sex, 2-way
print("aov(data=X,intake~sex+fat+sex:fat)")
print(summary(aov(data=X,intake~sex+fat+sex:fat)))

# Test fat and sex, 2-way
print("aov(data=X,intake~sex*fat)")
# Sincer reusing aov results, should save
# e.g., rat.2way.aov = aov(data=X,intake~sex*fat)
# but not now since I want to show full syntax
print(summary(aov(data=X,intake~sex*fat)))
print("model.tables(aov(data=X,intake~sex*fat))")
print(model.tables(aov(data=X,intake~sex*fat)))

# Graphics setup
ask.orig = par()$ask
par(ask=TRUE)

# interaction plot
interaction.plot(X$fat,X$sex,X$intake,main="aov(intake ~ fat*sex, data=X) interaction plot",xlab="FAT",ylab="mean INTAKE",trace.label="SEX")

# Repeated measures (two-way ANOVA)
girls <- read.table("22_paired_comparisons_SR_Box11.5.RData")
print(girls)
print(t.test(girls$five,girls$six))
print(t.test(girls$six-girls$five))
# Reshape to long format. See ?reshape
girls.long <- reshape( girls,
									 list(names(girls)),
									 direction="long",
									 idvar="girl",
									 v.names="facewidth",
									 times=c("five","six"))
girls.long$time = factor(girls.long$time)

print(girls.long)
with(girls.long,
	interaction.plot(
		time,girl,facewidth,
		main="girl-time interaction plot",
		xlab="AGE",
		ylab="mean FACEWIDTH",
		trace.label="GIRL"	)
)
with(girls.long,
	interaction.plot(
		girl,time,facewidth,
		main="girl-time interaction plot",
		xlab="GIRL",
		ylab="mean FACEWIDTH",
		trace.label="AGE"	
	)
)

print(summary(aov(facewidth ~ time*girl,data=girls.long)))
# Oops, girl (1-15) must be made a factor
girls.long$girl = factor(girls.long$girl)
print(summary(aov(facewidth ~ time*girl,data=girls.long)))
# Hmm, no significance test?
# There are no replicates within girls to compute residual error term
# Must assume no interaction
print(summary(aov(facewidth ~ time+girl,data=girls.long)))
# Compare
print(t.test(girls$six-girls$five))

# Graphics restore
par(ask=ask.orig)