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 adprodbcprod 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]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)