#########################################################
#########################################################
##########                                      #########
##########	R Scripts for MATH 204 Winter 2013  #########
##########                                      #########
#########################################################
#########################################################

#Randomized block design
#Soil data
#Table of data
xtabs(sulphur~solvent+soil,data=Soil)

summary(aov(sulphur~1,data=Soil))
summary(aov(sulphur~solvent+soil,data=Soil))

require(car)
#Not run
#leveneTest(sulphur~solvent+soil,center='mean',data=Soil)
#leveneTest(sulphur~solvent+soil,center='median',data=Soil)

#The need for blocking - example

Y<-c(-20.88,-4.76,-0.46,10.78,-10.47,
-15.75,0.11,5.64,16.98,-6.03,
-8.62,6.20,10.42,20.05,-1.29)

Group<-gl(3,5,labels=c("Group1","Group2","Group3"))
Block<-gl(5,1,15,labels=c(1:5))

library(gplots)
plotmeans(Y~Group,ylim=range(-40,40))
stripchart(Y~Group,vertical=T,pch=4,cex=0.75,lwd=2,add=T)

require(car)
leveneTest(Y~Group,center='mean')
leveneTest(Y~Group,center='median')

#Ignoring the blocking in ANOVA
summary(aov(Y~1))
summary(aov(Y~Group))

#Acknowledging the blocking in ANOVA
summary(aov(Y~Group+Block))

blocklabel<-paste(rep("Block",5),c(1:5))
plotmeans(Y~Group,ylim=range(-40,40))
points(rep(0.9,5),Y[Group=="Group1"],pch=c(1:5))
points(rep(1.9,5),Y[Group=="Group2"],pch=c(1:5))
points(rep(2.9,5),Y[Group=="Group3"],pch=c(1:5))
legend(0.5,40,blocklabel,pch=(1:5))

##################################################################
##################################################################
#Randomized block designs with replication
#Batteries data

#Table of data
#Counts of replicates
xtabs(~temp+material,data=Batteries)
#Means
xtabs(life~temp+material,data=Batteries)/xtabs(~temp+material,data=Batteries)

summary(aov(life~1,data=Batteries))
summary(aov(life~temp+material+temp:material,data=Batteries))
summary(aov(life~temp*material,data=Batteries))
leveneTest(life~temp*material,data=Batteries,center='mean')
leveneTest(life~temp*material,data=Batteries,center='median')


##################################################################
##################################################################
#Randomized complete factorial designs
#Butterfat data

leveneTest(butterfat~breed*age,data=Butterfat,center='mean')
leveneTest(butterfat~breed*age,data=Butterfat,center='median')

m3<-aov(butterfat~breed*age,data=Butterfat)
summary(m3)
m2<-aov(butterfat~breed+age,data=Butterfat)
summary(m2)
m1<-aov(butterfat~breed,data=Butterfat)
summary(m1)
m0<-aov(butterfat~1,data=Butterfat)
anova(m0,m1,m2,m3)

#Set-up means plot
plotmeans(0*butterfat~breed,data=Butterfat,ylim=range(3,7),ylab="Butterfat",pch=19)
plotmeans(butterfat~breed,data=Butterfat,subset=(age=="Two year old"),add=TRUE,n.label=FALSE,col="blue",ccol="blue",ylab="Butterfat",pch=19)
plotmeans(butterfat~breed,data=Butterfat,subset=(age=="Mature"),add=T,n.label=FALSE,lty=2,pch=15,col="red",barcol="red",ccol="red")
legend(0.5,7,c("Two year old","Mature"),pch=c(19,15),lty=c(1,2),col=c("blue","red"))


#Set-up means plot
plotmeans(0*butterfat~age,data=Butterfat,ylim=range(3,7),ylab="Butterfat",pch=19)
plotmeans(butterfat~age,data=Butterfat,subset=(breed=="Ayreshire"),add=TRUE,n.label=FALSE,ylim=range(3,6),col="blue",ccol="blue",ylab="Butterfat",pch=19)
plotmeans(butterfat~age,data=Butterfat,subset=(breed=="Canadian"),add=T,n.label=FALSE,lty=2,pch=15,col="red",barcol="red",ccol="red")
plotmeans(butterfat~age,data=Butterfat,subset=(breed=="Guernsey"),add=T,n.label=FALSE,lty=3,pch=17,col="green",barcol="green",ccol="green")
plotmeans(butterfat~age,data=Butterfat,subset=(breed=="Holstein"),add=T,n.label=FALSE,lty=4,pch=4,col="darkgreen",barcol="darkgreen",ccol="darkgreen")
plotmeans(butterfat~age,data=Butterfat,subset=(breed=="Jersey"),add=T,n.label=FALSE,lty=5,pch=3,col="purple",barcol="purple",ccol="purple")
legend(1.5,7,unique(Butterfat$breed),pch=c(19,15,17,4,3),lty=c(1:5),col=c("blue","red","green","darkgreen","purple"))


leveneTest(m1,center='mean')
leveneTest(m1,center='median')
leveneTest(m3,center='mean')
leveneTest(m3,center='median')

##################################################################
##################################################################
#Lyrics data

m3<-aov(SCORE~SONG*POOL,data=Lyrics)
summary(m3)
m2<-aov(SCORE~SONG+POOL,data=Lyrics)
summary(m2)
m1<-aov(SCORE~SONG,data=Lyrics)
summary(m1)
m0<-aov(SCORE~1,data=Lyrics)
anova(m0,m1,m2,m3)

#Set-up means plot
plotmeans(0*SCORE~SONG,data=Lyrics,ylim=range(2.5,4),ylab="SCORE",pch=19)
plotmeans(SCORE~SONG,subset=(POOL=="PSYCHOL "),data=Lyrics,ylim=range(2.5,4),n.label=FALSE,pch=19,col="blue",ccol="blue",)
plotmeans(SCORE~SONG,subset=(POOL=="VOLUNTEE"),data=Lyrics,add=T,lty=2,pch=15,n.label=FALSE,col="red",barcol="red",ccol="red")
legend(0.5,4,c("PSYCHOL ","VOLUNTEE"),pch=c(19,15),lty=c(1,2),col=c("blue","red"))

#Set-up means plot
plotmeans(0*SCORE~POOL,data=Lyrics,ylim=range(2.5,4),ylab="SCORE",pch=19)
plotmeans(SCORE~POOL,subset=(SONG=="NON-VIO     "),data=Lyrics,ylim=range(2.5,4),n.label=FALSE,pch=19,col="blue",ccol="blue",)
plotmeans(SCORE~POOL,subset=(SONG=="VIOLENT     "),data=Lyrics,add=T,lty=2,pch=15,n.label=FALSE,col="red",barcol="red",ccol="red")
legend(1.5,4,c("NON-VIO ","VIOLENT"),pch=c(19,15),lty=c(1,2),col=c("blue","red"))


leveneTest(m1,center='mean')
leveneTest(m1,center='median')
leveneTest(m3,center='mean')
leveneTest(m3,center='median')

plotmeans(Y~Group,ylim=range(-40,40))


##################################################################
##################################################################
#Gravel data

m3<-aov(amount~shift*quarry,data=Gravel)
summary(m3)
m2<-aov(amount~shift+quarry,data=Gravel)
summary(m2)
m1<-aov(amount~quarry,data=Gravel)
m0<-aov(amount~1,data=Gravel)
anova(m0,m1,m2,m3)

#Set-up means plot
plotmeans(0*amount~quarry,data=Gravel,ylim=range(5,35),ylab="Amount",pch=19)
plotmeans(amount~quarry,subset=(shift=="AM"),data=Gravel,n.label=FALSE,pch=19,col="blue",ccol="blue",add=T)
plotmeans(amount~quarry,subset=(shift=="PM"),data=Gravel,add=T,lty=2,pch=15,n.label=FALSE,col="red",barcol="red",ccol="red")
legend(3.5,35,c("AM","PM"),pch=c(19,15),lty=c(1,2),col=c("blue","red"))


plotmeans(0*amount~shift,data=Gravel,ylim=range(5,35),ylab="Amount",pch=19)
plotmeans(amount~shift,subset=(quarry=="A"),data=Gravel,n.label=FALSE,pch=19,col="blue",ccol="blue",add=T)
plotmeans(amount~shift,subset=(quarry=="B"),data=Gravel,add=T,lty=2,pch=15,n.label=FALSE,col="red",barcol="red",ccol="red")
plotmeans(amount~shift,subset=(quarry=="C"),data=Gravel,add=T,n.label=FALSE,lty=3,pch=17,col="green",barcol="green",ccol="green")
plotmeans(amount~shift,subset=(quarry=="D"),data=Gravel,add=T,n.label=FALSE,lty=4,pch=4,col="darkgreen",barcol="darkgreen",ccol="darkgreen")
legend(1.5,35,c("A","B","C","D"),pch=c(19,15,17,4),lty=c(1:4),col=c("blue","red","green","darkgreen"))

leveneTest(m1,center='mean')
leveneTest(m1,center='median')
leveneTest(m3,center='mean')
leveneTest(m3,center='median')

##################################################################
##################################################################
#Models: interactions without main effects

summary(aov(life~material*temp,data=Batteries))
summary(aov(life~material+material:temp,data=Batteries))
summary(aov(life~temp+material:temp,data=Batteries))