# meat storage data from page 41 of the text 

data <- scan(what = list(packaging = "", logbacteria = 0),multi.line=FALSE)
  plastic 7.66  plastic 6.98 plastic 7.80 
  vacuum 5.26  vacuum 5.44  vacuum 5.80 
  mixture 7.41  mixture 7.33  mixture 7.04 
  CO2  3.51  CO2 2.91  CO2 3.66  
   
data


meatpackage <- data.frame(data)
meatpackage 
rm(data)

meatpackage$packaging <- as.factor(meatpackage$packaging)


boxplot(logbacteria ~ packaging, data=meatpackage,main="Meat Packaging data")

anova(lm(logbacteria ~ packaging, data=meatpackage))



 # these expressions calculate the sum of squares for several contrasts

contrast1 <- c(-1,-1,3,-1)
l1hat <- contrast1 %*% tapply(meatpackage$logbacteria,meatpackage$packaging,mean)
denoml1hat <- contrast1^2 / tapply(meatpackage$logbacteria,meatpackage$packaging,length)
ssl1hat <- l1hat^2 / sum(denoml1hat)


contrast2 <- c(-1,-1,0,2)
l2hat <- contrast2 %*% tapply(meatpackage$logbacteria,meatpackage$packaging,mean)
denoml2hat <- contrast2^2 / tapply(meatpackage$logbacteria,meatpackage$packaging,length)
ssl2hat <- l2hat^2 / sum(denoml2hat)


contrast3 <- c(1,-1,0,0)
l3hat <- contrast3 %*% tapply(meatpackage$logbacteria,meatpackage$packaging,mean)
denoml3hat <- contrast3^2 / tapply(meatpackage$logbacteria,meatpackage$packaging,length)
ssl3hat <- l3hat^2 / sum(denoml3hat)


ssl1hat +ssl2hat +ssl3hat     # since sample sizes are equal and contrasts are orthogonal, 
                              # their sum equals the treatment sum of squares


 # Tukey pairwise comparisons

TukeyHSD(aov(logbacteria ~ packaging, data=meatpackage),"packaging" )

plot(TukeyHSD(aov(logbacteria ~ packaging, data=meatpackage),"packaging" ))