Options ls=75 Compress = No;
 Data;
      Infile Talent;
      input id size region age gender career weight c_type 
             plan info_1 info_2 english reading creat mech
             abstr math social physics office ses;
      if ses < 50 then physics = .;
      size=size-1;
      if size = 0 then size = 1;
      if region = 9 then region = 8;
      if size = 1 then do; s1=1; s2=-1; end;
      if size = 2 then do; s1=0; s2=2; end;
      if size = 3 then do; s1=-1; s2=-1; end;
      if plan = 1 then do; p1=-1; p2=1;  p3=-1; p4=2; end;
      if plan = 2 then do; p1=0;  p2=1;  p3=1;  p4=-3; end;
      if plan = 3 then do; p1=2;  p2=0;  p3=0;  p4=2; end;
      if plan = 4 then do; p1=0;  p2=-1; p3=-1; p4=-3; end;
      if plan = 5 then do; p1=-1; p2=-1; p3=1;  p4=2; end;
 Proc glm;
      class plan size gender region;
      model physics = math mech ses plan|size plan|gender
                     size|gender gender|region;
      lsmeans plan*size / cov out=means;
 proc IML;
        use means; reset noname; read all var _num_ into X;
        b=ncol(design(X[,1])); a=ncol(design(X[,2])); p=a-1; q=b-1;
        Sigma = X[,6:5+a*b]; mu=X[,3];
        Ha=I(a)-J(a,a,1/a); Ca=Ha[,1:p];
        Hb=I(b)-J(b,b,1/b); Cb=Hb[,1:q];
        Phi = (Cb@Ca)`*Sigma*(Cb@Ca); Psi = Ca`*shape(mu,b,a)`*Cb;
        call svd(U,D,V,Psi); wp=U[,1]; psi=shape(Psi`,p*q,1);
        Start als;
        wq=inv((I(q)@wp)`*Phi*(I(q)@wp))*(I(q)@wp)`*psi;
        wp=inv((wq@I(p))`*Phi*(wq@I(p)))*(wq@I(p))`*psi;
        epsi=psi`*(wq@wp)-R; R = R + epsi;
        finish;
        epsi = 1; R = 0;
        start iterate;
        do while(epsi >=.00001); run als; end;
        finish;
        run iterate;
        print "Maximal Contrast Coeff.: Treat. A" (Ca*wp/sqrt(wp`*Ca`*Ca*wp));
        print "Maximal Contrast Coeff.: Treat. B" (Cb*wq/sqrt(wq`*Cb`*Cb*wq));
        print "Maximal F Ratio for Product Contrast" R;