################################
## Lecture 3 Linear models II ##
## Two-Way Anova ##
################################
#############
## SLIDE 3 ## an alternative way to specify the full path to your file
#############
dir <- "C:\\My Documents\\Teaching\\Modules\\Adv. Research Skills\\Datasets\\"
my.file <- paste(dir,"2wayAnova.csv", sep="")
dataset <- read.table(my.file, header=T, sep=",")
attach(dataset)
names(dataset)
str(dataset)
dataset
tapply(plasma, list(sex, hormone), length)
#############
## SLIDE 4 ## boxplot not shown on slide, but given here if you want to try.
#############
boxplot(plasma ~ hormone, col="grey95", cex.lab=1.2, cex.axis=1.2,
ylab="Plasma Concentration (mg/100 ml)", xlab="Hormone Treatment", medcol="red")
stripchart(plasma ~ hormone, col=c("orange","blue"), cex.lab=1.2, cex.axis=1.2,
ylab="Plasma Concentration (mg/100 ml)", xlab="Hormone Treatment",
vertical=T, method="jitter", pch=16)
averages <- tapply(plasma, list(sex, hormone), mean)
segments(x0=0.8, y0=apply(averages,2,mean)[1],
x1=1.2, y1=apply(averages,2,mean)[1], col="orange", lwd=2)
segments(x0=1.8, y0=apply(averages,2,mean)[2],
x1=2.2, y1=apply(averages,2,mean)[2], col="blue", lwd=2)
h.ave <- tapply(plasma, hormone, mean)
h.ave
boxplot(plasma ~ sex, col="grey95", cex.lab=1.2, cex.axis=1.2,
ylab="Plasma Concentration (mg/100 ml)", xlab="Sex", medcol="red")
stripchart(plasma ~ sex, col=c("red","green"), cex.lab=1.2, cex.axis=1.2,
ylab="Plasma Concentration (mg/100 ml)", xlab="Sex", vertical=T,
method="jitter", pch=16)
segments(x0=0.8, y0=apply(averages,1,mean)[1],
x1=1.2, y1=apply(averages,1,mean)[1], col="red", lwd=2)
segments(x0=1.8, y0=apply(averages,1,mean)[2],
x1=2.2, y1=apply(averages,1,mean)[2], col="green", lwd=2)
s.ave <- tapply(plasma, sex, mean)
s.ave
#############
## SLIDE 5 ##
#############
h.ave
gd.ave <- mean(plasma); gd.ave
summary(aov(plasma ~ hormone))
SS.h <- 10*(21.825-13.50)^2 +
10*(21.825-30.15)^2
TSS <- sum((plasma-gd.ave)^2)
#############
## SLIDE 6 ##
#############
s.ave
gd.ave
summary(aov(plasma ~ sex))
SS.s <- 10*(21.825-23.70)^2 +
10*(21.825-19.95)^2
TSS <- sum((plasma-gd.ave)^2)
#############
## SLIDE 7 ##
#############
summary(aov(plasma ~ hormone))
summary(aov(plasma ~ sex))
summary(aov(plasma ~ hormone + sex))
#############
## SLIDE 8 ##
#############
summary(aov(plasma ~ hormone + sex))
summary(aov(plasma ~ hormone * sex))
#############
## SLIDE 9 ##
#############
gd.ave
h.diff <- h.ave - gd.ave
s.diff <- s.ave - gd.ave
predicted <- predict(aov(plasma~hormone+sex))
##############
## SLIDE 10 ##
##############
pred.ave <- tapply(predicted, list(sex, hormone), mean)
pred.ave
### averages <- tapply(plasma, list(sex, hormone), mean)
### already created for slide 4, so not necessary to repeat it again.
averages
sum((averages-pred.ave)^2)*5
##############
## SLIDE 11 ##
##############
averages
barplot(averages, beside=T, col=c("red","green"),
xlab="Hormone treatment", ylab="Plasma Concentration (mg/100 ml)",
cex.lab=1.2, cex.axis=1.2, cex.names=1.2 )
labs <- c("female","male") ## labels to be used in legend()
legend(locator(1), labs, fill=c("red","green"), cex=1.2 )
## locator(1) allows you to click on the graph the precise location
## of the upper left corner of your legend box. If it doesn't display
## in the right location, you need to start all over again with barplot
pred.ave
barplot(pred.ave, beside=T, col=c("red2","green3"),
xlab="Hormone treatment", ylab="Plasma Concentration (mg/100 ml)",
cex.lab=1.2, cex.axis=1.2, cex.names=1.2 )
legend(locator(1), labs, fill=c("red2","green3"), cex=1.2 )
##############
## SLIDE 12 ##
##############
interaction.plot(hormone, sex, plasma, col=c("red","green2"), lwd=2,
cex.lab=1.2, cex.axis=1.2, legend=F )
legend(locator(1), legend=c("female", "male"), cex=1.2,
lty=c(2,1), lwd=2, col=c("red","green2") )
interaction.plot(sex, hormone, plasma, col=c("orange","blue"), lwd=2,
cex.lab=1.2, cex.axis=1.2, legend=F )
legend(locator(1), legend=c("No Hormone", "Hormone"), cex=1.2,
lty=c(2,1), lwd=2, col=c("orange","blue") )
##############
## SLIDE 13 ##
##############
mod <- lm(plasma ~ sex + hormone)
summary(mod)
pred.ave
##############
## SLIDE 14 ##
##############
mod2 <- lm(plasma ~ sex * hormone)
summary(mod2)
averages
##############
## SLIDE 15 ##
##############
shapiro.test(mod$residuals)
plot(mod, which=2, lwd=2, pch=16, cex.axis=1.2,
cex.lab=1.2, cex.sub=1.2, cex.caption=1.2)
qqline(rstandard(mod), lwd=2, lty=1, col="red")
## adds a red line for extra clarity
##############
## SLIDE 16 ##
##############
plot(mod, which=1, lwd=2, pch=16, cex.axis=1.2,
cex.lab=1.2, cex.sub=1.2, cex.caption=1.2)
plot(mod, which=3, lwd=2, pch=16, cex.axis=1.2,
cex.lab=1.2, cex.sub=1.2, cex.caption=1.2)