###############################################
## Author: Joshua M. Tebbs
## Date: 11 December 2011
## Update: 31 December 2017
## STAT 509 course notes: R Code
## Chapter 9
###############################################

# Example 9.1
# Mortar strength data
# Page 136-137
OCM = c(51.45,42.96,41.11,48.06,38.27,38.88,42.74,49.62)
PIM = c(64.97,64.21,57.39,52.79,64.87,53.27,51.24,55.87,61.76,67.15)
RM = c(48.95,62.41,52.11,60.45,58.07,52.16,61.71,61.06,57.63,56.80)
PCM = c(35.28,38.59,48.64,50.99,51.52,52.85,46.75,48.31)
# Create side-by-side boxplots
boxplot(OCM,PIM,RM,PCM,xlab="",names=c("OCM","PIM","RM","PCM"),ylab="Strength (MPa)",col="lightblue")

# Perform F test/Create ANOVA table
# Page 143
# Create response with all data
strength = c(OCM,PIM,RM,PCM)
# Create a treatment indicator variable
mortar.type = c(
	rep("OCM",length(OCM)), 
 	rep("PIM",length(PIM)), 
 	rep("RM",length(RM)), 
 	rep("PCM",length(PCM))
	)
mortar.type = factor(mortar.type)
# Analysis of variance
anova(lm(strength~mortar.type))

# Plot F(3,32) pdf
# Page 143
f = seq(0,18,0.01)
plot(f,df(f,3,32),type="l",lty=1,xlab="F",ylab="PDF")
abline(h=0)
# Add points
points(x=16.848,y=0,pch=4,cex=1.5)

# Follow-up analysis
# Page 149-150
# Tukey confidence intervals
TukeyHSD(aov(lm(strength~mortar.type)),conf.level=0.95)