home *** CD-ROM | disk | FTP | other *** search
/ Oakland CPM Archive / oakcpm.iso / cpm / turbopas / mapstatf.lbr / MANOVA.PZS / MANOVA.ÐAS
Text File  |  1987-01-24  |  9KB  |  292 lines
  1. (* Multivariate Analysis Package - Multiple Analysis of Variance Module
  2.    Copyright 1986 Douglas L. Anderton.  This program may be freely
  3.    circulated so long as it is not sold for profit and any charge does
  4.    not exceed costs of reproduction *)
  5.  
  6. Program Manova(Input,Output);
  7. Const
  8.   N=30;
  9. Type
  10.   SUBS = 1..N;
  11.   RVEC = Array [SUBS] Of Real;
  12.   NBYN = Array [SUBS] Of RVEC;
  13.   IVEC = Array [SUBS] Of Integer;
  14.   S8 = Array [SUBS] Of String[8];
  15. Var
  16.   dfile, ofile : Text;
  17.   sel : IVEC;
  18.   miss, vars, t, u, v, w, x : RVEC;
  19.   varn : S8;
  20.   a, b, c, d : NBYN;
  21.   ng, i, j, k, l, nv : Integer;
  22.   dv, ot : SUBS;
  23.   en, kg, hg, ig, hlg, ga, fa, dt, f1, f2, f, a1, a2, df: Real;
  24.   err: Boolean;
  25.  
  26. Procedure openfiles(Var dfile, ofile:Text; Var ot:SUBS);
  27. Var
  28.   dfl, ofl:String[12];
  29. Begin
  30.   ClrScr; Writeln(' *** MANOVA: MULTIPLE ANALYSIS OF VARIANCE ***');
  31.   Writeln; Write('Name of the data file? ');
  32.   Readln(dfl); Assign(dfile,dfl); Reset(dfile);
  33.   Write('Name of the output file? ');
  34.   Readln(ofl); Assign(ofile,ofl); Rewrite(ofile);
  35.   ot:= 0;
  36.   If (ofl='CON:') Or (ofl='con:') Then ot:=1;
  37.   If (ofl='LST:') Or (ofl='lst:') Then ot:=2;
  38.   If (ot=2) Then
  39.     Begin
  40.     Writeln(ofile,'Multivariate Analysis Package (1.6) - ',
  41.       'Program: MANOVA, Datafile: ',dfl); Writeln(ofile);
  42.     End;
  43.   End; (* Of openfiles *)
  44.  
  45. Procedure wait(Var ofile:Text; ot:Integer);
  46. Begin
  47.   If ot=1 Then
  48.     Begin
  49.     Writeln; Write('- Press any key to continue -'); While Not KeyPressed Do;
  50.     ClrScr;
  51.     End
  52.   Else
  53.     Writeln(ofile);
  54.   End; (* Of wait *)
  55.  
  56. Procedure selcvar(Var sel:IVEC; Var varn:S8; Var miss:RVEC;
  57.                   Var nv:Integer; Var dv:SUBS);
  58. Var
  59.   cfile:Text;
  60.   cfl:String[12];
  61.   i,j,f:Integer;
  62.   mis:Real;
  63.   van:String[8];
  64. Begin
  65.   Write('Name of the codebook file (or NONE)? '); Readln(cfl);
  66.   If (cfl<>'NONE') And (cfl<>'none') Then f:=1 Else f:=0;
  67.   If f=1 Then Begin Assign(cfile,cfl); Reset(cfile); End;
  68.   Writeln;
  69.   Write('How many variables in data file? '); Readln(nv);
  70.   Write('Number of variables to use in MANOVA? '); Readln(dv);
  71.   For i:=1 To dv-1 Do
  72.    Begin
  73.     Write('Column number for dependent variable ',i,'? '); Readln(sel[i]);
  74.     Str(sel[i]:3,varn[i]); miss[i]:=-1E37;
  75.     End;
  76.   Writeln('Note: Data is presumed SORTED by treatment groups variable.');
  77.   Write('Column number for treatment groups variable? '); Readln(sel[dv]);
  78.   Str(sel[dv]:3,varn[dv]); miss[dv]:=-1E37;
  79.   If f=1 Then Begin
  80.     For j:=1 to nv Do Begin
  81.       mis:=-1E37;
  82.       Readln(cfile,f,van,mis);
  83.       For i:=1 to dv Do
  84.         If f=sel[i] Then Begin
  85.           varn[i]:=van; miss[i]:=mis;
  86.           Writeln('Col: ',sel[i],'  Name: ',varn[i],' Missing: ',miss[i]:6);
  87.           End;
  88.       End;
  89.     Close(cfile);
  90.     wait(cfile,1);
  91.     End;
  92.   End; (* Of selcvar *)
  93.  
  94. Procedure getcase(Var vars:RVEC; sel:IVEC; nv:Integer; dv:SUBS;
  95.                   Var dfile:Text);
  96. Var
  97.   i:Integer;
  98.   j:SUBS;
  99.   x:Real;
  100. Begin
  101.   For i:=1 To nv Do
  102.     Begin
  103.     Read(dfile,x);
  104.     For j:=1 To dv Do If (sel[j]=i) Then vars[j]:=x;
  105.     End;
  106.   End; (* Of getcase *)
  107.  
  108. Procedure vecout(Var x:RVEC; Var ofile:Text; varn:S8; dv:SUBS);
  109. Var
  110.   i:SUBS;
  111. Begin
  112.   For i:=1 To dv-1 Do
  113.     Begin
  114.     Write(ofile,varn[i]:8,' : ',x[i]:10,'   ');
  115.     If Frac(i/3)=0 Then Writeln(ofile);
  116.     End;
  117.   End; (* Of vecout *)
  118.  
  119. Procedure matout(Var x:NBYN; Var ofile:Text; r,c:SUBS);
  120. Var
  121.   i,j:SUBS;
  122. Begin
  123.   For i:=1 To r Do
  124.     Begin
  125.     Write(ofile,'row ',i:2,': ');
  126.     For j:=1 To c Do
  127.       Begin
  128.       Write(ofile,x[i,j]:10,'  ');
  129.       If Frac(i/7)=0 Then Writeln(ofile);
  130.       End;
  131.     If Frac(i/7)<>0 Then Writeln(ofile);
  132.     End;
  133.   Writeln(ofile);
  134.   End; (* Of matout *)
  135.  
  136. {$I MATINV.LIB}
  137.  
  138. Procedure groupout(Var a, b:NBYN; Var u, w:RVEC; dv:SUBS; ng:Real;
  139.                    Var ofile:Text; Var hlg, hg:Real);
  140. Var
  141.   i,j,k:SUBS;
  142.   det:Real;
  143. Begin
  144.   k:=dv-1;
  145.   For i:=1 To k Do
  146.     For j:=1 To k Do
  147.       Begin
  148.       a[i,j]:=a[i,j]-u[i]*u[j]/ng;
  149.       b[i,j]:=b[i,j]+a[i,j];
  150.       a[i,j]:=a[i,j]/(ng-1.0);
  151.       End;
  152.   For i:=1 To k Do
  153.     Begin
  154.     u[i]:=u[i]/ng; w[i]:=Sqrt(a[i,i]);
  155.     End;
  156.   ClrScr;
  157.   Writeln; Writeln(ofile,'Means for group: ',hg:6:0,' with ',ng:5:0,' cases');
  158.   vecout(u,ofile,varn,dv);
  159.   Writeln(ofile); Writeln(ofile,'Standard deviations for group: ',hg:6:0);
  160.   vecout(w,ofile,varn,dv);
  161.   Writeln(ofile); matinv(a,k,det,err); If err Then det:=0.0;
  162.   Writeln(ofile,'Dispersion Determinant: ',det:14);
  163.   If (det<=0) Then det:=1.0E-20;
  164.   hlg:=hlg+((ng-1.0)*ln(det));
  165.   End; (* Of groupout *)
  166.  
  167. Begin (* main *)
  168.   openfiles(dfile, ofile, ot);
  169.   selcvar(sel, varn, miss, nv, dv);
  170.   (* initialize *)
  171.   FillChar(a,6*N*N,0); FillChar(b,6*N*N,0); FillChar(c,6*N*N,0); 
  172.   FillChar(u,6*N,0); FillChar(t,6*N,0); FillChar(w,6*N,0); 
  173.   hlg:=0.0; ga:=0.0; fa:=0.0;
  174.   (* compute and output group level stats *)
  175.   Writeln;
  176.   en:=0.0; ng:=0; kg:=0;
  177.   While Not EOF(dfile) Do
  178.     Begin
  179.     getcase(vars, sel, nv, dv, dfile);
  180.     If Frac(en/10)=0.0 Then Write('+');
  181.     j:=0;
  182.     For i:=1 To dv Do
  183.       If vars[i]=miss[i] Then j:=1;
  184.     If (j=0) Or EOF(dfile) Then
  185.       Begin
  186.       ig:=vars[dv];
  187.       (* output and start new group *)
  188.       If ((hg<>ig) Or EOF(dfile)) And (en>0.0) Then
  189.         Begin
  190.         kg:=kg+1; dt:=ng;
  191.         groupout(a,b,u,w,dv,dt,ofile,hlg,hg);
  192.         fa:=fa+(1.0/(dt-1.0));
  193.         ga:=ga+(1.0/Sqr(dt-1.0));
  194.         wait(ofile,ot);
  195.         FillChar(a,6*N*N,0); FillChar(u,6*N,0);
  196.         ng:=0;
  197.         End;
  198.       (* accumulate group sscps *)
  199.       If Not EOF(dfile) Then
  200.         Begin
  201.         en:=en+1; ng:=ng+1;
  202.         hg:=ig;
  203.         For j:=1 To dv-1 Do
  204.           Begin
  205.           u[j]:=u[j]+vars[j]; t[j]:=t[j]+vars[j];
  206.           For i:=j To dv-1 Do
  207.             Begin
  208.             a[j,i]:=a[j,i]+vars[j]*vars[i]; a[i,j]:=a[j,i];
  209.             c[j,i]:=c[j,i]+vars[j]*vars[i]; c[i,j]:=c[j,i];
  210.             End;
  211.           End;
  212.         End;
  213.       End;
  214.     End;
  215.   (* compute And Output non group results *)
  216.   l:=dv-1;
  217.   For j:=1 To l Do
  218.     For k:=1 To l Do
  219.       Begin
  220.       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);
  221.       End;
  222.   For j:=1 To l Do
  223.     For k:=1 To l Do a[j,k]:=a[j,k]-b[j,k];
  224.   For j:=1 To l Do
  225.     Begin t[j]:=t[j]/en; u[j]:=Sqrt(c[j,j]); End;
  226.   matinv(c,l,dt,err);
  227.   ClrScr;
  228.   Writeln; Writeln(ofile,'Means for Total Sample:');
  229.   vecout(t,ofile,varn,dv);
  230.   Writeln(ofile); Writeln(ofile,'Pooled-Samples Standard Deviations: ');
  231.   vecout(u,ofile,varn,dv);
  232.   wait(ofile,ot);
  233.   Writeln(ofile); Writeln(ofile,'Between Groups SSCP Matrix: ');
  234.   matout(a,ofile,l,l);
  235.   wait(ofile,ot);
  236.   Writeln(ofile); Writeln(ofile,'Within Groups SSCP Matrix: ');
  237.   matout(b,ofile,l,l);
  238.   wait(ofile,ot);
  239.   Writeln(ofile); Writeln(ofile,'Inverse Pooled Dispersion Matrix: ');
  240.   matout(c,ofile,l,l);
  241.   Writeln(ofile); Writeln(ofile,'Dispersion Determinant: ',dt:14:4);
  242.   wait(ofile,ot);
  243.   dt:=(en-kg)*ln(dt)-hlg; ig:=l;
  244.   f1:=0.5*(kg-1.0)*ig*(ig+1.0);
  245.   fa:=(fa-(1.0/(en-kg)))*(2.0*Sqr(ig)+3.0*ig-1.0)/(6.0*(kg-1.0)*(ig+1));
  246.   ga:=(ga-(1.0/Sqr(en-kg)))*((ig-1.0)*(ig+2.0))/(6.0*(kg-1.0));
  247.   If (ga-Sqr(fa)) <= 0.0 Then
  248.     Begin f2:=(f1+2.0)/(Sqr(fa)-ga);
  249.     f:=(f2*dt)/(f1*((f2/(1.0-fa+(2.0/f2)))-dt)); End
  250.   Else
  251.     Begin f2:=(f1+2.0)/(ga-Sqr(fa)); f:=dt/(1.0-fa-(f1/f2)); End;
  252.   Writeln(ofile); Writeln(ofile,'Test for Equality of Dispersions:',f:14:4);
  253.   Writeln(ofile,'Degrees of Freedom df1:',f1:5:0,'  df2:',f2:9:0);
  254.   wait(ofile,ot);
  255.   Writeln(ofile,'Univariate F-tests:',(kg-1):4:0,' and',(en-kg):5:0,' d.f.');
  256.   Writeln(ofile);
  257.   Writeln(ofile,'variable    M. S. Between   M. S. Within         F-test');
  258.   For j:=1 To dv-1 Do
  259.     Begin
  260.     Write(ofile,varn[j]:8,'  ');
  261.     Write(ofile,(a[j,j]/(kg-1.0)):15:4);
  262.     Write(ofile,(b[j,j]/(en-kg)):15:4);
  263.     Writeln(ofile,((a[j,j]/(kg-1.0))/(b[j,j]/(en-kg))):15:4);
  264.     End;
  265.   wait(ofile,ot);
  266.   matinv(b,dv-1,fa,err);
  267.   matinv(d,dv-1,ga,err);
  268.   fa:=fa/ga;
  269.   Writeln(ofile,'Wilks Lambda:',fa:9:4,'  Eta Sq:',(1.-fa):9:4);
  270.   If ((dv-1)<3) Or (kg<4) Then
  271.     Begin
  272.     f1:=2.0; f2:=en-3.0; f:=((1.0-fa)/fa)*(f2/f1);
  273.     End
  274.   Else
  275.     Begin
  276.     ga:=Sqr(ig-1.0);
  277.     ga:=Sqrt((ga*Sqr(kg-1.0)-4.0)/(ga+Sqr(kg-1.0)-5.0));
  278.     f1:=(ig-1.0)*(kg-1.0);
  279.     f2:=(en-1.0)*ga-((ig-1.0)*(kg-1.0)-2.0)/2.0;
  280.     ga:=exp(fa*ln(1.0/ga));
  281.     f:=((1.0-ga)/ga)*(f2/f1);
  282.     End;
  283.   Writeln(ofile);
  284.   Writeln(ofile,'F-test for Overall Discrimination: ',f:11:4);
  285.   Writeln(ofile,'Degrees of Freedom df1:',f1:5:0,'  df2:',f2:5:0);
  286.   Close(dfile); Close(ofile);
  287.   Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
  288. End.
  289. 'variable    M. S. Between   M. S. Within         F-test');
  290.   For j:=1 To dv-1 Do
  291.     Begin
  292.