#
# This file contains R commands that reproduce the analysis of
# the undergraduate physics experiment described in Chapter 2.
#
# Anything not explained here will be explained in later chapters,
# as indicated.
#
# Read data from a named file. You will need to edit the text within " "
# below to identify the file in which you have stored the data. The optional
# argument header=T tells the software that the first line of the file
# gives the names of each of the variables whose values fill the
# remaining rows. The data can be accessed directly from the web-page:
#       http://www.lancs.ac.uk/staff/diggle/intro-stats-book/datasets/
# 
data<-read.table("../data/gravity.data",header=T)
names(data)
#
# assign components of data to named variables
#
distance<-data$d
time<-data$t
#
# draw scatterplot of untransformed gravity data and dsplay result
# on screen
#
plot(distance,time)
#
# draw scatterplot again, but now write results to file named within " "
# and customize the plot using some optional graphical arguments (Figure 2.3)
#
pdf(file="scatter1.pdf",width=8,height=6)
par(pty="s")
plot(distance,time,pch=19,xlim=c(0,100),ylim=c(0,0.6),
     cex.axis=1.5,cex.lab=1.5)
dev.off()
#
# scatterplot of transformed gravity data, written to named file (Figure 2.4)
#
pdf(file="scatter2.pdf",width=8,height=6)
par(pty="s")
y<-time
x<-sqrt(distance)
plot(x,y,pch=19,xlim=c(0,10),ylim=c(0,0.6),
     xlab="x = square-root of distance",
     ylab="y = time",
     cex.axis=1.5,cex.lab=1.5)
dev.off()
#
# approximate 95% confidence intervals for intercept and slope
# (will be explained in Chapters 6 and 7)
#
stuff<-lm(y~x)
ab<-stuff$coef
sigma<-summary(stuff)$sigma
cov.unscaled<-summary(stuff)$cov.unscaled
sd<-sigma*sqrt(diag(cov.unscaled))
ci<-cbind(ab-2*sd,ab+2*sd)
cibeta<-ci[2,]
cig<-2/(cibeta^2)[c(2,1)]
#
# approximate 95% confidence intervals for prediction of future 
# value of time at distance=200
#
xp<-sqrt(200)
yp<-ab[1]+ab[2]*xp
vp<-sigma*sigma*(1+c(1,xp)%*%cov.unscaled%*%c(1,xp))
cip<-c(yp-2*sqrt(vp),yp+2*sqrt(vp))
#
# approximate 95% confidence intervals for prediction of average 
# time over repeated runs with distance=200
#
vp.mean<-sigma*sigma*(c(1,xp)%*%cov.unscaled%*%c(1,xp))
cip.mean<-c(yp-2*sqrt(vp.mean),yp+2*sqrt(vp.mean))
#
# scatterplot of transformed gravity data and 95% prediction
# intervals for average time, written to file (Section 2.8)
#
pdf(file="scatter3.pdf",width=8,height=6)
par(pty="s")
y<-time
x<-sqrt(distance)
plot(x,y,pch=19,xlim=c(0,10),ylim=c(0,0.6),
     xlab="x = square-root of distance",
     ylab="y = time",
     cex.axis=1.5,cex.lab=1.5)
xmat<-cbind(rep(1,101),sqrt(0:100))
vpmat<-sigma*sigma*summary(stuff)$cov.unscaled
yp<-xmat%*%ab
sdp<-sqrt(diag(xmat%*%vpmat%*%t(xmat)))
low<-yp-2*sdp
high<-yp+2*sdp
lines(sqrt(0:100),low,lty=2)
lines(sqrt(0:100),high,lty=2)
dev.off()