/* Example 6.4.1.----Latin square design*/
proc iml;
a1={65 24, 66 20, 65 24, 65 26};
a2={68 21, 63 23, 67 25, 64 25};
a3={67 26, 64 25, 64 19, 64 25};
a4={68 27, 67 24, 63 20, 68 22};
y=a1//a2//a3//a4; y1=a1||a2||a3||a4;
J=4;m=J-1; f=(J-1)*(J-2);  p=2;
t=J*J(1,J**2,1/J**2)*y;
sst=y`*y-t`*t;
sa=(a1`*J(J,1,1/J))||(a2`*J(J,1,1/J))
  ||(a3`*J(J,1,1/J))||(a4`*J(J,1,1/J));
ssa=J*sa*sa`-t`*t;
sb=y1*{1 0 1 0 1 0 1 0, 0 1 0 1 0 1 0 1}`/4;
ssb=J*sb`*sb-t`*t;
c1={65 26, 67 25, 64 25, 68 27};
c2={65 24, 64 25, 64 19, 67 24};
c3={66 20, 68 21, 64 25, 63 20};
c4={65 24, 63 23, 67 26, 68 22};
sc=(c1`||c2`||c3`||c4`)*(I(J)@J(J,1,1/J));
ssc=J*sc*sc`-t`*t;
sse=sst-ssa-ssb-ssc;
u=det(sse)/det(sse+ssc);
F0=(f-1)*(1-sqrt(u))/(m*sqrt(u));
a0=probf(F0,2*m,2*(f-1));
Fa=finv(1-0.05,2*m,2*(f-1));
print F0 Fa;
if F0>Fa then print "F0 > Fa: reject";
         else print "F0 < Fa: accept";
quit;
             F0            FA
      0.9272115     3.2171754
             F0 < Fa: accept