/* Example 6.3.1.-- Randomized completel block design */
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};
J=4;I=5; p=2; n1=5;
m=J-1; f=(I-1)*(J-1);   alpha=0.05;
y=y1||y2||y3||y4;
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(4,5,1/4);
ybb=J(n1,4,1/4)*ybj;
sstr=n1*(ybj-yb[,1:4]`)`*(ybj-yb[,1:4]`);
sst=(y1-ybb)`*(y1-ybb)+(y2-ybb)`*(y2-ybb)
    +(y3-ybb)`*(y3-ybb)+(y4-ybb)`*(y4-ybb);
ssbl=J*(ybi-yb)*(ybi-yb)`;
sse=sst-sstr-ssbl;
u0=det(sse)/det(sse+sstr);
F0=(f-1)*(1-sqrt(u0))/(m*sqrt(u0));
pv=1-probf(F0, 2*m, 2*(f-1));
print u0 F0 pv;
if pv<alpha then print "pv < alpha: reject";
            else print "pv > alpha: accept";
a={1 0}`;  c={1 -1 0 0 }`;
m1= a`*ybj`*c;
e1=tinv(1-0.05/24, f)*sqrt(c`*c/I*a`*sse*a/f);
low=m1-e1;
upp=m1+e1;
print     low upp ;
quit;
           U0           F0           PV
    0.3419415    2.6037374     0.0463285
         pv < alpha:  reject
              LOW         UPP
       (-0.282129,  15.882179)