# Two-Way ANOVA is a technique for studying the relationship between a quantitative (numerical) dependent variable and two qualitative (categorical) independent variables.
# We want to assess if there is any interaction between them.  
# In one-way anova, there is only one categorical variable (or factor).
# In the donut example, the categorical variable was the type of fat.
# In the smokers example, the categorical variable was they type of smoker

# Assumptions:
# Population is normal; means may differ but sd same.
# the data points are relevant with respect to the scientific question under investigation;
# the mean of the response variable is influenced additively (if not interaction term) and linearly by the factors;
# the errors are independent;
# the errors have the same variance;
# the errors are normally distributed.

# http://courses.statistics.com/software/R/Rtwoway.htm
# what are the ind variables?  dep variables?

# Q: Is there any relationship between number of shots made out of 50, 
# and time of day or shoes worn?
# A: p-values are large, suggests no interaction.

basketball=read.table("http://courses.statistics.com/software/data/baskball.txt", header=T)

# 2X2 Table; rows are time; columns are shoes.
# cells are the sample sizes (or some other data)
attach(basketball)
tapply(Made,Time,mean)
tapply(Made,Shoes,mean)

# Why is one names int and noint?
int <- aov(Made ~ Time*Shoes)

# The two below return the same thing.
summary(int)
anova(int)

noint <- aov(Made~Time + Shoes)
summary(noint)

# lm and aov are not exactly the same but the following are the same
anova(aov(Made~Time+Shoes))
anova(lm(Made~Time+Shoes))
summary(aov(Made~Time+Shoes))

boxplot(Made ~ Time)
boxplot(Made ~ Shoes)

# Case Study 13.2.1
# there are k=3 treatment levels and b=5 blocks
# residuals=errors in case study 

block=rep(c("A","B","C","D","E"),3)
therapy=c(rep("cd",5),rep("dp",5),rep("lm",5))
scores=c(8,11,9,16,24,2,1,12,11,19,-2,0,6,2,11)
hat=data.frame(block,therapy,scores)
anova(aov(scores~therapy+block,data=hat))

# therapy: F=(260.93/(3-1))/(68.40/((5-1)*(3-1)))
# blocks:  F=(438/(5-1))/(68.40/((5-1)*(3-1)))





# Sunfish data; glm generalized linear model
tapply(sunfish$spinybase,sunfish$species,mean)
tapply(sunfish$length,sunfish$species,mean)
boxplot(length~species,data=sunfish)
mod1 <- aov(length~spinybase*species, data=sunfish)
mod2 <- aov(length~spinybase+species, data=sunfish)
anova(mod1,mod2)