/* Example 6.6.1.--- Analysis of Covariance  */
proc iml;
xx1={57  60  52  49  56  46  51  63  49  57,
     29  33  25  18  25  24  20  28  18  25,
     28  30  34  33  23  32  23  21  23  28,
     25  23  33  29  17  29  16  18  22  29,
     33  31  41  35  30  22  31  24  28  30};
xx2={59  54  56  59  57  52  52,
     26  17  19  26  15  21  18,
     36  19  33  31  25  24  35,
     35  20  43  32  23  19  33,
     35  28  38  29  24  24  33};
xx3={61  59  53  59  51  51  56  58  46  53,
     25  21  26  29  24  24  22  11  15  19,
     23  21  21  12  26  17  17  24  17  17,
     11  10   6  11  22   8   8  21  12  15,
      9  11  27  11  17  19   5  24  17  18};
n1=10; n2=7; n3=10; n=n1+n2+n3; t=3; q=1; p=5;
xbj=(xx1*J(n1,1,1/n1))||(xx2*J(n2,1,1/n2))||(xx3*J(n3,1,1/n3));
xb=(xx1||xx2||xx3)*J(n,1,1/n);
v1=(xx1-xbj[,1]*J(1,n1,1))||(xx2-xbj[,2]*J(1,n2,1))
    ||(xx3-xbj[,3]*J(1,n3,1));
v=v1*v1`;
w1=((xbj[,1]-xb)*sqrt(n1))||((xbj[,2]-xb)*sqrt(n2))
    ||((xbj[,3]-xb)*sqrt(n3));
w=w1*w1`;
u1=1/v[1,1]*(v[1,1]+w[1,1])*1/det(v+w)*det(v);
u2=det(v[2:5,2:5])/det(v[2:5, 2:5]+w[2:5, 2:5]);
m=4; p=2; f=21;
F0=(f-1)*(1-sqrt(u1))/(m*sqrt(u1));
F02=21*(1-sqrt(u2))/(m*sqrt(u2));
a0=probf(F0,2*m,2*(f-1));
a02=probf(F02,2*m,2*21);
Fa=finv(1-0.05,2*m,2*(f-1));
Fa2=finv(1-0.05,2*m,2*21);
print   F0 Fa;
if F0>Fa then print "F0 > Fa: reject";
         else print "F0 < Fa: accept";
print   F02 Fa2;
if F02>Fa2 then print "F02 > Fa2: reject";
         else print "F02 < Fa2: accept";
quit;
             F0           FA
      5.0416557    2.1801705
           F0 > Fa: reject
             F02         FA2
       4.9391656   2.1681167
          F02 > Fa2: reject