(* Multivariate Analysis Package - N-Way Crosstabulation Module
   A copyrighted program by Douglas L. Anderton 1985.  This
   program may be freely circulated so long as it is not sold for
   profit and any charge does not exceed costs of reproduction *)

{$A-}
Program Crosstab(Input,Output);
(* Up to NSUBS variables and up to NCODES unique codes for each, so
   long as the total number of cells is less than NCELLS.  You
   should use a variable with less than 8 unique codes for the
   column variable for the most attractive printout on a 80col.
   printer. *)
Const
   NSUBS=8;
   NCODES=25;
   NCELLS=3500;
Type
   SUBS=1..NSUBS;
   CODS=1..NCODES;
   TB = Array [1..NCELLS] Of Integer;
   RC = Array [CODS] Of Real;
   IX = Array [SUBS] Of RC;
   MC = Array [CODS] Of RC;
   I8 = Array [SUBS] Of Integer;
   R8 = Array [SUBS] Of Real;
   S8 = Array [SUBS] Of String[8];
Var
   dfile, ofile : Text;
   sel, m, c : I8;
   vars : R8;
   varn : S8;
   indx : IX;
   tabl : TB;
   i, j, k, nv, dv, ot : Integer;
   nc : Real;

Procedure openfiles(Var dfile, ofile:Text; Var ot:Integer);
Var
   dfl, ofl:String[12];
Begin
   ClrScr; Writeln(' *** CROSSTAB: N-WAY TABLES AND ASSOCIATION TESTS ***');
   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: CROSSTAB, Datafile: ',dfl); Writeln(ofile);
      End;
   End; (* Of openfiles *)

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

Procedure selectvar(Var sel, m:I8; Var varn:S8;Var nv, dv:Integer);
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 CROSSTAB? '); Readln(dv);
   Write('Column number for COLUMN variable? '); Read(sel[1]);
   Write('  Max. No. Categories? '); Readln(m[1]);
   Write('Column number for ROW variable? '); Read(sel[2]);
   Write('     Max. No. Categories? '); Readln(m[2]);
   For i:=3 To dv Do
      Begin
      Write('Column number for BREAKDOWN variable ',i-2,'? ');
      Read(sel[i]); Write('  Name? '); Readln(varn[i]);
      Write('     Max. No. Categories? '); Readln(m[i]);
      End;
   If f=1 Then Begin
     For j:=1 to nv Do Begin
       Readln(cfile,f,van,mis);
       For i:=1 to dv Do
         If f=sel[i] Then Begin
           varn[i]:=van; Writeln('Col: ',sel[i],'  Name: ',varn[i]);
           End;     
       End;
     Close(cfile);
     wait(1);
     End;
   End; (* Of selectvar *)

Procedure getcase(Var vars:R8;Var sel:I8;Var nv, dv:Integer; Var dfile:Text);
Var
   i, j:Integer;
   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 tables(Var vars:R8;Var dv:Integer; Var varn:S8; Var m, c:I8;
                 Var indx:IX; Var tabl:TB);
Var
   i,j,k: Integer;
   cell: I8;
Begin
   (* for each var search c[] codes stored in indx[] for match *)
   For i:=1 To dv Do
      Begin
      cell[i]:=0;
      For j:=c[i] downto 1 Do
         If(vars[i]=indx[i,j]) Then cell[i]:=j;
      (* no match found set up new code and identify cell *)
      If cell[i]=0 Then
         Begin
         c[i]:=c[i]+1;
         If c[i]>m[i] Then
            Begin
            ClrScr; gotoXY(5,10);
            Writeln('** Error: Over',m[i]:3,' values for ',varn[i],' **');
            Delay(2000); Bdos(0);
            End;
         indx[i,c[i]]:=vars[i]; cell[i]:=c[i];
         End;
      End;
   (* add to table *)
   j:=cell[1]; k:=1;
   For i:=2 To dv Do
     Begin
     k:=k*m[i-1];
     j:=j+cell[i]*k;
     End;
   tabl[j]:=tabl[j]+1;
   End; (* Of tables *)

Procedure tabstats(Var gt:Real;Var snr, snc:Integer;Var srtot, sctot:RC;
                   Var save:MC; ot:Integer; Var ofile:Text);
Var
   i, j, k, l: Integer;
   sr, sc, ar, ac, ex, chi, temp: Real;
Begin
   (* association statistics *)
   chi:=0.; temp:=0.; k:=0; l:=0;
   For i:=1 To snr Do
     Begin
     If srtot[i]>0 Then k:=k+1;
     For j:=1 To snc Do
       Begin
       If srtot[i]>0 Then l:=l+1;
       If gt>0. Then ex:=srtot[i]*sctot[j]/gt Else ex:=0.;
       If ex<>0. Then
         Begin
         chi:=chi+Sqr(save[i,j]-ex)/ex;
         temp:=temp+Sqr(Abs(save[i,j]-ex)-0.5)/ex;
         End;
       End;
     End;
   Writeln(ofile); i:=(k-1)*(l-1);
   Writeln(ofile,'Chi-Square:',chi:10:4,' Degrees of freedom:',i:4);
   If i=1 Then Writeln(ofile,'Yale''s correction for continuity:',temp:9:4);
   If(chi+gt)<>0. Then temp:=Sqrt(chi/(chi+gt)) Else temp:=0.;
   Writeln(ofile,'Contingency Coefficient:',temp:9:4);
   temp:=k; If k<l Then temp:=l;
   If(gt*(temp-1.0))<>0. Then chi:=sqrt(chi/(gt*(temp-1.0))) Else chi:=0.;
   Writeln(ofile,'Cramer''s V:',chi:7:4);
   sc:=0.;
   For i:=1 to snc Do
     Begin
     ar:=0.;
     For j:=1 to snr Do
       If save[i,j]>ar Then ar:=save[i,j];
     sc:=sc+ar;
     End;
   ar:=0.;
   For i:=1 to snr Do
     If srtot[i]>ar Then ar:=srtot[i];
   sr:=0.;
   For i:=1 to snr Do
     Begin
     ac:=0.;
     For j:=1 to snc Do
       If save[i,j]>ac Then ac:=save[i,j];
     sr:=sr+ac;
     End;
   ac:=0.;
   For i:=1 to snc Do
     If sctot[i]>ac Then ac:=sctot[i];
   If(gt-ar)<>0. Then
     Begin
     temp:=(sc-ar)/(gt-ar);
     Writeln(ofile,'Asymmetric Lambda with Row Dependent:',temp:9:4);
     End;
   If(gt-ac)<>0. Then
     Begin
     temp:=(sr-ac)/(gt-ac);
     Writeln(ofile,'Asymmetric Lambda with Column Dependent:',temp:9:4);
     End;
   If((gt+gt)-ar-ac)<>0. Then
     Begin
     temp:=(sr+sc-ac-ar)/((gt+gt)-ar-ac);
     Writeln(ofile,'Symmetric Lambda:',temp:9:4);
     End;
   Writeln(ofile);
   wait(ot);
   End; (* Of tabstats *)

Procedure tabout(dv:Integer;Var m, c:I8; Var varn:S8;
                 Var indx:IX; Var tabl:TB; Var ofile:Text; ot:Integer);
Var
   i,j,k:Integer;
   save: MC;
   srtot, sctot: RC;
   snr, snc: Integer;
   gt: Real;
Begin
   (* save 2-way dimensions *)
   snr:=c[2]; snc:=c[1];
   If dv=2 Then c[3]:=1;
   While c[3]>0 Do
      Begin
      (* write header *)
      If dv>2 Then
         Writeln(ofile,'Breakdown ',varn[3],' = ',indx[3,c[3]]:8:2);
      Writeln(ofile,'Table rows:',varn[2],' by columns:',varn[1]);
      Writeln(ofile); Write(ofile,'          ');
      For i:=1 To snc Do
        Write(ofile,indx[1,snc-i+1]:8:2);
      FillChar(sctot,6*NCODES,0); FillChar(srtot,6*NCODES,0);
      Writeln(ofile);
      c[2]:=snr;
      While c[2]>0 Do
         Begin (* loop over rows *)
         Write(ofile,indx[2,c[2]]:8:2,'  ');
         c[1]:=snc;
         While c[1]>0 Do
            Begin (* loop over cols *)
            j:=c[1]; k:=1;
            For i:=2 To dv Do
              Begin
              k:=k*m[i-1];
              j :=j+c[i]*k;
              End;
            Write(ofile,tabl[j]:8);
            (* save 2-way table *)
            save[c[2],c[1]]:=tabl[j];
            srtot[c[2]]:=srtot[c[2]]+tabl[j];
            sctot[c[1]]:=sctot[c[1]]+tabl[j];
            c[1]:=c[1]-1;
            End;
         Writeln(ofile,srtot[c[2]]:8:0);
         c[2]:=c[2]-1;
         End;
      (* write col totals and grand *)
      Write(ofile,'          ');
      gt:=0;
      For i:=1 To snc Do
         Begin
         Write(ofile,sctot[snc-i+1]:8:0);
         gt:=gt+sctot[i];
         End;
      Writeln(ofile,gt:8:0); Writeln(ofile);
      wait(ot);
      tabstats(gt, snr, snc, srtot, sctot, save, ot, ofile);
      c[3]:=c[3]-1;
      End;
   End; (* Of tabout *)

Procedure tablop(dv,k:Integer;Var m, c:I8; Var varn:S8;
                 Var indx:IX; Var tabl:TB; Var ofile:Text; ot:integer);
Var
   l: Integer;
Begin
   k:=k-1; l:=c[k];
   If k < 4 Then l:=1;
   While (l>0) And (k>1) Do
      Begin
      (* recursive loop to final 3 levels then write *)
      If k>3 Then
         Begin
         Writeln(ofile,'Breakdown ',varn[k],' = ',indx[k,c[k]]:8:2);
         tablop(dv,k,m,c,varn,indx,tabl,ofile,ot);
         c[k]:=c[k]-1;
         End;
      If k < 4 Then tabout(dv,m,c,varn,indx,tabl,ofile,ot);
      l:=l-1;
      End;
   k:=k+1;
   End; (* Of tablop *)

Begin (* main *)
   openfiles(dfile,ofile,ot);
   (* intialize *)
   nc:=0.;
   FillChar(tabl,2*NCELLS,0); FillChar(c,2*NSUBS,0); FillChar(m,2*NSUBS,0);
   selectvar(sel,m,varn,nv,dv);
   (* accumulate tables, c=#codes in indx for var *)
   Writeln;
   k:=0;
   While Not EOF(dfile) Do
      Begin
      k:=k+1;
      getcase(vars, sel, nv, dv, dfile);
      If Frac(k/10)=0.0 Then Write('+');
      If Not EOF(dfile) Then
         Begin
         nc:=nc+1;
         tables(vars,dv,varn,m,c,indx,tabl);
         End;
      End;
   (* compute for 2-ways And Output *)
   k:=dv+1;
   ClrScr;
   tablop(dv,k,m,c,varn,indx,tabl,ofile,ot);
   Close(dfile); Close(ofile);
   Assign(dfile,'MAPSTAT.COM'); Execute(dfile);
End.