*options nodate;
*This progam uses proc iml to test the multivariate general linear
hypothesis.;  
proc iml;
/* set up a 1xG vector of cell frequencies*/
N={6 6 6};
/* set up the Gxp mean matrix*/
barX ={2.516667 3.2333333 4.15500000,
       7.666667 6.7500000 4.65000000,
       9.400000 5.7166667 5.4333333};
/* set up the pxp pooled dispersion matrix*/
Sp={3.920111111 3.205777778 3.787000000,
    3.205777778 4.258444444 2.627444444,
    3.787000000 2.627444444 5.272222222};
G=ncol(N);                                 * number of groups;
TN=0;
do i=1 to G;                               * calculate total sample size;
TN=TN+N[1,i];
end;
N=diag(N);                                 * GxG diagonal matrix N;
print N, barX, Sp, G, TN;
START MGLH;                                * subroutine for the mglh
                                             calculations;
call gsorth(R,T,lindep,A);
H=R`*barX`*C`*(inv(C*(inv(N))*C`))*C*barX*R;
df=TN-G;
E=(df)*R`*Sp*R;

u=nrow(C);                                 * calculate s, m, and n (called
                                             nn);
q=ncol(A);
s=u><q;
m=.5*((abs(q-u))-1);
nn=.5*(df-q-1);
print df s m nn;
print C, A;
LH=Trace(H*(inv(E)));                      * Lawley-Hotelling trace;
df1LH=s*(2*m+s+1);
df2LH=2*(s*nn+1);
FLH=(df2Lh/(s*df1LH))*LH;
PLH=1-probf(FLH,df1LH,df2LH);
if q>1 then do;
print 'The following are Lawley-Hotelling Results';
print LH  FLH  df1LH  df2LH  PLH;
end;
PB=Trace(H*(inv(H+E)));                    * Pillai-Bartlett trace;
df1PB=s*(2*m+s+1);
df2PB=s*(2*nn+s+1);
FPB=(df2PB/(df1PB))*(PB/(S-pb));
PPB=1-probf(FPB,df1PB,df2PB);
if q>1 then do;
print 'The following are Pillai-Bartlet Results';
print PB  FPB  df1PB  df2PB  PPB;
end;
IF q=1 then do;
f=fpb;
df1=df1pb;
df2=df2pb;
pf=ppb;
print 'A has one column, so the MANOVA test on the between-subjects factor
 specializes to a univariate ANOVA';
print f df1 df2 pf;
end;
check=A[,+];                               *univariate test of within
                                            -subjects effects;
do j=1 to q;
if check >= 0 then goto stop;
if j=q then do;
print 'The following are univariate results for repeated measures designs';
print 'Ignore these results for a doubly multivariate design';
FU=(df*q)/(u*q)*(trace(H)/trace(E));
ehat=((Trace(E/df))**2)/(q*trace((E/df)*(E/df)));
etilde=(TN*q*ehat-2)/(q*(df-q*ehat));
df1U=etilde*u*q;
df2U=etilde*df*q;
PU=1-probf(FU,df1U,df2U);
print  FU  etilde df1U  df2U  PU;
end;end;
stop:
FINISH;
print 'following for interaction';
C={1 0 -1,
   0 1 -1};
A={ 1  0,
    0  1,
   -1 -1};
run mglh;
print 'following for defect main effect';
C={1 1  1};
A={ 1  0,
    0  1,
   -1 -1};
run mglh;
print 'following for method main effect';
C={1 0 -1,
   0 1 -1};
A={ 1,
    1,
    1};
run mglh;
quit;