#
# 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.
#
#
# EPage 505-508 
# 
#-----


# Ex 13.1:
#

# Part (a): 
xs = c(5,10,15,20,25)
n = length(xs)
sigma = 10 
xbar = mean(xs)
S_xx = sum( ( xs - xbar )^2 )

print( sqrt( sigma^2 * ( 1 - 1/n - ( xs - xbar )^2 / S_xx ) ) )

# Part (b): 
xs = c(5,10,15,20,50)
n = length(xs)
sigma = 10 
xbar = mean(xs)
S_xx = sum( ( xs - xbar )^2 )

print( sqrt( sigma^2 * ( 1 - 1/n - ( xs - xbar )^2 / S_xx ) ) )


# Ex 13.2:
#
DF = read.csv( "../../Data/CH13/ex13-02.txt", header=TRUE, quote="'" ) 
plot( DF$x, DF$'e.', ylim=c(-2.0,+2.0) )


# Ex 13.3:
#
DF = read.csv( "../../Data/CH12/exp12-06.txt", header=TRUE, quote="'" )
colnames(DF) = c("y","x")

plot( DF$x, DF$y )

# Fit a linear model: 
m = lm( y ~ x, data=DF )
s = summary(m)

# Plot the residuals against x:
#
plot( DF$x, m$residuals )

sres = rstandard(m) # this gives the standardized residuals
plot( m$residuals/s$sigma, sres )

print( sum( m$residuals ) ) # Verify that sum_i e_i = 0
print( sum( sres ) ) # Verify that sum_i e^*_i != 0 

# Plot the standardized residuals vs. x:
#
plot( DF$x, sres )


# Ex 13.4 (EPage 506)
# 
DF = read.csv( "../../Data/CH13/ex13-04.txt", header=TRUE, quote="'" )
plot( DF$x, DF$y )

res = DF$y - ( 20.6 - 0.047 * DF$x ) # compute the residuals 
plot( DF$x, res, main='residuals vs. x' )

s = 0.7921
n = length(DF$x)
S_xx = sum( ( DF$x - mean(DF$x) )^2 )
denom = s * ( 1 - 1 / n - ( DF$x - mean(DF$x) ) / S_xx )
standardized_residuals = res / denom 
plot( DF$x, standardized_residuals )


# Ex 13.5 (EPage 506)
# 
DF = read.csv( "../../Data/CH13/ex13-05.txt", header=TRUE, quote="'" )
plot( DF$Time, DF$'Ice.Thickness' )

m = lm( Ice.Thickness ~ Time, data=DF )
sm = summary(m)

plot( DF$Time, m$residuals )


# Ex 13.6 (EPage 506)
#
DF = data.frame( c(7,9.3,13.2,16.3,19.1,22.0), c(1046,1065,1094,1117,1130,1135) ) 
colnames(DF) = c('x','y')

plot( DF$x, DF$y )

m = lm( y ~ x, data=DF )
sm = summary(m)

plot( DF$x, m$residuals )

sres = rstandard(m) # this gives the standardized residuals
plot( DF$x, sres )


# Ex 13.7 (EPage 506)
#
DF = data.frame( c(0,2,4,6,8), c(110,123,119,86,62) )
colnames(DF) = c('x','y')

plot( DF$x, DF$y )

res = DF$y - ( 127 - 6.65 * DF$x )
n = length(DF$x)
S_xx = sum( ( DF$x - mean(DF$x) )^2 )
s = 16.94
denom = s * ( 1 - 1/n - ( DF$x - mean(DF$x) ) / S_xx )
plot( DF$x, res/denom ) # the standardized residuals 


# Ex 13.8 (EPage 506/507)
#
DF = read.csv( "../../Data/CH13/ex13-08.txt", header=TRUE, quote="'" )
plot( DF$HR, DF$VO2 )

m = lm( VO2 ~ HR, data=DF )
sm = summary(m)

# Lets look for outliers using standard plots for linear regression:
#
plot( DF$HR, m$residuals )
plot( DF$HR, rstandard(m) )
plot( DF$VO2, m$fitted.values ) 
plot( DF$HR, DF$VO2 ); abline( m ) 

# Ex 13.9 (EPage 507)
#
DF = read.csv( "../../Data/CH13/ex13-09.txt", header=TRUE, quote="'" )
colnames( DF ) = c("X13","Y1","Y2","Y3","X4","Y4")

# Lets put all plots on one "figure":
#
#postscript("../../WriteUp/Graphics/Chapter13/chap_13_ex_9_anscombe_plots.eps", onefile=FALSE, horizontal=FALSE)

old_par = par(mfrow=c(4,2))

plot( DF$X13, DF$Y1 )
m = lm( Y1 ~ X13, data=DF )
#plot( DF$X13, rstandard(m) ) 
plot( m$fitted.values, rstandard(m) )

plot( DF$X13, DF$Y2 )
m = lm( Y2 ~ X13, data=DF )
#plot( DF$X13, rstandard(m) )
plot( m$fitted.values, rstandard(m) )

plot( DF$X13, DF$Y3 )
m = lm( Y3 ~ X13, data=DF )
#plot( DF$X13, rstandard(m) )
plot( m$fitted.values, rstandard(m) )

plot( DF$X4, DF$Y4 )
m = lm( Y4 ~ X4, data=DF )
#plot( DF$X4, rstandard(m) )

# Compute the standardized residuals for this case myself:
#
n = length(DF$X4)
S_xx = sum( ( DF$X4 - mean(DF$X4) )^2 )
s = summary(m)$sigma
denom = s * ( 1 - 1/n - ( DF$X4 - mean(DF$X4) )/S_xx ) 
sres = ( DF$Y4 - m$fitted.values )/denom 
plot( m$fitted.values, sres ) 

par(old_par)

#dev.off()


# Ex 13.12 (EPage 507)
#
res = c(23,-27,5,17,-8,9,15)
sum(res)

res = c(23,-27,5,17,-8,-12,2)
xs = c(3,-4,8,12,-14,-20,25)
sum( xs * res )


# Ex 13.13 (EPage 507)
#
2 * pt( -2.5, 25-2 )


# Ex 13.14 (EPage 507)
#
source('H0_linear_model.R')
DF = read.csv( "../../Data/CH13/ex13-14.txt", header=TRUE, quote="'" )
H0_linear_model( DF$x, DF$y )