(* Multivariate Analysis Package - Multiple Analysis of Variance Module
   Copyright 1986 Douglas L. Anderton.  This program may be freely
   circulated so long as it is not sold for profit and any charge does
   not exceed costs of reproduction *)

Program Manova(Input,Output);
Const
  N=30;
Type
  SUBS = 1..N;
  RVEC = Array [SUBS] Of Real;
  NBYN = Array [SUBS] Of RVEC;
  IVEC = Array [SUBS] Of Integer;
  S8 = Array [SUBS] Of String[8];
Var
  dfile, ofile : Text;
  sel : IVEC;
  miss, vars, t, u, v, w, x : RVEC;
  varn : S8;
  a, b, c, d : NBYN;
  ng, i, j, k, l, nv : Integer;
  dv, ot : SUBS;
  en, kg, hg, ig, hlg, ga, fa, dt, f1, f2, f, a1, a2, df: Real;
  err: Boolean;

Procedure openfiles(Var dfile, ofile:Text; Var ot:SUBS);
Var
  dfl, ofl:String[12];
Begin
  ClrScr; Writeln(' *** MANOVA: MULTIPLE ANALYSIS OF VARIANCE ***');
  Writeln; Write('Name of the data file? ');
  Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  Write('Name of the output file? ');
  Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  ot:= 0;
  If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  If (ot=2) Then
    Begin
    Writeln(ofile,'Multivariate Analysis Package (1.6) - ',
      'Program: MANOVA, Datafile: ',dfl); Writeln(ofile);
    End;
  End; (* Of openfiles *)

Procedure wait(Var ofile:Text; ot:Integer);
Begin
  If ot=1 Then
    Begin
    Writeln; Write('- Press any key to continue -'); While Not KeyPressed Do;
    ClrScr;
    End
  Else
    Writeln(ofile);
  End; (* Of wait *)

Procedure selcvar(Var sel:IVEC; Var varn:S8; Var miss:RVEC;
                  Var nv:Integer; Var dv:SUBS);
Var
  cfile:Text;
  cfl:String[12];
  i,j,f:Integer;
  mis:Real;
  van:String[8];
Begin
  Write('Name of the codebook file (or NONE)? '); Readln(cfl);
  If (cfl<>'NONE') And (cfl<>'none') Then f:=1 Else f:=0;
  If f=1 Then Begin Assign(cfile,cfl); Reset(cfile); End;
  Writeln;
  Write('How many variables in data file? '); Readln(nv);
  Write('Number of variables to use in MANOVA? '); Readln(dv);
  For i:=1 To dv-1 Do
   Begin
    Write('Column number for dependent variable ',i,'? '); Readln(sel[i]);
    Str(sel[i]:3,varn[i]); miss[i]:=-1E37;
    End;
  Writeln('Note: Data is presumed SORTED by treatment groups variable.');
  Write('Column number for treatment groups variable? '); Readln(sel[dv]);
  Str(sel[dv]:3,varn[dv]); miss[dv]:=-1E37;
  If f=1 Then Begin
    For j:=1 to nv Do Begin
      mis:=-1E37;
      Readln(cfile,f,van,mis);
      For i:=1 to dv Do
        If f=sel[i] Then Begin
          varn[i]:=van; miss[i]:=mis;
          Writeln('Col: ',sel[i],'  Name: ',varn[i],' Missing: ',miss[i]:6);
          End;
      End;
    Close(cfile);
    wait(cfile,1);
    End;
  End; (* Of selcvar *)

Procedure getcase(Var vars:RVEC; sel:IVEC; nv:Integer; dv:SUBS;
                  Var dfile:Text);
Var
  i:Integer;
  j:SUBS;
  x:Real;
Begin
  For i:=1 To nv Do
    Begin
    Read(dfile,x);
    For j:=1 To dv Do If (sel[j]=i) Then vars[j]:=x;
    End;
  End; (* Of getcase *)

Procedure vecout(Var x:RVEC; Var ofile:Text; varn:S8; dv:SUBS);
Var
  i:SUBS;
Begin
  For i:=1 To dv-1 Do
    Begin
    Write(ofile,varn[i]:8,' : ',x[i]:10,'   ');
    If Frac(i/3)=0 Then Writeln(ofile);
    End;
  End; (* Of vecout *)

Procedure matout(Var x:NBYN; Var ofile:Text; r,c:SUBS);
Var
  i,j:SUBS;
Begin
  For i:=1 To r Do
    Begin
    Write(ofile,'row ',i:2,': ');
    For j:=1 To c Do
      Begin
      Write(ofile,x[i,j]:10,'  ');
      If Frac(i/7)=0 Then Writeln(ofile);
      End;
    If Frac(i/7)<>0 Then Writeln(ofile);
    End;
  Writeln(ofile);
  End; (* Of matout *)

{$I MATINV.LIB}

Procedure groupout(Var a, b:NBYN; Var u, w:RVEC; dv:SUBS; ng:Real;
                   Var ofile:Text; Var hlg, hg:Real);
Var
  i,j,k:SUBS;
  det:Real;
Begin
  k:=dv-1;
  For i:=1 To k Do
    For j:=1 To k Do
      Begin
      a[i,j]:=a[i,j]-u[i]*u[j]/ng;
      b[i,j]:=b[i,j]+a[i,j];
      a[i,j]:=a[i,j]/(ng-1.0);
      End;
  For i:=1 To k Do
    Begin
    u[i]:=u[i]/ng; w[i]:=Sqrt(a[i,i]);
    End;
  ClrScr;
  Writeln; Writeln(ofile,'Means for group: ',hg:6:0,' with ',ng:5:0,' cases');
  vecout(u,ofile,varn,dv);
  Writeln(ofile); Writeln(ofile,'Standard deviations for group: ',hg:6:0);
  vecout(w,ofile,varn,dv);
  Writeln(ofile); matinv(a,k,det,err); If err Then det:=0.0;
  Writeln(ofile,'Dispersion Determinant: ',det:14);
  If (det<=0) Then det:=1.0E-20;
  hlg:=hlg+((ng-1.0)*ln(det));
  End; (* Of groupout *)

Begin (* main *)
  openfiles(dfile, ofile, ot);
  selcvar(sel, varn, miss, nv, dv);
  (* initialize *)
  FillChar(a,6*N*N,0); FillChar(b,6*N*N,0); FillChar(c,6*N*N,0); 
  FillChar(u,6*N,0); FillChar(t,6*N,0); FillChar(w,6*N,0); 
  hlg:=0.0; ga:=0.0; fa:=0.0;
  (* compute and output group level stats *)
  Writeln;
  en:=0.0; ng:=0; kg:=0;
  While Not EOF(dfile) Do
    Begin
    getcase(vars, sel, nv, dv, dfile);
    If Frac(en/10)=0.0 Then Write('+');
    j:=0;
    For i:=1 To dv Do
      If vars[i]=miss[i] Then j:=1;
    If (j=0) Or EOF(dfile) Then
      Begin
      ig:=vars[dv];
      (* output and start new group *)
      If ((hg<>ig) Or EOF(dfile)) And (en>0.0) Then
        Begin
        kg:=kg+1; dt:=ng;
        groupout(a,b,u,w,dv,dt,ofile,hlg,hg);
        fa:=fa+(1.0/(dt-1.0));
        ga:=ga+(1.0/Sqr(dt-1.0));
        wait(ofile,ot);
        FillChar(a,6*N*N,0); FillChar(u,6*N,0);
        ng:=0;
        End;
      (* accumulate group sscps *)
      If Not EOF(dfile) Then
        Begin
        en:=en+1; ng:=ng+1;
        hg:=ig;
        For j:=1 To dv-1 Do
          Begin
          u[j]:=u[j]+vars[j]; t[j]:=t[j]+vars[j];
          For i:=j To dv-1 Do
            Begin
            a[j,i]:=a[j,i]+vars[j]*vars[i]; a[i,j]:=a[j,i];
            c[j,i]:=c[j,i]+vars[j]*vars[i]; c[i,j]:=c[j,i];
            End;
          End;
        End;
      End;
    End;
  (* compute And Output non group results *)
  l:=dv-1;
  For j:=1 To l Do
    For k:=1 To l Do
      Begin
      a[j,k]:=c[j,k]-t[j]*t[k]/en; d[j,k]:=a[j,k]; c[j,k]:=b[j,k]/(en-kg);
      End;
  For j:=1 To l Do
    For k:=1 To l Do a[j,k]:=a[j,k]-b[j,k];
  For j:=1 To l Do
    Begin t[j]:=t[j]/en; u[j]:=Sqrt(c[j,j]); End;
  matinv(c,l,dt,err);
  ClrScr;
  Writeln; Writeln(ofile,'Means for Total Sample:');
  vecout(t,ofile,varn,dv);
  Writeln(ofile); Writeln(ofile,'Pooled-Samples Standard Deviations: ');
  vecout(u,ofile,varn,dv);
  wait(ofile,ot);
  Writeln(ofile); Writeln(ofile,'Between Groups SSCP Matrix: ');
  matout(a,ofile,l,l);
  wait(ofile,ot);
  Writeln(ofile); Writeln(ofile,'Within Groups SSCP Matrix: ');
  matout(b,ofile,l,l);
  wait(ofile,ot);
  Writeln(ofile); Writeln(ofile,'Inverse Pooled Dispersion Matrix: ');
  matout(c,ofile,l,l);
  Writeln(ofile); Writeln(ofile,'Dispersion Determinant: ',dt:14:4);
  wait(ofile,ot);
  dt:=(en-kg)*ln(dt)-hlg; ig:=l;
  f1:=0.5*(kg-1.0)*ig*(ig+1.0);
  fa:=(fa-(1.0/(en-kg)))*(2.0*Sqr(ig)+3.0*ig-1.0)/(6.0*(kg-1.0)*(ig+1));
  ga:=(ga-(1.0/Sqr(en-kg)))*((ig-1.0)*(ig+2.0))/(6.0*(kg-1.0));
  If (ga-Sqr(fa)) <= 0.0 Then
    Begin f2:=(f1+2.0)/(Sqr(fa)-ga);
    f:=(f2*dt)/(f1*((f2/(1.0-fa+(2.0/f2)))-dt)); End
  Else
    Begin f2:=(f1+2.0)/(ga-Sqr(fa)); f:=dt/(1.0-fa-(f1/f2)); End;
  Writeln(ofile); Writeln(ofile,'Test for Equality of Dispersions:',f:14:4);
  Writeln(ofile,'Degrees of Freedom df1:',f1:5:0,'  df2:',f2:9:0);
  wait(ofile,ot);
  Writeln(ofile,'Univariate F-tests:',(kg-1):4:0,' and',(en-kg):5:0,' d.f.');
  Writeln(ofile);
  Writeln(ofile,'variable    M. S. Between   M. S. Within         F-test');
  For j:=1 To dv-1 Do
    Begin
    Write(ofile,varn[j]:8,'  ');
    Write(ofile,(a[j,j]/(kg-1.0)):15:4);
    Write(ofile,(b[j,j]/(en-kg)):15:4);
    Writeln(ofile,((a[j,j]/(kg-1.0))/(b[j,j]/(en-kg))):15:4);
    End;
  wait(ofile,ot);
  matinv(b,dv-1,fa,err);
  matinv(d,dv-1,ga,err);
  fa:=fa/ga;
  Writeln(ofile,'Wilks Lambda:',fa:9:4,'  Eta Sq:',(1.-fa):9:4);
  If ((dv-1)<3) Or (kg<4) Then
    Begin
    f1:=2.0; f2:=en-3.0; f:=((1.0-fa)/fa)*(f2/f1);
    End
  Else
    Begin
    ga:=Sqr(ig-1.0);
    ga:=Sqrt((ga*Sqr(kg-1.0)-4.0)/(ga+Sqr(kg-1.0)-5.0));
    f1:=(ig-1.0)*(kg-1.0);
    f2:=(en-1.0)*ga-((ig-1.0)*(kg-1.0)-2.0)/2.0;
    ga:=exp(fa*ln(1.0/ga));
    f:=((1.0-ga)/ga)*(f2/f1);
    End;
  Writeln(ofile);
  Writeln(ofile,'F-test for Overall Discrimination: ',f:11:4);
  Writeln(ofile,'Degrees of Freedom df1:',f1:5:0,'  df2:',f2:5:0);
  Close(dfile); Close(ofile);
  Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
End.