###############################################
## Author: Joshua M. Tebbs
## Date: 3 November 2025
## Update: 10 November 2025
## STAT 509 course notes: R Code Chapter 12
###############################################

# Example 12.2
# Page 222-229
# Enter data
a1b1 = c(35.4, 26.9, 25.2, 33.3, 31.5)
a1b2 = c(39.0, 33.3, 41.2, 31.1, 36.3)
a2b1 = c(37.2, 27.9, 35.5, 27.4, 34.3)
a2b2 = c(49.3, 39.8, 39.2, 47.6, 46.7)

# Side-by-side boxplots
# Page 223
boxplot(a1b1,a1b2,a2b1,a2b2,xlab="Treatment combination",
	names=c("a1b1","a1b2","a2b1","a2b2"),ylab="Yield (bushels/plot)",col="lightblue")

# One way classification (ANOVA)
# Page 223
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)
fit = lm(Yield ~ Treatment)
anova(fit)

# Two way ANOVA with interaction
# Page 224
Yield = c(a1b1,a1b2,a2b1,a2b2)
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)

# Interaction plot
# Page 226
interaction.plot(Nitrogen,Phosphorus,Yield,fun=mean,type="l",xtick=TRUE,
	ylab="Mean yield (bushels/plot)",ylim=c(25,50))

# Two way ANOVA without interaction
# Page 227
Yield = c(a1b1,a1b2,a2b1,a2b2)
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.2 = lm(Yield ~ Nitrogen + Phosphorus)
anova(fit.2)

# Individual factor side-by-side boxplots
# Page 228
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")

# Confidence intervals for difference of factor means
# Page 228
# MSE = mean squared residuals from no interaction model
MSE = 21.21
# 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)

####################################################################################################
####################################################################################################

# Example 12.3
# Page 229-232
# Enter 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)

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

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

# Interaction plots
# Page 232
# AC interaction
Temperature = factor(A)
Conc.Form = factor(C)
interaction.plot(Temperature,Conc.Form,filtration,fun=mean,type="l",xtick=TRUE,
	ylab="Filtration rate (gpm/sq.ft)",ylim=c(40,100))
# AD interaction
Temperature = factor(A)
Stir.Rate = factor(D)
interaction.plot(Temperature,Stir.Rate,filtration,fun=mean,type="l",xtick=TRUE,
	ylab="Filtration rate (gpm/sq.ft)",ylim=c(40,100))

# Multiple linear regression
# Page 232
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