Program pfish(input,output);
{}
{This program is coded in Pascal which is standard (ANSI/ISO) except}
{for the use of type  longreal.}
{}
{This program has been made available as part of the facilities of the}
{list-server facility  exact-stats  operated at the internet node}
{ mailbase@mailbase.ac.uk .}
{}
{This program undertakes the Fisher Exact Test upon 2x2 tables of counts.}
{It undertakes to provide a high standard of completeness and}
{appropriateness of reporting, using a basic 24-line scrolling text}
{terminal.  It will cope with data-counts which are very large.  It}
{provides information to support several alternative methods for}
{2-tailed testing.}
{}
{The Pascal source of this program may be distributed and used freely,}
{subject to 5 simple conditions : }
{}
{1      The source code, including this initial commentary, must not be}
{       altered in any way (except see Point 2 here below.}
{}
{2      If the code needs to be modified to comply with the requirements}
{       of a particular compiler, then an additional comment should be}
{       added at the very beginning (line 2) to explain this.}
{}
{3      The source and authorship should be acknowledged in publications}
{       wherever such acknowledgement is appropriate.}
{}
{4      Copies of the program must not be sold for substantial profit}
{       (although a reasonable handling charge is allowable).}
{}
{5      Copyright remains with the author named below.}
{}
{Fisher Exact Probability Test}
{Version 26-Jul-94}
{Norman Marsh, University of Liverpool, Liverpool L69 3BX, UK, e-mail:}
{jw34@liverpool.ac.uk .  The author would welcome comments and advice}
{regarding use of this program.}
  label 99;
  const version = 'N. Marsh, Univ. of Liverpool, 26-Jul-94';
  var n,r0,r1,r2,c0,c1,c2,ccadj,b0,f0,adprod,bcprod : integer;
      t0,tmem,lfn,m,p,p0,p1,p2,q,q0,q1,orh,orl,orthis : longreal;
      clab : char;
      fpused,opprpt,orfinite : Boolean;
      cellcount : array [1..2,1..2] of integer;
      celllabels : array [1..2,1..2] of char;
      rowtotal,coltotal : array [1..2] of integer;
  function kbyesno : Boolean;
  {Get answer to yes/no question, from keyboard}
    label 9;
    var ch : char;
  begin
    write(' ? ');
    if eoln then readln;
  9:
    read(ch);  readln;
    if (ch='y') or (ch='Y') then
      kbyesno:=true
    else
    begin
      if (ch='n') or (ch='N') then
        kbyesno:=false
      else
      begin
        write(chr(7));
        write('Please respond with ''Y(es)'' or ''N(o)''');
        goto 9
      end
    end
  end;
  function lnfactorial(f0 : integer) : longreal;
  {Evaluate ln(f0!)}
    var i : integer;
        g : longreal;
  begin
    if f0>1 then
    begin
      g:=0.0;
      for i:=2 to f0 do g:=g+ln(i);
      lnfactorial:=g
    end
    else lnfactorial:=0.0
  end;
  procedure pwrite9(pvalue : longreal);
  {Print a probability value in 9 columns of text characters}
  begin
    if ((pvalue<0.0) or (pvalue>=1.0)) then
    begin
      writeln;
      write(chr(7));
      writeln('**ILLEGAL P-VALUE ENCOUNTERED - NOTIFY PROGRAM AUTHOR!**')
    end
    else
    begin
      if pvalue<0.000010 then
      begin
        write(pvalue:9);
        fpused:=true
      end
      else
        write(pvalue:9:6)
      end
  end;
  procedure acreport;
  begin
    write('Odds Ratio (AD/BC)  = ');
    if ((cellcount[1,2]=0) or (cellcount[2,1]=0)) then
      write(' INFINITE')
    else
      write(1.0*cellcount[1,1]*cellcount[2,2]/
                (cellcount[1,2]*cellcount[2,1])
              :9:3);
    write('  ;   ');
    write('p.(conf+tail) = ');  pwrite9(q);  writeln
  end;
  function oracalc(default : longreal) : longreal;
  {calculate upper odds ratio for alternate tail, infinity -> default}
  begin
    if ((cellcount[r1,c2]=0) or (cellcount[r2,c1]=0)) then
      oracalc:=default
    else
      oracalc:=1.0*cellcount[r1,c1]*cellcount[r2,c2]/
                   (cellcount[r1,c2]*cellcount[r2,c1])
  end;

{Beginning of main program action}
begin
  {Setup Single-Character Cell Labels}
  clab:='A';
  for r0:=1 to 2 do
    for c0:=1 to 2 do
    begin
      celllabels[r0,c0]:=clab;
      clab:=succ(clab)
    end;

  {Opening Announcements, and Question about alternate tail exploration}
  write('PFISH');  write('  [');  write(version);  write(']');
  writeln;  writeln;
  writeln(
'This program accepts a set of 4 observed counts arising from a 2 x 2'
           );
  writeln(
'cross-classification, and computes the Fisher Exact Probability Test.'
         );
  writeln(
'This is a full (exact) randomisation test of departures of the given'
         );
  writeln(
'counts from the values predictable from supposed independence of the two'
         );
  writeln(
'binary distributions for the crossed classifications.  The parameters of'
         );
  writeln(
'the two binary distributions are estimated from the observed frequencies.'
         );
  writeln;
  writeln(
'This program can optionally generate information to support two-tail'
         );
  writeln(
'testing.  This involves exploring the alternate tail of the randomisation'
         );
  writeln(
'distribution, which can entail significant additional computation time.'
         );
  writeln(
'A method for 2-tail testing without this extra computation is to test the'
         );
  writeln(
'single tail with alpha set at half of the required 2-tail alpha value.'
         );
  writeln;
  writeln(
'Do you require exploration of the alternate tail of the randomisation'
         );
  write('distribution');
  opprpt:=kbyesno;
  writeln;

  {Get data from user}
  repeat
    fpused:=false;
    writeln;
    writeln(
'Type in the observed counts, each followed by a press of the  RETURN  key'
           );

    for b0:=1 to 2 do begin rowtotal[b0]:=0; coltotal[b0]:=0 end;

    for r0:=1 to 2 do
    begin
      for c0:=1 to 2 do
      repeat
        write(celllabels[r0,c0]);
        write(' (');
        write(' (Row '); write(r0:1); write(', Col '); write(c0:1);
        write('), Freq = ');
        read(f0);
        if f0<0 then
        begin
          write(chr(7));
          writeln('FAIL - Negative Frequency Specified');
          goto 99

        end;

        cellcount[r0,c0]:=f0
      until f0>=0

    end;

    {Computation for given tail}
    writeln;
    writeln('Computing - wait! ...');

    n:=0; t0:=0.0;

    for r0:=1 to 2 do
      for c0:=1 to 2 do
      begin
        f0:=cellcount[r0,c0];
        rowtotal[r0]:=rowtotal[r0]+f0;
        coltotal[c0]:=coltotal[c0]+f0;
        n:=n+f0;
        t0:=t0-lnfactorial(f0)
      end;

    tmem:=-lnfactorial(n);

    for b0:=1 to 2 do
    begin
      f0:=rowtotal[b0];
      if f0=0 then
      begin
        write(chr(7)); writeln('FAIL - Empty Row');
        goto 99

      end;

      tmem:=tmem+lnfactorial(f0);

      f0:=coltotal[b0];
      if f0=0 then
      begin
        write(chr(7)); writeln('FAIL - Empty Column');
        goto 99

      end;

      tmem:=tmem+lnfactorial(f0)
    end;

    t0:=t0+tmem;

    write('Row Proportions : ');
    for r0:=1 to 2 do write(rowtotal[r0]/n:8:5);
    writeln;

    write('Col Proportions : ');
    for c0:=1 to 2 do write(coltotal[c0]/n:8:5);
    writeln;

    adprod:=cellcount[1,1]*cellcount[2,2];
    bcprod:=cellcount[1,2]*cellcount[2,1];

    if adprod=bcprod then
    begin
      writeln(
'Data conform exactly to independence of row and col- effects'
             );
      writeln;
      goto 99
    end;

    if adprod<bcprod then
    begin
      r1:=1;  c1:=1;  r2:=2;  c2:=2
    end
    else
    begin
      if adprod>bcprod then
      begin
        r1:=1;  c1:=2;  r2:=2;  c2:=1
      end
    end;

    {Report for given tail}
    write('Low frequencies are in cells ');
    write(celllabels[r1,c1]);
    write(' and ');
    write(celllabels[r2,c2]);
    writeln;

    orfinite:=(adprod*bcprod)>0;

    if orfinite then
    begin
      if adprod>bcprod then
        orh:=1.0*adprod/bcprod
      else
        orh:=1.0*bcprod/adprod;

      orl:=1.0/orh
    end;

    write('Odds-Ratio ');
    if orfinite then
    begin
      write('AD/BC = ');
      write(1.0*adprod/bcprod:9:3);
      write(' ;   BC/AD = ');
      writeln(1.0*bcprod/adprod:9:3)
    end
    else
      writeln(' is zero/infinite');

    p:=exp(t0);
    p0:=p;

    while ((cellcount[r1,c1]>0) and (cellcount[r2,c2]>0)) do
    begin
      t0:=t0+ln(cellcount[r1,c1])+ln(cellcount[r2,c2]);
      cellcount[r1,c1]:=cellcount[r1,c1]-1;
      cellcount[r2,c2]:=cellcount[r2,c2]-1;
      cellcount[r1,c2]:=cellcount[r1,c2]+1;
      cellcount[r2,c1]:=cellcount[r2,c1]+1;
      t0:=t0-ln(cellcount[r1,c2])-ln(cellcount[r2,c1]);
      p:=p+exp(t0)
    end;

    p1:=p-p0;

    write('Prop. of Randomisations for this exact arrangement = ');
    pwrite9(p0);  writeln;

    write('Prop. of Randomisations more extreme in this direction = ');
    pwrite9(p1);  writeln;

    write(
'Prop. of Randomisations at least this extreme in this direction = '
         );
    pwrite9(p);  writeln;

    if not opprpt then goto 99;

    if p>=0.25 then
    begin
      writeln('(No account offered for the alternate tail)');
      goto 99

    end;

    {Exploration of the alternate tail}
    writeln;
    write('Alternate tail : low counts in cells ');
    write(celllabels[r1,c2]);
    write(' and ');  write(celllabels[r2,c1]:1);
    writeln;
    writeln(
'Contiguous sub-sequence of configurations for the alternate tail :'
           );
    writeln('Computing - wait! ...');

    {Generate extreme configuration of alternate tail}
    ccadj:=cellcount[r1,c2];
    if cellcount[r2,c1]<ccadj then
      ccadj:=cellcount[r2,c1];
    cellcount[r1,c1]:=cellcount[r1,c1]+ccadj;
    cellcount[r2,c2]:=cellcount[r2,c2]+ccadj;
    cellcount[r1,c2]:=cellcount[r1,c2]-ccadj;
    cellcount[r2,c1]:=cellcount[r2,c1]-ccadj;

    t0:=tmem;
    for r0:=1 to 2 do
      for c0:=1 to 2 do
        t0:=t0-lnfactorial(cellcount[r0,c0]);

    q0:=exp(t0);  q1:=0.0;  q:=q0;

    {Explore a sufficient proportion of the alternate tail}
    while ((cellcount[r1,c1]>0) and
           (cellcount[r2,c2]>0) and
           ((q<=p) or (q0<=p0) or (q1<=p1))) do
    begin
      t0:=t0+ln(cellcount[r1,c1])+ln(cellcount[r2,c2]);
      cellcount[r1,c1]:=cellcount[r1,c1]-1;
      cellcount[r2,c2]:=cellcount[r2,c2]-1;
      cellcount[r1,c2]:=cellcount[r1,c2]+1;
      cellcount[r2,c1]:=cellcount[r2,c1]+1;
      t0:=t0-ln(cellcount[r1,c2])-ln(cellcount[r2,c1]);
      q0:=exp(t0);
      q1:=q;
      q:=q+q0
    end;

    {Explore alternate tail further to include inverse odds-ratio}
    if orfinite then
      while oracalc(1.0)>orh do
      begin
        t0:=t0+ln(cellcount[r1,c1])+ln(cellcount[r2,c2]);
        cellcount[r1,c1]:=cellcount[r1,c1]-1;
        cellcount[r2,c2]:=cellcount[r2,c2]-1;
        cellcount[r1,c2]:=cellcount[r1,c2]+1;
        cellcount[r2,c1]:=cellcount[r2,c1]+1;
        t0:=t0-ln(cellcount[r1,c2])-ln(cellcount[r2,c1]);
        q0:=exp(t0);
        q1:=q;
        q:=q+q0
      end;

    {Generate and Report sub-tails of the alternate tail}
    while ((cellcount[r1,c2]>0) and
           (cellcount[r2,c1]>0) and
           (q>p1)) do
    begin
      acreport;
      t0:=t0+ln(cellcount[r1,c2])+ln(cellcount[r2,c1]);
      cellcount[r1,c2]:=cellcount[r1,c2]-1;
      cellcount[r2,c1]:=cellcount[r2,c1]-1;
      cellcount[r1,c1]:=cellcount[r1,c1]+1;
      cellcount[r2,c2]:=cellcount[r2,c2]+1;
      t0:=t0-ln(cellcount[r1,c1])-ln(cellcount[r2,c2]);
      q:=q1;  q0:=exp(t0);  q1:=q-q0
    end;

    {explore alternate tail further to include inverse odds-ratio}
    if orfinite then
      while oracalc(2*orh)<orh do
      begin
        acreport;
        t0:=t0+ln(cellcount[r1,c2])+ln(cellcount[r2,c1]);
        cellcount[r1,c2]:=cellcount[r1,c2]-1;
        cellcount[r2,c1]:=cellcount[r2,c1]-1;
        cellcount[r1,c1]:=cellcount[r1,c1]+1;
        cellcount[r2,c2]:=cellcount[r2,c2]+1;
        t0:=t0-ln(cellcount[r1,c1])-ln(cellcount[r2,c2]);
        q:=q1;  q0:=exp(t0);  q1:=q-q0
      end;

    acreport;

    write('(Number of smaller tails unreported = ');
    if cellcount[r1,c2]<cellcount[r2,c1]then
      write(cellcount[r1,c2]:1)
    else
      write(cellcount[r2,c1]:1);

    writeln(')');
    writeln;

  99:
    if fpused then
    begin
      writeln(
'NB: some probabilities have been reported in floating-point format!'
             );
      writeln
    end;

    write('Do you have further table(s) to process')

  until not kbyesno

end.