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

# Example 12.2
# Corn yield data
# Page 199-200
# Enter the data
a1b1 = c(35,26,25,33,31)
a1b2 = c(39,33,41,31,36)
a2b1 = c(37,27,35,27,34)
a2b2 = c(49,39,39,47,46)

# Page 200
# One way ANOVA/Overall F test
# Create response with all data
yield = c(a1b1,a1b2,a2b1,a2b2)
# Create treatment indicator variable
treatment = c(
	rep("a1b1",length(a1b1)), 
	rep("a1b2",length(a1b2)), 
	rep("a2b1",length(a2b1)), 
	rep("a2b2",length(a2b2))
	)
treatment = factor(treatment)
# One way analysis
anova(lm(yield~treatment))

# Page 201
# Construct side by side boxplots
boxplot(a1b1,a1b2,a2b1,a2b2,range=1,xlab="Treatment combination",names=c("a1b1","a1b2","a2b1","a2b2"),ylab="Yield (bushels/plot)",col="lightblue")

# Page 202
# Two way ANOVA with interaction
nitrogen = c(
	rep("a1",length(a1b1)),
	rep("a1",length(a1b2)),
	rep("a2",length(a2b1)),
	rep("a2",length(a2b1))
	)
phosphorus = c(
	rep("b1",length(a1b1)),
	rep("b2",length(a1b2)),
	rep("b1",length(a2b1)),
	rep("b2",length(a2b1))	
	)
nitrogen = factor(nitrogen)
phosphorus = factor(phosphorus)
fit = lm(yield~nitrogen+phosphorus+nitrogen*phosphorus)
anova(fit)

# Page 203
# Plot F(1,16) pdf
f = seq(0,8,0.01)
pdf = df(f,1,16)
plot(f,pdf,type="l",lty=1,xlab="F",ylab="PDF")
abline(h=0)
abline(v=0,lty=2)
# Add points
points(x=2.25,y=0,pch=4,cex=1.5)

# Page 204
# Interaction plot
interaction.plot(nitrogen,phosphorus,yield,fun=mean,type="l",ylab="Mean yield (bushels/plot)",ylim=c(25,50))

# Page 206
# Two way ANOVA with no interaction
fit.2 = lm(yield~nitrogen+phosphorus)
anova(fit.2)

# Page 206
# Confidence intervals for difference of factor means
# MSE = mean squared residuals from no interaction model
MSE = 21.47
# Aggregate Factor A data
nitrogen.1 = c(a1b1,a1b2)
nitrogen.2 = c(a2b1,a2b2)
# Aggregate Factor B data
phosphorus.1 = c(a1b1,a2b1)
phosphorus.2 = c(a1b2,a2b2)
# Nitrogen CI
est = mean(nitrogen.1)-mean(nitrogen.2)
moe = qt(0.975,17)*sqrt(MSE*0.2)
c(est-moe,est+moe)
# Phospherus CI
est = mean(phosphorus.1)-mean(phosphorus.2)
moe = qt(0.975,17)*sqrt(MSE*0.2)
c(est-moe,est+moe)

# Page 207
# Individual factor side by side boxplots
# Plots were constructed separately
boxplot(nitrogen.1,nitrogen.2,xlab="",names=c("Nitrogen.1","Nitrogen.2"),ylab="Yield (bushels/plot)",ylim=c(20,50),col="lightblue")
boxplot(phosphorus.1,phosphorus.2,xlab="",names=c("Phosphorus.1","Phosphorus.2"),ylab="Yield (bushels/plot)",ylim=c(20,50),col="lightblue")


##############
# Example 12.3
##############

# Filtration rate data
# Page 208
# Enter the data
filtration = c(45,71,48,65,68,60,80,65,43,100,45,104,75,86,70,96)
A = c("a1","a2","a1","a2","a1","a2","a1","a2","a1","a2","a1","a2","a1","a2","a1","a2")
B = c("b1","b1","b2","b2","b1","b1","b2","b2","b1","b1","b2","b2","b1","b1","b2","b2")
C = c("c1","c1","c1","c1","c2","c2","c2","c2","c1","c1","c1","c1","c2","c2","c2","c2")
D = c("d1","d1","d1","d1","d1","d1","d1","d1","d2","d2","d2","d2","d2","d2","d2","d2")
A = factor(A)
B = factor(B)
C = factor(C)
D = factor(D)

# Page 210
# Fit full model
fit = lm(filtration~A*B*C*D)
anova(fit)

# Page 211
# Fit smaller model
fit = lm(filtration~A+C+D+A*C+A*D)
anova(fit)

# Page 211
# Interaction plots
# Plots were constructed separately
# AC Interaction plot
Factor.A = factor(A)
Factor.C = factor(C)
interaction.plot(Factor.A,Factor.C,filtration,fun=mean,type="l",ylab="Filtration rate (gal/hr)",ylim=c(40,100))
# AD Interaction plot
Factor.A = factor(A)
Factor.D = factor(D)
interaction.plot(Factor.A,Factor.D,filtration,fun=mean,type="l",ylab="Filtration rate (gal/hr)",ylim=c(40,100))

# Page 212
# Multiple linear regression model
filtration = c(45,71,48,65,68,60,80,65,43,100,45,104,75,86,70,96)
temp = c(-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1,-1,1)
conc = c(-1,-1,-1,-1,1,1,1,1,-1,-1,-1,-1,1,1,1,1)
stir = c(-1,-1,-1,-1,-1,-1,-1,-1,1,1,1,1,1,1,1,1)
fit = lm(filtration~temp+conc+stir+temp:conc+temp:stir)
fit