# April 24, 2008
# ANOVA Demonstrations

######################
data(coagulation)
attach(coagulation)
names(coagulation)
is.factor(diet)
boxplot(coag ~ diet, data = coagulation)

# Use lm to set up a linear model. Syntax has the following form:
#	lm(dependent variable ~ independent variable, data = dataframe)
#

coag.lm = lm(coag ~ diet, data = coagulation)
anova(coag.lm)

# if we wanted to get p value ourselves
pf(13.571, 3, 20, lower.tail = FALSE)

# or we could do it all at once
anova(lm(coag ~ diet, data = coagulation))

# See what happens when we don't make sure our categorical variable is a factor
anova(lm(coag ~ as.numeric(diet), data = coagulation))

# checking assumptions
#par(mfrow = c(2,2))
#plot(aov(coag ~ diet))
qqnorm(coag.lm$res)
qqline(residuals(coag.lm))
plot(coag.lm$fitted, coag.lm$res,xlab="Fitted",ylab="Residuals") 
shapiro.test(residuals(coag.lm))
bartlett.test(coagulation$coag, coagulation$diet)

# To find which levels are significantly different
tukey <- TukeyHSD(aov(coag ~ diet, data = coagulation), conf.level = 0.95)
plot(tukey)

###### 2 Way ANOVA

library(rats)
attach(rats)
# 48 rats were subjected to three poisons and 4 treatments. The outcome variable was survival time in tens of hours
rats
par(mfrow=c(2,1))
plot(time ~ treat + poison, data=rats)

# write out the full model including main effects for each factor as well as interaction
g <- lm(time ~ poison + treat + poison:treat, data = rats)
anova(g)

# short hand notation
g <- lm(time ~ poison*treat, data = rats)
anova(g)

# checking interaction effects
interaction.plot(rats$poison,rats$treat,rats$time)

# check assumptions
qqnorm(g$res)
plot(g$fitted,g$res,xlab="Fitted",ylab="Residuals") 

# transform
g <- lm(1/time ~ poison*treat,rats) 
qqnorm(g$res)
anova(g)
g <- lm(1/time ~ poison+treat, rats)
anova(g)

###########
# nonparametric analog
kruskal.test(coagulation$coag, coagulation$diet)