airperm <- read.csv("http://www.stat.ufl.edu/~winner/sta4210/mydata/airperm_woven_reg.csv",
       header=TRUE)

attach(airperm)
names(airperm)

Y <- mean_ap
X1 <- warp
X2 <- weft
X3 <- mass 
n <- length(Y)
X0 <- rep(1,n)
X <- as.matrix(cbind(X0,X1,X2,X3))    # Form the X-matrix (n=30 rows, 4 Cols)
Y <- as.matrix(Y,ncol=1)       # Form the Y-vector (n=30 rows, 1 col)
p <- ncol(X)

# Notes:   t(X) = transpose of X, %*% = matrix multiplication, solve(A) = A^(-1) 

(XX <- t(X) %*% X)             # Obtain X'X matrix (4 rows, 4 cols)

(XY <- t(X) %*% Y)             # Obtain X'Y vector (4 rows, 1 col)

(XXI <- solve(XX))             # Obtain (X'X)^(-1) matrix (4 rows, 4 cols)

(b <- XXI %*% XY)              # Obtain b-vector (4 rows, 1 col)

Y_hat <- X %*% b               # Obtain the vector of fitted values (n=30 rows, 1 col)

e <- Y - Y_hat                 # Obtain the vector of residuals (n=30 rows, 1 col)

print(cbind(Y_hat,e))

H <- X %*% XXI %*% t(X)              # Obtain the Hat matrix

J_n <- matrix(rep(1/n,n^2),ncol=n)   # Obtain the (1/n)J matrix (n=30 rows, n=30 cols)   

I_n <- diag(n)                       # Obtain the identity matrix (n=30 rows, n=30 cols)  

(SSTO <- t(Y) %*% (I_n - J_n) %*% Y) # Obtain Total Sum of Squares (SSTO) 
# SSTO can also be computed as:
# (SSTO <- (t(Y) %*% Y) - (t(Y) %*% (I_n - J_n) %*% Y))

(SSE <- t(Y) %*% (I_n - H) %*% Y)    # Obtain Error Sum of Squares (SSE)             
# SSE can also be computed as:
# (SSE <- (t(Y) %*% Y) - (t(b) %*% XY))

(SSR <- t(Y) %*% (H - J_n) %*% Y)    # Obtain Regression Sum of Squares (SSR)             
# SSR can also be computed as:
# (SSR <- (t(b) %*% XY) - (t(Y) %*% J_n %*% Y))

(MSE <- SSE/(n-p))                   # Obtain MSE = s^2

(s2_b <- MSE[1,1] * XXI)             # Obtain s^2{b}, must use MSE[1,1] and * to do scalar multiplication

se_b <- sqrt(diag(s2_b))            # Obtain SE's of individual regression coefficients

print(cbind((b-qt(.975,n-p)*se_b),(b+qt(.975,n-p)*se_b)))  # Print CI's for Beta coefficients

(X_h <- matrix(c(1,60,35,125),ncol=1))      # Create X_h vector, for case where X1=60,X2=35,X3=125

(Y_hat_h <- t(X_h) %*% b)            # Obtain the fitted value when budget=20

(s2_yhat_h <- t(X_h) %*% s2_b %*% X_h)        # Obtain s^2{Y_hat_h}

(s2_pred <- MSE + (t(X_h) %*% s2_b %*% X_h))  # Obtain s^2{pred}

### Print 95% CI for Mean at X_h and 95% PI for Individual Observation
print(cbind((Y_hat_h-qt(.975,n-p)*sqrt(s2_yhat_h)),(Y_hat_h+qt(.975,n-p)*sqrt(s2_yhat_h))))
print(cbind((Y_hat_h-qt(.975,n-p)*sqrt(s2_pred)),(Y_hat_h+qt(.975,n-p)*sqrt(s2_pred))))