# amend first line to reflect exact file name and path
amh<-read.table('d:/feb08/UNC/s664/classes/amherst2.txt',header=T)
attach(amh)
x<-year-mean(year)
nreg<-lm(temp~x)
print(summary(nreg))
#
# Compute predictions and standard error for confidence intervals
#
y1<-c(165,170,175)
x1<-y1-mean(year)
pr1<-predict(nreg,data.frame(x=x1),se.fit=T)
print(pr1)
#
# Amend standard error calculation to prediction intervals
#
se1<-pr1$residual.scale*sqrt(1+(pr1$se.fit/pr1$residual.scale)^2)
#
# Print lower and upper limits for confidence and prediction intervals
#
pr1$fit-qt(0.975,160)*pr1$se.fit
pr1$fit+qt(0.975,160)*pr1$se.fit
pr1$fit-qt(0.975,160)*se1
pr1$fit+qt(0.975,160)*se1


# plot of residuals
# Fig 2.6a
plot(year,nreg$resid,xlab='Year',ylab='Residual')
# Fig 2.6b but vary "breaks" for different results
hist(nreg$resid)
hist(nreg$resid,breaks=15)

# probability plot
# 
# note we first divide the residuals by the residual
# standard deviation (summary(nreg)$sigma))
# and then compute QQ plot
qqnorm(nreg$resid/summary(nreg)$sigma)
lines(c(-1000,1000),c(-1000,1000))