pdf("muscle1.pdf")


XY1 <- matrix(c(43.7,19,177,    43.7,43,279,    43.7,56,346,    54.6,13,160,    54.6,19,193,    54.6,43,280,    54.6,56,335,    
55.7,13,169,    55.7,26,212,    55.7,34.5,244,  55.7,43,285,    58.8,13,181,    58.8,43,298,    60.5,19,212,    60.5,43,317,    
60.5,56,347,    61.9,13,186,    61.9,19,216,    61.9,34.5,265,  61.9,43,306,    61.9,56,348,    66.7,13,209,    66.7,43,324,    
66.7,56,352),   byrow=T,ncol=3)

# XY1 is a matrix containing n=24 rows and 3 columns (X1,X2,Y) we must form the X and Y matrices from it
# This is similar to when we input a raw data frame (dataset) and create a matrix

X0 <- matrix(rep(1,nrow(XY1)),ncol=1) #  Create a Column of 1s for the intercept (same number of rows as XY1)
X <- cbind(X0,XY1[,1],XY1[,2])   # Bind the columns for Intercept (X0), X1 (XY[,1]), X2 (XY[,2]) into X
Y <- as.matrix(XY1[,3])                            # Create Y from the third column of XY1

P <- X %*% solve(t(X) %*% X) %*% t(X)               # Compute P = X((X'X)^(-1))X
I_N <- diag(nrow(X))                  # Compute I(n)
J_N <- matrix(rep(1/nrow(X),nrow(X)^2),ncol=nrow(X))   # Compute (1/n)J

SS_TOT_U <- t(Y) %*% I_N %*% Y        # Compute SS(Total Uncorrected)
DF_TOT_U <- round(sum(diag(I_N)),0)         # Compute df(Total Uncorrected)

SS_MODEL  <- t(Y) %*% P %*% Y         # Compute SS(Model)
DF_MODEL <- round(sum(diag(P)),0)         # Compute df(Model)

SS_RESIDUAL  <- t(Y) %*% (I_N-P) %*% Y         # Compute SS(Residual)
DF_RESIDUAL <- round(sum(diag(I_N-P)),0)         # Compute df(Residual)

SS_MEAN  <- t(Y) %*% J_N %*% Y         # Compute SS(Model)
DF_MEAN <- round(sum(diag(J_N)),0)         # Compute df(Model)

SS_REGRESSION <- t(Y) %*% (P-J_N) %*% Y         # Compute SS(Regression)
DF_REGRESSION <- round(sum(diag(P-J_N)),0)         # Compute df(Regression)

print(cbind(SS_TOT_U,DF_TOT_U))
print(cbind(SS_MODEL,DF_MODEL))
print(cbind(SS_RESIDUAL,DF_RESIDUAL))
print(cbind(SS_MEAN,DF_MEAN))
print(cbind(SS_REGRESSION,DF_REGRESSION))

# F-test for H0: beta1=beta2=0  (Test Statistic, Critical Value (alpha=0.05,P-value)  
F_OBS <- (SS_REGRESSION/DF_REGRESSION)/(SS_RESIDUAL/DF_RESIDUAL)  
F_CRIT <- qf(.95,DF_REGRESSION,DF_RESIDUAL)
F_P <- 1-pf(F_OBS,DF_REGRESSION,DF_RESIDUAL)

print(cbind(F_OBS,F_CRIT,F_P))

dev.off()