###############################################
## Author: Joshua M. Tebbs
## Date: 11 December 2011
## Update: 31 December 2017
## STAT 509 course notes: R Code
## Chapter 10
###############################################

# Example 10.1
# Sewage sludge data
# Page 154-155
filtration.rate=c(125.3,98.2,201.4,147.3,145.9,124.7,112.2,120.2,161.2,178.9,159.5,145.8,75.1,151.4,144.2,125.0,198.8,132.5,159.6,110.7)
moisture=c(77.9,76.8,81.5,79.8,78.2,78.3,77.5,77.0,80.1,80.2,79.9,79.0,76.7,78.2,79.5,78.1,81.5,77.0,79.0,78.6)
# Construct scatterplot
plot(filtration.rate,moisture,xlab="Filtration rate (kg-DS/m/hr)",ylab="Moisture (Percentage)",pch=16)

# Page 156-157
# Fit the model
fit = lm(moisture~filtration.rate)
fit
# Page 157
# Construct scatterplot (with least squares line superimposed)
plot(filtration.rate,moisture,xlab="Filtration rate (kg-DS/m/hr)",ylab="Moisture (Percentage)",pch=16)
abline(fit,lwd=2,col="red")


# Page 160-162
# Calculate fitted values and residuals
fitted.values = predict(fit) 
residuals = moisture-fitted.values
# Show residuals sum to 0
sum(residuals)
# Calculate MSres
MSres = sum(residuals^2)/(length(moisture)-2)
MSres
# Calculate residual standard error
resid.std.error = sqrt(MSres)
resid.std.error

# Page 163
# Confidence interval for beta0 and beta1
fit = lm(moisture~filtration.rate)
confint(fit,level=0.95)

# Page 164-165
# Regression output for inference
fit = lm(moisture~filtration.rate)
summary(fit)

# Page 165
# Plot t(18) pdf
t = seq(-9,9,0.01)
pdf = dt(t,18)
plot(t,pdf,type="l",lty=1,xlab="t",ylab="PDF")
abline(h=0)
points(x=8.484,y=0,pch=4,cex=1.5)

# Page 168
# Calculate confidence/prediction intervals
fit = lm(moisture~filtration.rate)
predict(fit,data.frame(filtration.rate=150),level=0.95,interval="confidence")
predict(fit,data.frame(filtration.rate=150),level=0.95,interval="prediction")

# Page 169
# Confidence/Prediction bands
fit = lm(moisture~filtration.rate)
x.not = seq(75,202,1)
conf = predict(fit,data.frame(filtration.rate=x.not),level=0.95,interval="confidence")
pred = predict(fit,data.frame(filtration.rate=x.not),level=0.95,interval="prediction")
plot(filtration.rate,moisture,xlab="Filtration rate (kg-DS/m/hr)",ylab="Moisture (Percentage)",pch=16,
	xlim=c(75,202),
	ylim=c(min(pred[,2]),max(pred[,3])))
abline(fit,lwd=2,col="red")
lines(x.not,conf[,2],type="l",lty=2)
lines(x.not,conf[,3],type="l",lty=2)
lines(x.not,pred[,2],type="l",lty=4)
lines(x.not,pred[,3],type="l",lty=4)
# Add legend
legend(80,82,lty = c(2,4), 
	c(
	expression(paste("95% confidence")),
	expression(paste("95% prediction"))
	))