#########################################################
# Hamster Example: source code file hamster-limits.txt  #
# Compute confidence and prediction limits.             #
#########################################################

# Read and print data frame.
hamster = read.table("hamster.txt", header=T)
cat("hamster data frame\n")
print(hamster)

# Obtain and show regression information.
reg = lm(Lifetime ~ Hiber, data=hamster)
cat("Regression summary:\n")
print(summary(reg))

# Obtain and print confidence 
# and prediction limits.
indep = data.frame(Hiber=hamster$Hiber)
pred = predict(reg, indep, interval="predict")
conf = predict(reg, indep, interval="confidence")
cat("\nConfidence limits\n")
print(conf)
cat("\nPrediction limits:\n")
print(pred)

# Plot everything.
pdf("conf-pred.pdf")
x = hamster$Hiber
y = hamster$Lifetime
plot(x, y, main="Hamster Dataset",
	       xlab="Percent of Life Hibernated",
           ylab="Lifetime in Days")
lines(x, pred[,"fit"], lty="solid", col="red")
lines(x, pred[,"lwr"], lty="dashed", col="green")
lines(x, pred[,"upr"], lty="dashed", col="green")
lines(x, conf[,"lwr"], lty="dotted", col="blue")
lines(x, conf[,"upr"], lty="dotted", col="blue")

dev.off( )