###############################################
## Author: Joshua M. Tebbs
## Date: 31 August 2025
## Update: 9 September 2025
## STAT 509 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
# See https://cran.r-project.org/web/packages/car/car.pdf for documentation

# Figure 8.1
# Page 131
x = seq(35,75,0.01)
plot(x,dnorm(x,50,3),type="l",lty=1,xlab="x",ylab=expression(f[X](x)),xaxt='n',ylim=c(0,0.15),cex.lab=1.25)
axis(side = 1, at = c(40,50,60,70), labels = FALSE, tck = -0.015)
lines(x,dnorm(x,60,3),lty=4)
abline(h=0)
# Add legend
legend(63,0.15,lty = c(1,4), 
	c(
	expression(paste("Population 1")),
	expression(paste("Population 2"))
	))


# Figure 8.2
# Page 132
# This figure is the same as Figure 7.1.


# Example 8.1
# Data analysis
# Page 134-136
wool = c(116.6,123.8,148.8,127.7,141.3,129.1,139.9,149.4,148.8,130.1,
	135.9,125.5,128.8,154.7,155.1,155.9,132.9,165.5,150.9,139.7)
synthetic = c(157.2,144.5,157.9,161.3,144.4,165.8,169.3,163.4,151.3,140.2,
	157.1,162.6,158.7,161.9,173.3,148.3,149.0,158.8,149.0,162.2,132.9,170.3,134.1,151.4,138.6)
quantile(wool,type=2) # 5-number summary
quantile(synthetic,type=2)
t.test(wool,synthetic,conf.level=0.95,var.equal=TRUE)$conf.int


# Figure 8.3
# Page 134
wool = c(116.6,123.8,148.8,127.7,141.3,129.1,139.9,149.4,148.8,130.1,
	135.9,125.5,128.8,154.7,155.1,155.9,132.9,165.5,150.9,139.7)
synthetic = c(157.2,144.5,157.9,161.3,144.4,165.8,169.3,163.4,151.3,140.2,
	157.1,162.6,158.7,161.9,173.3,148.3,149.0,158.8,149.0,162.2,132.9,170.3,134.1,151.4,138.6)
boxplot(wool,synthetic,xlab="",names=c("Wool","Synthetic"),ylab="Breaking strength (in MPa)",
	ylim=c(100,200),col="lightblue")


# Figure 8.4
# Page 136
# Left
qqPlot(wool,distribution="norm",mean=mean(wool),sd=sd(wool),
	xlab="Normal quantiles",ylab="Breaking strength (in MPa)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)
# Right
qqPlot(synthetic,distribution="norm",mean=mean(synthetic),sd=sd(synthetic),
	xlab="Normal quantiles",ylab="Breaking strength (in MPa)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 8.5
# Page 137
x = seq(30,80,0.01)
plot(x,dnorm(x,40,2),type="l",lty=1,xlab="x",ylab=expression(f[X](x)),xaxt='n',ylim=c(0,0.20),cex.lab=1.25)
axis(side = 1, at = c(30,40,50,60,70,80), labels = FALSE, tck = -0.015)
lines(x,dnorm(x,55,6),lty=4)
abline(h=0)
# Add legend
legend(60,0.175,lty = c(1,4), 
	c(
	expression(paste("Population 1")),
	expression(paste("Population 2"))
	))


# Example 8.2
# Data analysis
# Page 138-140
transgenic = c(38.8,39.0,39.7,40.0,40.8,40.9,41.0,41.0,41.0,42.5,
	42.6,43.0,43.0,43.4,43.5,43.5,43.8,44.4,44.7,44.7,
	44.7,45.3,45.7,45.8,46.4,46.5,46.6,46.7,46.7,46.8,
	46.9,47.1,47.1,47.1,47.3,47.6,47.7,48.1,48.3,49.3,
	49.3,49.8,50.3,50.9,52.1)
commercial = c(36.7,37.1,38.9,39.5,39.5,39.8,40.0,40.2,40.3,40.5,
	40.5,40.7,41.1,41.2,41.5,41.5,41.6,41.6,41.7,42.4,
	43.1,43.3,43.3,43.4,43.7,44.1,44.2,45.2,45.3,45.4,
	46.0,46.1,46.4,46.6,46.6,46.9,47.3,47.5,48.1,48.2,
	48.4,48.6,49.0,49.1,49.3,49.6,50.1,50.2,50.4,50.6,
	52.2,53.0,55.5,56.4)
t.test(transgenic,commercial,conf.level=0.99,var.equal=FALSE)$conf.int


# Figure 8.6
# Page 139
boxplot(transgenic,commercial,xlab="",names=c("Transgenic","Commercial"),ylab="Hatching weight (in grams)",
	ylim=c(30,60),col="lightblue")


# Figure 8.7
# Page 140
# Left
qqPlot(transgenic,distribution="norm",mean=mean(transgenic),sd=sd(transgenic),
	xlab="Normal quantiles",ylab="Hatching weight (in grams)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)
# Right
qqPlot(commercial,distribution="norm",mean=mean(commercial),sd=sd(commercial),
	xlab="Normal quantiles",ylab="Hatching weight (in grams)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Example 8.3
# Data analysis
# Page 141-144
plasma = c(2.5,3.1,2.1,3.5,3.1,1.8,6.0,3.0,36.0,4.7,6.9,3.9,4.2,1.6,7.2,1.8,20.0,2.0,2.5,4.1)
fat = c(4.9,5.9,4.4,6.9,9.0,4.2,10.0,5.5,41.0,4.4,7.0,2.9,4.6,1.4,7.7,1.1,11.0,2.5,2.3,1.5)
diff = plasma-fat
# Correct analysis (matched pairs)
t.test(diff,conf.level=0.95)$conf.int
# Incorrect analysis (two independent samples)
t.test(plasma,fat,conf.level=0.95,var.equal=TRUE)$conf.int
# Analysis without outlier
t.test(diff[diff<5],conf.level=0.95)$conf.int


# Figure 8.8
# Page 144
qqPlot(diff,distribution="norm",mean=mean(diff),sd=sd(diff),
	xlab="Normal quantiles",ylab="Observed data differences",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 8.9
# Page 146
q = seq(0,6,0.01)
plot(q,df(q,5,10),type="l",lty=1,xlab="q",ylab=expression(f[Q](q)),cex.lab=1.25,ylim=c(0,1.1))
lines(q,df(q,10,20),lty=4)
lines(q,df(q,20,40),lty=8)
abline(h=0)
# Add legend
legend(3,0.8,lty = c(1,4,8),
	c(expression(paste(nu[1]," = 5, ",nu[2]," = 10")),
	expression(paste(nu[1]," = 10, ",nu[2], " = 20")),
	expression(paste(nu[1]," = 20, ", nu[2]," = 40"))
	))


# Figure 8.10
# Page 147
# This figure is the same as Figure 7.4.


# Example 8.4
# Data analysis
# Page 148-150
rotary = c(127.75,127.87,127.86,127.92,128.03,127.94,127.91,128.10,128.01,128.11,127.79,127.93,
           127.89,127.96,127.80,127.94,128.02,127.82,128.11,127.92,127.74,127.78,127.85,127.96)
inline = c(127.90,127.90,127.74,127.93,127.62,127.76,127.63,127.93,127.86,127.73,127.82,127.84,
           128.06,127.88,127.85,127.60,128.02,128.05,127.95,127.89,127.82,127.92,127.71,127.78)
var.test(inline,rotary,conf.level=0.95)$conf.int   # inline = population 2


# Figure 8.11
# Page 149
boxplot(rotary,inline,xlab="",names=c("Rotary","Inline"),ylab="Fill volume (in fluid ounces)",
	ylim=c(127.5,128.2),col="lightblue")


# Figure 8.12
# Page 150
# Left
qqPlot(rotary,distribution="norm",mean=mean(rotary),sd=sd(rotary),
	xlab="Normal quantiles",ylab="Fill volume (in fluid ounces)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)
# Right
qqPlot(inline,distribution="norm",mean=mean(inline),sd=sd(inline),
	xlab="Normal quantiles",ylab="Fill volume (in fluid ounces)",pch=16,
	envelope=list(border=TRUE,style="lines"),id=FALSE)


# Figure 8.13
# Page 152
# This figure is the same as Figure 7.1.


# Example 8.5
# Data analysis
# Page 153
prop.test(c(68,46),c(165,159),conf.level=0.95,correct=FALSE)$conf.int