/*Example 6.3.2-randomized complete block design with interaction*/
proc iml;
y1={58 4, 62 6, 60 7, 54 6, 58 6};
y2={49 7, 55 6, 48 4, 52 4, 49 7};
y3={54 6, 52 6, 51 7, 58 5, 61 6};
y4={58 4, 55 4, 51 5, 49 6, 52 4};
y=y1||y2||y3||y4;
J=4;I=5; p=2; n1=5;
m=I*J-J-I-p+1;  f=(I-1)*(J-1);   alpha=0.05;
yyb=J(1,I,1/I)*y;
ybj=yyb[1,{1 2}]`||yyb[1,{3 4}]`||yyb[1,{5 6}]`||yyb[1,{7 8}]`;
d={1 0 1 0 1 0 1 0,
   0 1 0 1 0 1 0 1};
ybi=d*Y`/J;
yb=ybj*J(J,1,1/J);
yb1=yb*J(1,I, 1);
yb2=yb*J(1,J, 1);
beta=ybi-yb1;
tao=ybj-yb2;
b1=beta*(I(5)@J(1,4,1));
ff=b1[,1:4]#tao||b1[,5:8]#tao||b1[,9:12]#tao||b1[,13:16]#tao
    ||b1[,17:20]#tao;
a1={1 0 0 0, 0 0 0 0, 0 1 0 0, 0 0 0 0, 0 0 1 0, 0 0 0 0,
    0 0 0 1, 0 0 0 0};
a2={0 0 0 0, 1 0 0 0, 0 0 0 0, 0 1 0 0, 0 0 0 0, 0 0 1 0,
    0 0 0 0, 0 0 0 1 };
z= (y[1,]*a1//y[1,]*a2)-ybj-ybi[,1]*J(1,4,1)+yb*J(1,4,1);
do k=2 to I;
z=z|| ((y[k,]*a1//y[k,]*a2)-ybj-ybi[,k]*J(1,4,1)+yb*J(1,4,1));
end;
sse=z*z`;
sh=z*ff`*inv(ff*ff`)*ff*z`;
se=sse-sh;
u=det(se)/det(sse);
F0=(m-1)*(1-sqrt(u))/(p*sqrt(u));
a0=probf(F0,2*p,2*(m-1));
Fa=finv(1-0.05,2*p,2*(m-1));
print  F0 Fa;
if F0>Fa then print "F0 > Fa: reject";
         else print "F0 < Fa: accept";
quit;
              F0          Fa
       1.3950202   2.9277442
           F0 < Fa: accept