# RPD -- Examples 4.2-4.5 -- X = Ozone, y = Yield

X <- cbind(1, c(0.02, 0.07, 0.11, 0.15))

y <- rbind(242, 237, 231, 201)

XtX <- t(X) %*% X

Xty <- t(X) %*% y

XtXinv <- solve(XtX)

betahat <- XtXinv %*% Xty

SSTotalU <- drop(t(y) %*% y)  # Alternative:  sum(y^2)

  # drop() removes unwanted matrix attributes

SSModel <- drop(t(betahat) %*% Xty)

SSResid <- SSTotalU - SSModel

SSMu <- drop(t(y) %*% matrix(1/4,4,4) %*% y)

  # More efficient alternative:  4 * mean(y)^2

SSTotalC <- SSTotalU - SSMu

SSReg <- SSModel - SSMu

print(SSTotalU)
print(SSMu)
print(SSTotalC)
print(SSReg)
print(SSResid)


s2 <- SSResid/2

s2betahat <- XtXinv*s2

print(s2)
print(s2betahat)



yld.oz <- data.frame(ozone = X[,2], yield = y)

yld.oz.model <- lm(yield ~ ozone, data=yld.oz)

vcov(yld.oz.model)  # the variance-covariance matrix of betahat