\documentclass[notitlepage]{article}
\usepackage{Math533}
\usepackage{listings}
\usepackage{numprint}


\lstloadlanguages{R}

\definecolor{keywordcolor}{rgb}{0,0.6,0.6}
\definecolor{delimcolor}{rgb}{0.461,0.039,0.102}
\definecolor{Rcommentcolor}{rgb}{0.101,0.043,0.432}

\lstdefinestyle{Rsettings}{
  basicstyle=\ttfamily,
  breaklines=true,
  showstringspaces=false,
  keywords={if, else, function, theFunction, tmp}, % Write as many keywords
  otherkeywords={},
  commentstyle=\itshape\color{Rcommentcolor},
  keywordstyle=\color{keywordcolor},
  moredelim=[s][\color{delimcolor}]{"}{"},
}

\lstset{basicstyle=\ttfamily, numbers=none, literate={~} {$\sim$}{2}}

\begin{document}
\begin{center}
{\textsclarge{Prediction for the Rocket Propulsion Data}}
\end{center}

<<predict1,comment='+', fig.width=7.25, fig.height=4.25, fig.align='center'>>=
data.source<-"http://www.math.mcgill.ca/dstephens/Regression/Data/2-1-RocketProp.csv"
RocketProp<-read.csv(file=data.source)
names(RocketProp)<-c('i','Strength','Age')
x<-RocketProp$Age
y<-RocketProp$Strength
n<-length(x)
par(mar=c(4,4,0,1))
plot(x,y,pch=19,xlab='Age',ylab='Shear Strength')
fit.RP<-lm(y~x)
summary(fit.RP)
@

\pagebreak

Predictions on a grid of new x values using the \texttt{predict} function:
<<predict2,comment='+', fig.width=8.25, fig.height=6.25, fig.align='center'>>=
xnew<-seq(0,25,by=0.1)

#Confidence interval
ynew.interval<-predict(fit.RP,newdata=data.frame(x=xnew),interval='confidence')
plot(x,y,pch=19,xlab='Age',ylab='Shear Strength')
abline(coef(fit.RP),col='red')
title('Prediction for Rocket Propulsion Data')
lines(xnew,ynew.interval[,2],lty=2,col='red')
lines(xnew,ynew.interval[,3],lty=2,col='red')

#Prediction interval
yonew.interval<-predict(fit.RP,newdata=data.frame(x=xnew),interval='prediction')
lines(xnew,yonew.interval[,2],lty=2,col='blue')
lines(xnew,yonew.interval[,3],lty=2,col='blue')
legend(20,2600,c('Prediction','Conf. Interv.','Pred. Interv.'),col=c('red','red','blue'),lty=c(1,2,2))
@

Using the matrix formulae:

<<predict4,comment='+', fig.width=5.25, fig.height=4.25, fig.align='center'>>=
X<-cbind(rep(1,n),x)
(XtX<-t(X)%*%X)
Xty<-t(X) %*% y
beta.hat<-solve(XtX,Xty)
Hat.matrix<-X %*% solve(XtX) %*% t(X)

SS.Res<-(t(y)%*%(diag(1,n)-Hat.matrix)%*%y)[1,1]
MS.Res<-SS.Res/(n-2)
sigma.hat<-sqrt(MS.Res)

nnew<-length(xnew)
Xnew<-cbind(rep(1,nnew),xnew)

Ynew<-Xnew%*%beta.hat[,1]
Hnew<-diag(Xnew %*% solve(XtX) %*% t(Xnew))
Var.Ynew<-sigma.hat^2*Hnew
Var.Yonew<-sigma.hat^2*(1+Hnew)

#Look up the quantile
alpha<-0.05
q.val<-qt(1-alpha/2,n-2)

#Intervals
Conf.vals<-cbind(Ynew-q.val*sqrt(Var.Ynew),Ynew+q.val*sqrt(Var.Ynew))
Pred.vals<-cbind(Ynew-q.val*sqrt(Var.Yonew),Ynew+q.val*sqrt(Var.Yonew))

#Check computation - should be zero
sum(ynew.interval[,2:3]-Conf.vals)
sum(yonew.interval[,2:3]-Pred.vals)
@
\end{document}