#
# Epage 40
#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

library(alr3)
data(forbes)
attach(forbes)

Ktemp <- 255.37 + (5/9) * Temp
u1    <- 1 / Ktemp 

# 2.2.1: plot u1 vs. Lpres and the LS fit
plot(u1,Lpres)
abline(lm(Lpres~u1))

# 2.2.2: use the R command "lm" and display a summary: 
#
m_CC <- lm(Lpres~u1)
summary(m_CC)

# get the old (Forbes) linear model:
#
m_F <- lm(Lpres ~ Temp)

# plot Lpres vs. Temp for both models 
#
p_F  <- predict(m_F)
p_CC <- predict(m_CC)
postscript("../../WriteUp/Graphics/Chapter2/prob_2_F_vs_CC.eps", onefile=FALSE, horizontal=FALSE)
plot( p_F, p_CC, xlab="Lpres from Forbes model", ylab="Lpres from 2.8" )
dev.off()


# 2.2.4:
#
data(hooker)
u1    <- 1 / ( (5/9) * hooker$Temp + 255.37 );
Lpres <- 100 * log10( hooker$Pressure )                 # the actual pressure from the hooker data set

# fit a model based on this new data:
m_hooker <- lm(Lpres~u1)

# get the predictions under this model: 
y_hat <- predict(m_hooker)

# get the standard error of prediction ... for each of the samples from Hooker's data set using the model fitted on the Hooker data set:
n         <- length( m_hooker$residuals ) 
RSS       <- sum( m_hooker$residuals ^ 2 )
sigma_hat <- sqrt( RSS/(n-2) ) 
SXX       <- sum( ( u1 - mean(u1) )^2 )
sepred    <- sigma_hat * sqrt( 1 + 1/n + ( u1 - mean(u1) )^2 / SXX )

# compure the z scores:
z <- ( Lpres - y_hat ) / sepred 

# these may have a mean of zero and a variance near one: 
mean(z)
sqrt(var(z))

# 2.2.6:
#
u1_temp <- 1 / ( (5/9) * forbes$Temp + 255.37 )
forbes_hooker_predictions <- predict( m_hooker, newdata=data.frame(u1=u1_temp) )
sepred <- sigma_hat * sqrt( 1 + 1/n + ( u1_temp - mean(u1) )^2/ SXX )

z <- ( forbes$Lpres - forbes_hooker_predictions ) / sepred

# the mean does not seem to be zero while the variance/standard deviation does seem to be one: 
mean(z)
sqrt(var(z))