source("http://www.uvm.edu/~rsingle/stat221/data/scripts-221.R")

 dat <- bookdata("ch05q02.txt")
 head(dat)
 
 y.x1x2 <- lm(SBP~QUET+AGE, data=dat)
 anova(y.x1x2) 
  
#-1 Recode the data in terms of Y & X
    Y  <- dat$SBP
    X1 <- dat$QUET
    X2 <- dat$AGE
    X0 <- rep(1, times=length(X1))
     
#-2 Create the design matrix X
    X <- cbind(X0,X1,X2)  # bind 3 column vectors together
    head(X)
 
#-3 X' X matrix
    t(X) %*% X
 
#-4 beta_hat = (X'X)^-1 X' Y
    bhat <- inverse(t(X) %*% X) %*% t(X) %*% Y
    bhat      

#-5 Compare to estimates from summary()
    summary(y.x1x2)
 
#-6 Fitted values from the regression model
#   Yhat = H Y = X(X'X)^-1 X' Y  [where H is the "hat" matrix]
    y.x1x2$fit
  
#-7 Perspective Plot of the Regression Plane
    xx1 <- seq( min(X1), max(X1), length=20 )
    xx2 <- seq( min(X2), max(X2), length=20 )
 
    f <- function(xx1,xx2) { 55.32 + 9.75*xx1 + 1.05*xx2}
    y <- outer(xx1, xx2, f)
    persp(xx1, xx2, y, theta = 30, phi = 30, expand = 0.5, 
          col = "lightblue", xlab = "X1", ylab = "X2", zlab = "Y") 

#-8 Variance-Covariance matrix for Beta_i
    anova(y.x1x2)
    mse <- 79.49574
    VCV = mse * inverse(t(X) %*% X)      #VCV matrix
    round(VCV,2)
    summary(y.x1x2)$coeff[,2:3]
    summary(y.x1x2)$coeff[,2:3]^2