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

# 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 10.1
# Page 177
time = c(1218,1105,1191,1262,1268,1207,1070,1043,1345,1110,1227,1113,1158,1160,1060,1047,
	1043,1103,983,1144,1055,1000,1048,1075)
max.O2 = c(27.33,38.10,27.08,35.06,27.45,27.46,32.82,34.92,21.23,34.66,26.49,
	31.17,31.18,28.21,36.81,38.28,38.29,32.18,41.91,32.80,33.65,38.67,45.62,41.76)
plot(time,max.O2,xlab="Time (in seconds)",ylab="Maximum oxygen uptake (in mL/kg/min)",
	xlim=c(900,1400),ylim=c(20,50),pch=16)


# Figure 10.2
# Page 180
fit = lm(max.O2 ~ time)
fit
plot(time,max.O2,xlab="Time (in seconds)",ylab="Maximum oxygen uptake (in mL/kg/min)",
	xlim=c(900,1400),ylim=c(20,50),pch=16)
# Superimpose least-squares line
abline(fit,lwd=2,col="red") 


# Figure 10.3
# Page 182
fit = lm(max.O2 ~ time)
plot(time,max.O2,xlab="Time (in seconds)",ylab="Maximum oxygen uptake (in mL/kg/min)",
	xlim=c(1200,1300),ylim=c(24,36),pch=16)
abline(fit,col="red",lwd=2)
# Add vertical distance lines
lines(c(1218,1218),c(27.33,fit$coefficients[1] + fit$coefficients[2]*1218),col="blue",lwd=2)
lines(c(1262,1262),c(35.06,fit$coefficients[1] + fit$coefficients[2]*1262),col="blue",lwd=2)
lines(c(1268,1268),c(27.45,fit$coefficients[1] + fit$coefficients[2]*1268),col="blue",lwd=2)
lines(c(1207,1207),c(27.46,fit$coefficients[1] + fit$coefficients[2]*1207),col="blue",lwd=2)
lines(c(1227,1227),c(26.49,fit$coefficients[1] + fit$coefficients[2]*1227),col="blue",lwd=2)
# Calculate residuals and fitted values
residuals(fit)
predict(fit)


# Example 10.1
# Data analysis
# Estimate the model
fit = lm(max.O2 ~ time)
# Confidence intervals for the regression parameters (Page 185)
options(digits=2)
confint(fit,level=0.95)
# Hypothesis tests for the regression parameters (Page 186)
summary(fit)


# Figure 10.4
# Page 187
t = seq(-7,7,0.01)
pdf = dt(t,22)
plot(t,pdf,type="l",lty=1,xlab="t",ylab="t(22) pdf",cex.lab=1.25)
abline(h=0)
points(x=-6.52,y=0,pch=19,cex=1.5)


# Figure 10.5
# Page 188
fit = lm(max.O2 ~ time)
plot(time,max.O2,xlab="Time (in seconds)",ylab="Maximum oxygen uptake (in mL/kg/min)",
	xlim=c(900,1400),ylim=c(20,50),pch=16)
abline(fit,lwd=2,col="red")
abline(v=1200,lty=2)


# Figure 10.6
# Page 190
fit = lm(max.O2 ~ time)
x.0 = seq(900,1400,1)
conf = predict(fit,data.frame(time=x.0),level=0.95,interval="confidence")
pred = predict(fit,data.frame(time=x.0),level=0.95,interval="prediction")
plot(time,max.O2,xlab="Time (in seconds)",ylab="Maximum oxygen uptake (in mL/kg/min)",
	xlim=c(900,1400),ylim=c(min(pred[,2]),max(pred[,3])),pch=16)
abline(fit,lwd=2,col="red")
# Add bands for confidence intervals and prediction intervals
lines(x.0,conf[,2],type="l",lty=2)
lines(x.0,conf[,3],type="l",lty=2)
lines(x.0,pred[,2],type="l",lty=4)
lines(x.0,pred[,3],type="l",lty=4)
# Add legend
legend(1200,50,lty = c(2,4), 
	c(
	expression(paste("95% confidence")),
	expression(paste("95% prediction"))
	))


# Example 10.1
# Data analysis
# Page 191
# Estimate the model
fit = lm(max.O2 ~ time)
# Confidence interval for population mean
predict(fit,data.frame(time=1200),level=0.95,interval="confidence")
# Prediction interval
predict(fit,data.frame(time=1200),level=0.95,interval="prediction")