# 
# 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(htwt)
attach(htwt)

m = lm( Wt ~ Ht )

# 2.1.1:
#
postscript("../../WriteUp/Graphics/Chapter2/prob_1.eps", onefile=FALSE, horizontal=FALSE)
plot(Ht,Wt,xlab="Wt = weight in kilograms",ylab="Ht = height in centimeters")
abline(m)
dev.off()

# 2.1.2:
#
hmeans <- mean(htwt)
n      <- 10 
hcov   <- (n-1)*var(htwt)

x_bar <- hmeans[1]
y_bar <- hmeans[2]

SXX <- hcov[1,1]
SYY <- hcov[2,2]
SXY <- hcov[1,2]

beta_hat_1 <- SXY/SXX
beta_hat_0 <- y_bar - beta_hat_1 * x_bar

# 2.1.3:
#
RSS <- SYY - SXY^2 / SXX; RSS 

sigma_hat2 <- RSS/(n-2)

sigma_hat <- sqrt( sigma_hat2 )

# compute the estimated standard errors of beta:
#
se_beta_0 = sqrt( sigma_hat2 * ( 1 / n + x_bar^2 / SXX ) )
se_beta_1 = sqrt( sigma_hat2 * ( 1 / SXX ) )

# the estimated covariance between beta_0 and beta_1:
#
cov_beta_0_beta_1 = -sigma_hat2 * x_bar / SXX

# t-tests for beta_0 = 0 and beta_1 = 0
#
t_0 = beta_hat_0 / se_beta_0
t_1 = beta_hat_1 / se_beta_1

# compute the p-values for these estimates ... compare to the t(n-2) = t(8) distribution
2 * pt(-abs(t_0),n-2) # 0.58305

2 * pt(-abs(t_1),n-2) # 0.17310

# get the ANOVA table:
#
m1 <- lm( Wt ~ Ht )
anova(m1)

#
# > anova(m1)
# Analysis of Variance Table

# Response: Wt
#           Df Sum Sq Mean Sq F value Pr(>F)
# Ht         1 159.95  159.95   2.237 0.1731
# Residuals  8 572.01   71.50

# get the F-test for regression v.s. t_1^2  :
#
F <- ( ( SYY - RSS )/1 ) / sigma_hat2

t_1^2