###############################################
## Author: Joshua M. Tebbs
## Date: 8 September 2026
## Update: 16 September 2026
## STAT 515 course notes: R Code Chapter 6
###############################################

# 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 6.1
# Page 93
x = seq(35,75,0.01)
plot(x,dnorm(x,50,3),type="l",lty=1,xlab="y",ylab=expression(f[Y](y)),xaxt='n',ylim=c(0,0.15),cex.lab=1.25)
axis(side = 1, at = c(50,60), labels = c(expression(mu[1]),expression(mu[2])), tck = -0.01)
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 6.2
# Page 94
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.025,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.025,10)),c(0,y,0),col="lightblue")
points(x=qt(0.025,10),y=0,pch=19,cex=1)
x = seq(qt(0.975,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.975,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.975,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3.5,0.05,expression(alpha/2),cex=1.25)
text(3.5,0.05,expression(alpha/2),cex=1.25)
text(2,0.3,expression(t(n[1]+n[2]-2)),cex=1.25)


# Example 6.1
# Data analysis
# Page 96-98
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) # 5-number summary
t.test(wool,synthetic,conf.level=0.95,var.equal=TRUE)$conf.int


# Figure 6.3
# Page xx
# Boxplots
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 6.4
# Page 98
# Left
library(car)
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 6.5
# Page 99
# Top--this is the same as Figure 6.2.
# Bottom left
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-5,qt(0.05,10),0.001)
y = dt(x,10)
polygon(c(-5,x,qt(0.05,10)),c(0,y,0),col="lightblue")
points(x=qt(0.05,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(-3,0.05,expression(alpha),cex=1.25)
text(2,0.3,expression(t(n[1]+n[2]-2)),cex=1.25)
# Bottom right
x = seq(-5,5,0.001)
pdf = dt(x,10)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(qt(0.95,10),5,0.001)
y = dt(x,10)
polygon(c(qt(0.95,10),x,5),c(0,y,0),col="lightblue")
points(x=qt(0.95,10),y=0,pch=19,cex=1)
text(-0.025,0.1,expression(1-alpha),cex=1.25)
text(3,0.05,expression(alpha),cex=1.25)
text(2,0.3,expression(t(n[1]+n[2]-2)),cex=1.25)


# Example 6.2
# Data analysis
# Page 100-102
treated = c(18,43,28,50,16,32,13,35,38,33,6,7)
untreated = c(40,54,26,63,21,37,39,23,48,58,28,39)
t.test(treated,untreated,conf.level=0.95,alternative="less",var.equal=TRUE)

# Figure 6.6
# Page 100
boxplot(treated,untreated,xlab="",names=c("Treated","Untreated"),ylab="Worm count",
	ylim=c(0,70),col="lightblue")

# One-sided rejection region (lower tail) using alpha = 0.05
# Page 101
x = seq(-4,4,0.001)
pdf = dt(x,22)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-4,qt(0.05,22),0.001)
y = dt(x,22)
polygon(c(-4,x,qt(0.05,22)),c(0,y,0),col="lightblue")
points(x=qt(0.05,22),y=0,pch=19,cex=1)
text(-0.025,0.1,"0.95",cex=1.25)
text(-3,0.05,"0.05",cex=1.25)
text(1.5,0.3,"t(22)",cex=1.25)
text(qt(0.05,22),-0.011,-1.717,cex=1.1)

# p-value
# Page 102
x = seq(-4,4,0.001)
pdf = dt(x,22)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-4,-2.27,0.001)
y = dt(x,22)
polygon(c(-4,x,-2.27),c(0,y,0),col="lightblue")
points(x=-2.27,y=0,pch=19,cex=1)
text(-3.2,0.03,"0.017",cex=1.25)
text(1.5,0.3,"t(22)",cex=1.25)
text(-2.27,-0.011,-2.27,cex=1.1)


# Figure 6.7
# Page 103
x = seq(30,80,0.01)
plot(x,dnorm(x,40,2),type="l",lty=1,xlab="y",ylab=expression(f[Y](y)),xaxt='n',ylim=c(0,0.20),cex.lab=1.25)
axis(side = 1, at = c(40,55), labels = c(expression(mu[1]),expression(mu[2])), tck = -0.01)
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 6.3
# Data analysis
# Page 104-106
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 6.8
# Page 105
boxplot(transgenic,commercial,xlab="",names=c("Transgenic","Commercial"),ylab="Hatching weight (in grams)",
	ylim=c(30,60),col="lightblue")

# Figure 6.9
# Page 106
# Left
library(car)
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 6.4
# Data analysis
# Page 107-111
north = c(27.8,14.5,39.1,3.2,58.8,55.5,25,5.4,19,30.6,15.1,3.6,28.4,15,2.2,14.2,44.2,
	    25.7,11.2,46.8,36.9,54.1,10.2,2.5,13.8,43.5,13.8,39.7,6.4,4.8)
south = c(44.4,26.1,50.4,23.3,39.5,51,48.1,47.2,40.3,37.4,36.8,21.7,35.7,32,40.4,12.8,
	     5.6,44.3,52.9,38,2.6,44.6,45.5,29.1,18.7,7,43.8,28.3,36.9,51.6)
t.test(north,south,conf.level=0.90,alternative="two.sided",var.equal=FALSE)

# Figure 6.9
# Page 108
boxplot(north,south,xlab="",names=c("North","South"),ylab="Tree diameter (in cm)",
	col="lightblue")

# Rejection region
# Page 109
x = seq(-4,4,0.001)
pdf = dt(x,55.7)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-4,qt(0.05,55.7),0.001)
y = dt(x,55.7)
polygon(c(-4,x,qt(0.05,55.7)),c(0,y,0),col="lightblue")
points(x=qt(0.05,55.7),y=0,pch=19,cex=1)
x = seq(qt(0.95,55.7),4,0.001)
y = dt(x,55.7)
polygon(c(qt(0.95,55.7),x,4),c(0,y,0),col="lightblue")
points(x=qt(0.95,55.7),y=0,pch=19,cex=1)
text(-0.025,0.1,"0.90",cex=1.25)
text(-3.2,0.04,"0.05",cex=1.25)
text(3.2,0.04,"0.05",cex=1.25)
text(1.5,0.3,expression(t(nu)),cex=1.25)
text(qt(0.05,55.7),-0.011,-1.673,cex=1.1)
text(qt(0.95,55.7),-0.011,1.673,cex=1.1)

# p-value
# Page 110
x = seq(-4,4,0.001)
pdf = dt(x,55.7)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",ylim=c(0,0.4))
abline(h=0)
x = seq(-4,-2.63,0.001)
y = dt(x,55.7)
polygon(c(-4,x,-2.63),c(0,y,0),col="lightblue")
points(x=-2.63,y=0,pch=19,cex=1)
x = seq(2.63,4,0.001)
y = dt(x,55.7)
polygon(c(2.63,x,4),c(0,y,0),col="lightblue")
points(x=2.63,y=0,pch=19,cex=1)
text(1.5,0.3,expression(t(nu)),cex=1.25)
text(-2.63,-0.011,-2.63,cex=1.1)
text(2.63,-0.011,2.63,cex=1.1)


# Example 6.5
# Data analysis
# Page 112-117
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 # data differences
t.test(plasma,fat,paired=TRUE,conf.level=0.95)$conf.int
# Incorrect analysis (two independent samples)
t.test(plasma,fat,conf.level=0.95,var.equal=FALSE)$conf.int
# Analysis without outlier
t.test(diff[diff<5],conf.level=0.95)$conf.int

# Figure 6.12
# Page 117
diff = plasma-fat # data differences
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 6.13
# Page 120
x = seq(-3.5,3.5,0.001)
pdf = dnorm(x,0,1)
plot(x,pdf,type="l",lty=1,xlab="",xaxt="n",yaxt="n",bty="n",ylab="",xlim = c(-3.5,6),ylim=c(0,0.4))
abline(h=0)
x = seq(qnorm(0.95,0,1),3.5,0.001)
y = dnorm(x,0,1)
polygon(c(qnorm(0.95,0,1),x,3.5),c(0,y,0),col="lightblue")
points(x=qnorm(0.95,0,1),y=0,pch=19,cex=1)
text(2.5,0.035,expression(alpha),cex=1.25)
text(-1.25,0.3,expression(H[0]),cex=1.25)
x = seq(-2.5,5.5,0.001)
lines(x,dnorm(x,2,1),lty=4)
text(3.25,0.3,expression(H[a]),cex=1.25)
axis(side=1,at=c(0,2),labels=c(expression(0),expression(Delta)),cex.axis=1.25)