options nodate nonumber ls=76 ps=54;

proc iml;
Y={87,86,94,96,84,94,97,93,96,97,100,98};

XF={
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};


XM={
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};

XR={
1 89,
1 82,
1 88,
1 94,
1 89,
1 90,
1 91,
1 92,
1 89,
1 99,
1 84,
1 87};

I_n = i(nrow(Y));

PF = XF*inv(XF`*XF)*XF`;
PM = XM*inv(XM`*XM)*XM`;
PR = XR*inv(XR`*XR)*XR`;

SSE_F = Y`*(I_n - PF)*Y;   dfE_F = nrow(Y)-ncol(XF);
SSE_M = Y`*(I_n - PM)*Y;   dfE_M = nrow(Y)-ncol(XM);
SSE_R = Y`*(I_n - PR)*Y;   dfE_R = nrow(Y)-ncol(XR);

print "Full Model    SSE_F = " SSE_F " dfE_F = " dfE_F;
print "Middle Model    SSE_M = " SSE_M " dfE_M = " dfE_M;
print "Reduced Model    SSE_R = " SSE_R " dfE_R = " dfE_R;

/* 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 = finv(.95,dfE_M-dfE_F,dfE_F);
probF1 = 1-probf(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 = finv(.95,dfE_R-dfE_M,dfE_M);
probF2 = 1-probf(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 = finv(.95,dfE_R-dfE_F,dfE_F);
probF3 = 1-probf(F3,dfE_R-dfE_F,dfE_F);

print "Full vs Middle   " F1  F1_05  probF1;
print "Middle vs Reduced   " F2  F2_05  probF2;
print "Full vs Reduced   " F3  F3_05  probF3;


/* Bartlett's Test for Equal Variances   */

e_F = (I_n - PF)*Y;

SSE1 = e_F(|1:4,1|)`*e_F(|1:4,1|); dfE1 = 4-2;    s2_1=SSE1/dfE1;
SSE2 = e_F(|5:8,1|)`*e_F(|5:8,1|); dfE2 = 4-2;    s2_2=SSE2/dfE2;
SSE3 = 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 = cinv(.95,t-1);
probB = 1-probchi(B,t-1);

print "Baartlett's Test B = " B "X2_05 = " X2B_05 "probB = " probB;
run;

quit;