###############################################
## Author: Joshua M. Tebbs
## Date: 10 September 2025
## Update: 17 September 2025
## STAT 509 course notes: R Code Chapter 9
###############################################

# Use the following command as necessary
par(mar=c(4,4.5,4,3.5)) # adjust plotting region to avoid cutting off vertical axis label 

# Figure 9.1
# Page 155
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)
boxplot(OCM,PIM,RM,PCM,xlab="",names=c("OCM","PIM","RM","PCM"),ylab="Strength (in psi)",
	ylim=c(30,70),col="lightblue")


# Figure 9.2
# Page 157
x = seq(35,95,0.01)
# Left
plot(x,dnorm(x,65,3),type="l",lty=1,xlab="",ylab=expression(f[X](x)),xaxt='n',ylim=c(0,0.175),cex.lab=1.25)
axis(side = 1, at = c(40,50,60,70,80,90), labels = FALSE, tck = -0.015)
lines(x,dnorm(x,65,3),lty=4)
lines(x,dnorm(x,65,3),lty=8)
lines(x,dnorm(x,65,3),lty=8)
abline(h=0)
# Add legend
legend(75,0.175,lty = c(1,4,8,16), 
	c(
	expression(paste("Population 1")),
	expression(paste("Population 2")),
	expression(paste("Population 3")),
	expression(paste("Population 4"))
	))
# Right
plot(x,dnorm(x,50,3),type="l",lty=1,xlab="",ylab=expression(f[X](x)),xaxt='n',ylim=c(0,0.175),cex.lab=1.25)
axis(side = 1, at = c(40,50,60,70,80,90), labels = FALSE, tck = -0.015)
lines(x,dnorm(x,55,3),lty=4)
lines(x,dnorm(x,65,3),lty=8)
lines(x,dnorm(x,80,3),lty=8)
abline(h=0)
# Add legend
legend(75,0.175,lty = c(1,4,8,16), 
	c(
	expression(paste("Population 1")),
	expression(paste("Population 2")),
	expression(paste("Population 3")),
	expression(paste("Population 4"))
	))


# Figure 9.3
# Page 160
q = seq(0,6,0.01)
plot(q,df(q,6,12),type="l",lty=1,xlab="",ylab="",xaxp=c(1,1,1),cex.lab=1.25)
abline(h=0)


# Example 9.1
# Data analysis
# Page 161
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)
# Concatenate all the data
Strength = c(OCM,PIM,RM,PCM)
# Create a treatment indicator variable
Mortar = c(
	rep("OCM",length(OCM)), 
 	rep("PIM",length(PIM)), 
 	rep("RM",length(RM)), 
 	rep("PCM",length(PCM))
	)
# Inform R that Mortar is a factor variable
Mortar = factor(Mortar)
# Analysis of variance
fit = lm(Strength ~ Mortar)
anova(fit)


# Figure 9.4
# Page 162
f = seq(0,18,0.01)
plot(f,df(f,3,32),type="l",lty=1,xlab="F",ylab="F(3,32) pdf",cex.lab=1.25)
abline(h=0)
points(x=16.848,y=0,pch=19,cex=1.5)


# Figure 9.5
# Page 163
boxplot(OCM,PIM,RM,PCM,xlab="",names=c("OCM","PIM","RM","PCM"),ylab="Strength (in psi)",
	ylim=c(30,70),col="lightblue")
points(x=rep(0.375,length(Strength)),y=Strength,pch=19,cex=1,col="black")


# Figure 9.6
# Page 167
f = seq(0,6,0.01)
plot(f,df(f,3,32),type="l",lty=1,xlab="F",ylab="F(3,32) pdf",cex.lab=1.25)
abline(h=0)
# Add shaded area
x = seq(qf(0.95,3,32),6,0.01)
y = df(x,3,32)
polygon(c(qf(0.95,3,32),x,6),c(0,y,0),col="lightblue")
points(x=qf(0.95,3,32),y=0,pch=19,cex=1)
# Add text
text(3.75,0.06,0.05,cex=1)


# Figure 9.7
# Page 172
# This is the same as Figure 9.1.


# Example 9.1
# Data analysis
# Page 171-173
# Tukey pairwise confidence intervals
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)
# Concatenate all the data
Strength = c(OCM,PIM,RM,PCM)
# Create a treatment indicator variable
Mortar = c(
	rep("OCM",length(OCM)), 
 	rep("PIM",length(PIM)), 
 	rep("RM",length(RM)), 
 	rep("PCM",length(PCM))
	)
# Inform R that Mortar is a factor variable
Mortar = factor(Mortar)
# Analysis of variance
fit = lm(Strength ~ Mortar)
aov(fit)
TukeyHSD(aov(fit),conf.level=0.95)