> Y <- matrix(c(87,86,94,96,84,94,97,93,96,97,100,98),byrow=T,ncol=1)
> 
> XF <- matrix(c(
+ 1, 89, 0,  0, 0,  0,
+ 1, 82, 0,  0, 0,  0,
+ 1, 88, 0,  0, 0,  0,
+ 1, 94, 0,  0, 0,  0,
+ 0,  0, 1, 89, 0,  0,
+ 0,  0, 1, 90, 0,  0,
+ 0,  0, 1, 91, 0,  0,
+ 0,  0, 1, 92, 0,  0,
+ 0,  0, 0,  0, 1, 89,
+ 0,  0, 0,  0, 1, 99,
+ 0,  0, 0,  0, 1, 84,
+ 0,  0, 0,  0, 1, 87), byrow=T,ncol=6)
> 
> XM <- matrix(c(
+ 1, 0, 0, 89,
+ 1, 0, 0, 82,
+ 1, 0, 0, 88,
+ 1, 0, 0, 94,
+ 0, 1, 0, 89,
+ 0, 1, 0, 90,
+ 0, 1, 0, 91,
+ 0, 1, 0, 92,
+ 0, 0, 1, 89,
+ 0, 0, 1, 99,
+ 0, 0, 1, 84,
+ 0, 0, 1, 87), byrow=T, ncol=4)
> 
> XR <- matrix(c(
+ 1, 89,
+ 1, 82,
+ 1, 88,
+ 1, 94,
+ 1, 89,
+ 1, 90,
+ 1, 91,
+ 1, 92,
+ 1, 89,
+ 1, 99,
+ 1, 84,
+ 1, 87), byrow=T, ncol=2)
> 
> I_n <- diag(nrow(Y))
> 
> PF <- XF %*% solve(t(XF) %*% XF) %*% t(XF)
> PM <- XM %*% solve(t(XM) %*% XM) %*% t(XM)
> PR <- XR %*% solve(t(XR) %*% XR) %*% t(XR)
> 
> SSE_F <- t(Y) %*% (I_n - PF) %*% Y;    dfE_F <- nrow(Y) - ncol(XF)
> SSE_M <- t(Y) %*% (I_n - PM) %*% Y;    dfE_M <- nrow(Y) - ncol(XM)
> SSE_R <- t(Y) %*% (I_n - PR) %*% Y;    dfE_R <- nrow(Y) - ncol(XR)
> 
> cat("Full Model   SSE_F =", SSE_F, "dfE_F =", dfE_F, "\n")
Full Model   SSE_F = 86.08416 dfE_F = 6 
> cat("Middle Model   SSE_M =", SSE_M, "dfE_M =", dfE_M, "\n")
Middle Model   SSE_M = 164.2775 dfE_M = 8 
> cat("Reduced Model   SSE_R =", SSE_R, "dfE_R =", dfE_R, "\n") 
Reduced Model   SSE_R = 269.9488 dfE_R = 10 
> 
> #### Test for Equality of Slopes   Full Model vs Middle Model
> 
> F1 <- ((SSE_M-SSE_F)/(dfE_M-dfE_F))/(SSE_F/dfE_F)
> F1_05 <-  qf(.95,dfE_M-dfE_F,dfE_F)
> probF1  <-  1-pf(F1,dfE_M-dfE_F,dfE_F)
> 
> #### Test for Equality of Intercepts, Given Equal Slopes - Mid vs Reduced 
> 
> F2 <- ((SSE_R-SSE_M)/(dfE_R-dfE_M))/(SSE_M/dfE_M)
> F2_05 <- qf(.95,dfE_R-dfE_M,dfE_M)
> probF2 <- 1-pf(F2,dfE_R-dfE_M,dfE_M)
> 
> ### Joint Test for Equality of Slopes and Intercepts - Full vs Reduced  
> 
> F3 <- ((SSE_R-SSE_F)/(dfE_R-dfE_F))/(SSE_F/dfE_F)
> F3_05 <- qf(.95,dfE_R-dfE_F,dfE_F)
> probF3 <- 1-pf(F3,dfE_R-dfE_F,dfE_F)
> 
> 
> cat("Full vs Middle:  F1 =", F1, "F1_05 =", F1_05, "probF1 =", probF1, "\n")
Full vs Middle:  F1 = 2.725008 F1_05 = 5.143253 probF1 = 0.1438916 
> cat("Middle vs Reduced:  F2 =", F2, "F2_05 =", F2_05, "probF2 =", probF2, "\n")
Middle vs Reduced:  F2 = 2.572996 F2_05 = 4.45897 probF2 = 0.1371471 
> cat("Full vs Reduced:  F3 =", F3, "F3_05 =", F3_05, "probF3 =", probF3, "\n")
Full vs Reduced:  F3 = 3.203807 F3_05 = 4.533677 probF3 = 0.09869023 
> 
> ##### Bartlett's Test for equal variances
> 
> e_F <- (I_n - PF) %*% Y
> 
> SSE1  <-  t(e_F[1:4,1]) %*% e_F[1:4,1];   dfE1  <-  4-2;    s2_1 <- SSE1/dfE1;
> SSE2  <-  t(e_F[5:8,1]) %*% e_F[5:8,1];   dfE2  <-  4-2;    s2_2 <- SSE2/dfE2;
> SSE3  <-  t(e_F[9:12,1]) %*% e_F[9:12,1]; dfE3  <-  4-2;    s3_2 <- SSE3/dfE3;
> dfE <- dfE1+dfE2+dfE3
> MSE  <-  (SSE1 + SSE2 + SSE3)/dfE
> 
> t <- 3
> C  <-  1 + (1/(3*(t-1)))*(1/dfE1 + 1/dfE2 + 1/dfE3 - 1/dfE)
> B  <-  (1/C)*(dfE*log(MSE) - (dfE1*log(s2_1)+dfE2*log(s2_2)+dfE3*log(s3_2)))
> X2B_05  <-  qchisq(.95,t-1);
> probB  <-  1-pchisq(B,t-1);
> 
> cat("Bartlett's Test:  B=", B, "X2_05 =", X2B_05, "probB =", probB, "\n")
Bartlett's Test:  B= 1.593242 X2_05 = 5.991465 probB = 0.4508498 
> 
> 
>