###############################################
## Author: Joshua M. Tebbs
## Date: 11 October 2026
## Update: 21 October 2026
## STAT 515 course notes: R Code Chapter 8
###############################################

# 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 

# Install and load "car" R package
# For normal qq plots
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 8.1
# Page 146
# Diabetes data
# data file housed on my desktop
# read.csv imports the data into R
diabetes = read.csv(file = "C:\\Users\\tebbs\\OneDrive - University of South Carolina\\Desktop\\STAT515_data\\diabetes_gym.csv",header = T)
Calories = diabetes[,36] # extracts column for Calories
Exercise = diabetes[,37] # extracts column for Exercise
Exercise = factor(Exercise) # assures Exercise is a factor variable
summary(Exercise) # provides counts
boxplot(Calories ~ Exercise,data=diabetes,ylab="Calories burned", xlab="",ylim=c(0,2000),col="lightblue")


# Figure 8.2
# Page 148
# Sugar beet yield data
Control = c(814.8,813.2,974.9,862.0,750.8,769.0,1026.0,849.4,946.3,997.9)
Organic = c(1235.3,1185.9,1117.0,1171.8,1284.7,1211.5,1288.9,1001.4,1428.4,1373.6)
I.A.1 = c(1157.5,1236.1,1074.3,1171.5,1031.3,1015.9,950.1,1108.5,1275.8,999.4)
I.A.2 = c(955.0,1039.4,1318.6,926.9,1230.1,835.3,1013.8,1128.3,1023.7,1353.5) 
I.N.1 = c(1070.0,1153.1,940.1,998.5,1264.3,1351.1,1117.5,1389.3,1037.1,1047.3)
I.N.2 = c(1077.2,1137.7,1187.4,1335.8,1262.6,1126.7,1081.6,1134.6,1272.0,1231.3)
# Create side-by-side boxplots
boxplot(Control,Organic,I.A.1,I.A.2,I.N.1,I.N.2,
	xlab="",names=c("Control","Organic","I-A.1","I-A.2","I-N.1","I-N.2"),
	ylab="Yield (kg/acre)",ylim=c(700,1500),col="lightblue")


# Figure 8.3
# Page 149
# Mortar strength data
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 (in psi)",
	ylim=c(30,72),col="lightblue")


# Figure 8.4
# Page 151
x = seq(35,95,0.01)
# Left
plot(x,dnorm(x,65,3),type="l",lty=1,xlab="",ylab=expression(f[Y](y)),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[Y](y)),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 8.5
# Page 155 
x = seq(0,6,0.01)
pdf = df(x,3,32)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",yaxt="n",xaxt="n",bty="n")
abline(h=0)
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)
axis(side = 1, at = c(1),pos=0)
text(1,0.125,expression(1-alpha),cex=1.25)
text(3.75,0.06,expression(alpha),cex=1.25)
text(1.75,0.55,expression(F(t-1,N-t)),cex=1.25)


# Figure 8.6
# Page 156
# Left: This is the same as Figure 8.3.
# Right
x = seq(0,6,0.01)
pdf = df(x,3,32)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",yaxt="n",xaxt="n",bty="n")
abline(h=0)
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)
axis(side = 1, at = c(1),pos=0)
text(1.6,0.55,expression(F(3,32)),cex=1.25)
text(qf(0.955,3,32),-0.02,2.901,cex=1)


# Example 8.3
# Page 156-157
options(digits=6)
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 8.7
# Page 157
# P-value
x = seq(0,18,0.01)
pdf = df(x,3,32)
plot(x,pdf,type="l",lty=1,xlab="",ylab="",yaxt="n",xaxt="n",bty="n")
abline(h=0)
axis(side = 1, at = c(1),pos=0)
text(2.8,0.55,expression(F(3,32)),cex=1.25)
points(x=16.85,y=0,pch=19,cex=1.25)


# Figure 8.8
# Page 158
# Left: This is the same as Figure 8.1.
# Right: This is the same as Figure 8.2.


# Example 8.1
# Page 158
# Read in data (code above)
fit = lm(Calories ~ Exercise)
anova(fit)


# Example 8.2
# Page 159
Control = c(814.8,813.2,974.9,862.0,750.8,769.0,1026.0,849.4,946.3,997.9)
Organic = c(1235.3,1185.9,1117.0,1171.8,1284.7,1211.5,1288.9,1001.4,1428.4,1373.6)
I.A.1 = c(1157.5,1236.1,1074.3,1171.5,1031.3,1015.9,950.1,1108.5,1275.8,999.4)
I.A.2 = c(955.0,1039.4,1318.6,926.9,1230.1,835.3,1013.8,1128.3,1023.7,1353.5) 
I.N.1 = c(1070.0,1153.1,940.1,998.5,1264.3,1351.1,1117.5,1389.3,1037.1,1047.3)
I.N.2 = c(1077.2,1137.7,1187.4,1335.8,1262.6,1126.7,1081.6,1134.6,1272.0,1231.3)
# Concatenate all the data
Yield = c(Control,Organic,I.A.1,I.A.2,I.N.1,I.N.2)
# Create a treatment indicator variable
Source = c(
    rep("Control",length(Control)), 
    rep("Organic",length(Organic)), 
    rep("Inorg.A1",length(I.A.1)), 
    rep("Inorg.A2",length(I.A.2)),
    rep("Inorg.N1",length(I.N.1)), 
    rep("Inorg.N2",length(I.N.2))
    )
# Inform R that Source is a factor variable
Source = factor(Source)
# Analysis of variance
fit = lm(Yield ~ Source)
anova(fit)


# Figure 8.9
# Page 160
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)
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 8.10
# Page 164
library(car)
# Use fit object from Example 8.3 analysis (above)
residuals = residuals(fit)
library(car)
qqPlot(residuals,distribution="norm",mean=mean(residuals),sd=sd(residuals),
	xlab="Normal quantiles",ylab="Residuals",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Example 8.4
# Page 166-169
breast = c(4,22,576,342,567,86,590,8,35,1644,173)
colon = c(57,32,201,1267,144,170,428,157,58,123,861,15,152)
kidney = c(176,89,147,58,659,293,3,24,1554)
lung = c(39,427,17,460,90,187,58,52,100,200,42,167,33,31,138)
stomach = c(121,12,9,18,258,43,142,36,149,182,82,64,27)

# Figure 8.11
# Page 167
# Left
boxplot(breast,colon,kidney,lung,stomach,range=1.75,xlab="",
               names=c("Breast","Colon","Kidney","Lung","Stomach"),ylab="Survival time (in days)",
               col="lightblue")
# Right
boxplot(log(breast),log(colon),log(kidney),log(lung),log(stomach),range=1.5,xlab="",
               names=c("Breast","Colon","Kidney","Lung","Stomach"),ylab="Log survival time (in days)",
               col="lightblue")

# Concatenate all the data
Survival.time = c(breast,colon,kidney,lung,stomach)
# Create a treatment indicator variable
Cancer.type = c(
	rep("breast",length(breast)), 
 	rep("colon",length(colon)), 
 	rep("kidney",length(kidney)), 
 	rep("lung",length(lung)),
	rep("stomach",length(stomach))
	)
# Inform R that Cancer.type is a factor variable
Cancer.type = factor(Cancer.type)
# Analysis of variance
fit = lm(log(Survival.time) ~ Cancer.type) # log-transformed analysis
anova(fit)

# Figure 8.12
# Page 169
# Use fit object from Example 8.4 analysis (above)
residuals = residuals(fit)
library(car)
qqPlot(residuals,distribution="norm",mean=mean(residuals),sd=sd(residuals),
	xlab="Normal quantiles",ylab="Residuals",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Example 8.2 (Tukey analysis)
# Page 173-177

# Figure 8.13
# Page 173
# This is the same as Figure 8.2.

Control = c(814.8,813.2,974.9,862.0,750.8,769.0,1026.0,849.4,946.3,997.9)
Organic = c(1235.3,1185.9,1117.0,1171.8,1284.7,1211.5,1288.9,1001.4,1428.4,1373.6)
I.A.1 = c(1157.5,1236.1,1074.3,1171.5,1031.3,1015.9,950.1,1108.5,1275.8,999.4)
I.A.2 = c(955.0,1039.4,1318.6,926.9,1230.1,835.3,1013.8,1128.3,1023.7,1353.5) 
I.N.1 = c(1070.0,1153.1,940.1,998.5,1264.3,1351.1,1117.5,1389.3,1037.1,1047.3)
I.N.2 = c(1077.2,1137.7,1187.4,1335.8,1262.6,1126.7,1081.6,1134.6,1272.0,1231.3)
# Concatenate all the data
Yield = c(Control,Organic,I.A.1,I.A.2,I.N.1,I.N.2)
# Create a treatment indicator variable
Source = c(
    rep("Control",length(Control)), 
    rep("Organic",length(Organic)), 
    rep("Inorg.A1",length(I.A.1)), 
    rep("Inorg.A2",length(I.A.2)),
    rep("Inorg.N1",length(I.N.1)), 
    rep("Inorg.N2",length(I.N.2))
    )
# Inform R that Source is a factor variable
Source = factor(Source)
# Analysis of variance
fit = lm(Yield ~ Source)
anova(fit)
TukeyHSD(aov(fit),conf.level=0.95)


# Example 8.3 (Tukey analysis)
# Page 178-179
options(digits=6)
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)
TukeyHSD(aov(fit),conf.level=0.95)