/* Example 6.2.1.----Completely randomized 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};
n1=5;n2=5;n3=5; n4=5; n=20; p=2; J=4;
m=J-1; f=n-J;   alpha=0.05;
y=y1||y2||y3||y4;
yyb=J(1,n1,1/n1)*y;
ybj=yyb[1,{1 2}]//yyb[1,{3 4}]//yyb[1,{5 6}]//yyb[1,{7 8}];
yb=J(4,4, 1/4)*ybj;
ybb=J(n1,4,1/4)*ybj;
sstr=n1*(ybj-yb)`*(ybj-yb);
sst=(y1-ybb)`*(y1-ybb)+(y2-ybb)`*(y2-ybb)
    +(y3-ybb)`*(y3-ybb)+(y4-ybb)`*(y4-ybb);
sse=sst-sstr;
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}`; b={1 -1 0 0 }`;
m1= a`*ybj`*b;
e1=tinv(1-0.05/24, f)*sqrt(b`*b/n1*a`*sse*a/f);
low=m1-e1;
upp=m1+e1;
print  low upp;
quit;
           U0           F0           PV
    0.3838435    3.0703595     0.018312
         pv < alpha:  reject
          LOW            UPP
   (0.5391323,     15.060868)