Private Sub Command1_Click()
cols = Int(Val(Text1.Text))
rows = Int(Val(text2.Text))
label9.Caption = " "
label10.Caption = " "
label11.Caption = " "
label12.Caption = " "
allpvals = false: Rem change to true to specify that all P values of
possible tables are to be listed
detailed = false: Rem change to true to specify that a detailed table-by-table
report is wanted
ReDim n(rows), m(cols), ob(rows, cols), ex(rows, cols), A(rows, cols)
Open "fisher.prn" For Output As 2
If detailed Or allpvals Then Open "detail.prn" For Output As 3
Print #2, rows; "rows and "; cols; "columns"
If detailed Then Print #3, "Observed frequencies:"
Print #2, "Observed Frequencies:"
Rem read the data into array ob, i.e. observed counts
g = 0: Rem Grand total of cell contents
For i = 0 To rows - 1: grid1.Row = i
For j = 0 To cols - 1: grid1.Col = j
x = Val(grid1.Text)
ob(i + 1, j + 1) = x: A(i + 1, j + 1) = x: n(i + 1) = n(i + 1) + x:
m(j + 1) = m(j + 1) + x: g = g + x
Print #2, Format$(x, " ####.##");: If detailed Then Print #3,
x;
Next j: Print #2,: If detailed Then Print #3,
Next i
rechecking = False
recheck:
'Check for any empty rows
condensed = False
i = 1
nextrow:
If n(i) = 0 Then
If Not (condensed Or rechecking) Then MsgBox "Condensing table to remove
empty rows and/or columns"
If i < rows Then
For k = i To rows - 1: n(k) = n(k + 1)
For j = 1 To cols: ob(k, j) = ob(k + 1, j): A(k,
j) = A(k + 1, j)
Next j: Next k: n(k) = 0
End If
rows = rows - 1: If rows > 1 Then vscroll2.Value = rows
condensed = True
End If
i = i + 1: If i < rows + 1 Then GoTo nextrow
'Check for any empty columns
i = 1
nextcolumn:
If m(i) = 0 Then
If Not (condensed Or rechecking) Then MsgBox "Condensing table to remove
empty column(s)"
If i < cols Then
For k = i To cols - 1: m(k) = m(k + 1)
For j = 1 To rows: ob(j, k) = ob(j, k + 1): A(j,
k) = A(j, k + 1)
Next j: Next k: m(k) = 0
End If
cols = cols - 1: If cols > 1 Then vscroll1.Value = cols
condensed = True
End If
i = i + 1: If i < cols + 1 Then GoTo nextcolumn
'Update grid if necessary
If condensed Then
For i = 0 To rows - 1: grid1.Row = i
For j = 0 To cols - 1: grid1.Col = j
grid1.Text = Str$(ob(i + 1, j + 1))
Next j: Next i
End If
If condensed And (rows > 2 Or cols > 2) Then rechecking = True: GoTo recheck 'Make sure there are no more empty rows or columns
If rows < 2 Or cols < 2 Then MsgBox "There must be at least two rows and two columns in the contingency table": GoTo skipcalc
Form2.Hide
Cls
Rem Now calculate and print expected frequencies and observed chisquare
Rem and the exact probability of the observed table, for use in recognising
Rem more extreme tables.
Print #2,: Print #2, "Expected frequencies:"
If detailed Then Print #3, "Expected frequencies:"
For Row = 1 To rows: For Col = 1 To cols
ex(Row, Col) = n(Row) * m(Col) / g: Print #2, Format$(ex(Row, Col),
" ####.##");: If detailed Then Print #3, ex(Row, Col);
Pchis = Pchis + (ob(Row, Col) - ex(Row, Col)) ^ 2 / ex(Row, Col)
Rem Pchis is Pearson's chi-squared; Pchisy is same with Yates' correction
Pchisy = Pchisy + (Abs(ob(Row, Col) - ex(Row, Col)) - 0.5) ^ 2 / ex(Row,
Col)
Next Col: Print #2,: If detailed Then Print #3,
Next Row: Print #2,: If detailed Then Print #3,: Print #3,:
Print #2, "Chi squared ="; Pchis; "or, with Yates' correction,"; Pchisy
If detailed Then Print #3, "Chi squared ="; Pchis; "or, with Yates'
correction,"; Pchisy
Rem Prepare a look-up table of logarithmic factorials to avoid repetitive
calculation
ReDim LF(g)
For i = 1 To g: LF(i) = LF(i - 1) + Log(i): Next i
Rem and calculate the log, L1, of the factors in the P formula that
are the same for all possible contingency tables with the given marginal
totals.
L1 = 0
For i = 1 To rows: L1 = L1 + LF(n(i)): Next i
For i = 1 To cols: L1 = L1 + LF(m(i)): Next i
L1 = L1 - LF(g)
Rem Now calculate exact probability of the observed table, so that more
extreme tables can be recognised:
P = 0: Pobs = 1: Lobs = 0: Call hypergeometric(A()): Pobs = P: Lobs
= L2
Peq = 0: Rem Total Probability of all tables of probability exactly
Pobs.
Print #2, "Exact prob. of the observed table, given the marginal totals
="; P
If detailed Then Print #3,: Print #3, "Possible tables:"
Rem Estimate approximate no. of possible contingency tables that have
the given marginal totals, and hence estimate the time required to complete
the calculation
ReDim marg(rows + cols), edge(rows + cols): Rem Marg are marginal totals
plus 1;
Rem edge indicates by the no. of rows or the no. of cols which edge
of the
Rem table a particular marginal total was obtained from.
GoTo L3
subr:
For i = 1 To rows: marg(i) = n(i) + 1: edge(i) = cols - 1: Next i
For i = 1 To cols: marg(i + rows) = m(i) + 1: edge(i + rows) = rows
- 1: Next i
Return
subc:
For i = 1 To cols: marg(i) = m(i) + 1: edge(i) = rows - 1: Next i
For i = 1 To rows: marg(i + cols) = n(i) + 1: edge(i + cols) = cols
- 1: Next i
Return
L3:
If rows > cols Then GoSub subr Else GoSub subc
product = 1: factors = (rows - 1) * (cols - 1)
L4:
Min = marg(1): fact = edge(1): element = 1
For i = 2 To rows + cols: If marg(i) < Min Then Min = marg(i): fact
= edge(i): element = i
Next i
marg(element) = g: Rem to make sure this element is not used again
uses = 0
L5:
product = product * Min: uses = uses + 1: factors = factors - 1
If factors > 0 And uses < fact Then GoTo L5
product = product / uses
If factors > 0 Then GoTo L4
product = Int(product)
Print #2, "Estimated no. of possible tables with same marginal totals
="; product
text3.Text = Str$(product)
If detailed Then Print #3, "Estimated no. of possible tables with same
marginal totals ="; product
Rem Main part of program
T1 = Now: text5.Text = "calculating"
x = 1: y = 1: Rem Defines element of table at which minimisation of
L is to start
GoSub minimise
tables = 1: Rem No. of possible tables so far analysed
P = 0: PT = 0
loop1:
If detailed Then GoSub printL
Call hypergeometric(A())
If detailed Or allpvals Then Print #3, "P for this table ="; P2
If tables / 1000 = Int(tables / 1000) Then GoSub progress
GoSub increment
increased = success
x = x + 1: If x = cols Then x = 1: y = y + 1
If increased Then tables = tables + 1: GoSub minimise: GoTo loop1
'Report the results
Pans = P + Peq / 2 ' P value as defined by Anscombe
Print #2,: Print #2, "Actual no. of tables with same marginal totals
="; tables
Print #2, "Probability of a more extreme table
= Pex ="; P
Print #2, "Probability of an equally extreme table = Peq ="; Peq
Print #2, "Probability defined by Anscombe = Pex + Peq/2 ="; Pans
If detailed Then Print #3, "Probability of a more extreme table
="; P: Print #3, "Probability of an equally extreme table ="; Peq
Print #2, "Total P for all tables (should be 1)="; PT
GoSub progress
Print #2,: Print #2, "As a test of the null hypothesis that the row
categories are independent of the"
Print #2, "column categories, the author recommends that you use the
P value as defined by Anscombe,"
Print #2, "namely P ="; Pans
Print #2, "Earlier statisticians tended to use the P value giving the
probability of an equally extreme"
Print #2, "or more extreme table, namely P ="; P + Peq
If Pans > 0.000001 Then text5.Text = Format$(Pans, "0.############")
Else text5.Text = Str$(Pans)
Close
GoTo skipcalc
'Subroutine to update the display of the number of tables analysed,
and time taken so far
progress:
text4.Text = Str$(tables): T2 = Now - T1
label9.Caption = Str$(Int(T2)) + " days": T2 = 24 * (T2 - Int(T2))
label10.Caption = Str$(Int(T2)) + " hours": T2 = 60 * (T2 - Int(T2))
label11.Caption = Str$(Int(T2)) + " minutes": T2 = 60 * (T2 - Int(T2))
label12.Caption = Str$(Int(T2)) + " seconds"
label12.Refresh
DoEvents 'To allow user to shrink Fisher to an icon or run other programs
Return
Rem Subroutine to increment by 1 the lowest-order element of A(,) that
can be incremented. This method depends on the following theorem:
Rem an element A(y,x) can be incremented without altering any of the
higher-order elements if and only if there exist two non-zero elements
A(y+v,x) and A(y,x=u) where u, v > 0.
increment:
y = rows - 1: x = cols - 1: Rem Define the element that we first wish
to try to increment
success = 0
start1:
For Row = rows To y + 1 Step -1
For Col = cols To x + 1 Step -1
If A(Row, x) = 0 Or A(y, Col) = 0 Then GoTo skip1
A(y, x) = A(y, x) + 1: A(Row, Col) = A(Row, Col) + 1
A(y, Col) = A(y, Col) - 1: A(Row, x) = A(Row, x) - 1
success = -1: Row = y - 1: Col = x - 1: Rem to avoid further looping
skip1:
Next Col: Next Row
If (x = 1 And y = 1) Or success Then Return
x = x - 1: If x = 0 Then x = cols: y = y - 1
GoTo start1
Rem Subroutine to minimise the value of the lower-order elements of
L, starting from row y and column x, without affecting the higher-order
elements. This method depends on the following theorem:
Rem The element A(y,x) can be reduced without altering any of the higher-order
elements if and only if there exists a non-zero element A(y+v,x+u) where
u, v > 0.
minimise:
start3:
If y = rows Or x = cols Then GoTo skip3
start2:
success = 0
For Row = rows To y + 1 Step -1
For Col = cols To x + 1 Step -1
If A(Row, Col) = 0 Then GoTo skip2
D = A(Row, Col): If A(y, x) < D Then D = A(y, x)
A(y, x) = A(y, x) - D: A(Row, Col) = A(Row, Col) - D
A(y, Col) = A(y, Col) + D: A(Row, x) = A(Row, x) + D
success = -1
skip2:
Next Col: Next Row
If success And A(y, x) > 0 Then GoTo start2
x = x + 1: If x = cols Then x = 1: y = y + 1
GoTo start3
skip3:
Return
Rem Subroutine to print to Detail.prn the elements of array A(,)
printL:
Print #3,
For i = 1 To rows: For j = 1 To cols: Print #3, A(i, j);: Next j: Print
#3,: Next i
Return
skipcalc:
Close
End Sub
Private Sub Command2_Click()
End
End Sub
Private Sub Command3_Click()
If text5.Text = "calculating" Then GoTo skipcommand3
Open "fisher.prn" For Input As 3
While Not EOF(3)
r$ = Input$(1, 3): r1$ = r1$ + r$
Wend
Form2.Show
Form2.Text1.Text = r1$
Close 3
skipcommand3:
End Sub
Private Sub Command4_Click()
Rem Clear the data
cols = Int(Val(Text1.Text))
rows = Int(Val(text2.Text))
ReDim n(rows), m(cols), ob(rows, cols), A(rows, cols)
g = 0: Rem Grand total of cell contents
For i = 0 To rows - 1: grid1.Row = i
For j = 0 To cols - 1: grid1.Col = j
grid1.Text = ""
ob(i + 1, j + 1) = 0: A(i + 1, j + 1) = 0: n(i + 1) = 0: m(j + 1) =
0
Next j: Next i
End Sub
Private Sub Form_Load()
spacing1 = Form1.Height - command1.Top
spacing2 = Form1.Height - frame2.Top
'spacing2 = form1.Height - label5.Top
'spacing3 = form1.Height - label6.Top
spacing4 = Form1.Height - text5.Top
tolerance = 0.0000000000001
'This determines how equally extreme two tables must be in order to
be counted as equal. Tolerance = 0.01 to .00000000000001 gives the same
answers.
End Sub
Private Sub Grid1_KeyUp(keycode As Integer, Shift As Integer)
If keycode > 47 And keycode < 58 Then
grid1.Text = grid1.Text + Chr$(keycode)
End If
If keycode > 95 And keycode < 106 Then
grid1.Text = grid1.Text + Chr$(keycode - 48)
End If
If keycode = 109 Or keycode = 189 Then grid1.Text = "-" & grid1.Text
If keycode = 8 And grid1.Text <> "" Then
grid1.Text = Left$(grid1.Text, Len(grid1.Text) - 1)
End If
If keycode = 46 Then grid1.Text = ""
'If grid1.Col = 1 Then x1$(grid1.Row) = grid1.Text
'If grid1.Col = 2 Then y1$(grid1.Row) = grid1.Text
End Sub
Private Sub hypergeometric(A())
Rem subroutine to calculate exact probability of this particular table
and add it to the cumulative prob. of all tables more extreme than the
observed table, i.e. more extreme as judged by its improbability.
L2 = L1
For Row = 1 To rows: For Col = 1 To cols: L2 = L2 - LF(A(Row, Col)):
Next Col: Next Row
P2 = Exp(L2)
If detailed Then
If Abs(L2 - Lobs) < tolerance And Pobs < 1
Then
Print #3, "This is an equally
extreme table "
Else
If L2 < Lobs And Pobs
< 1 Then Print #3, "This is a table more extreme than the observed one"
End If
End If
'If P2 < Pobs Then P = P + P2
'If P2 = Pobs Then Peq = Peq + P2
If Abs(L2 - Lobs) < tolerance Then
Peq = Peq + P2
Else
If L2 < Lobs Then P = P + P2
End If
PT = PT + P2
End Sub
Private Sub Text1_Change()
cols = Int(Val(Text1.Text))
If cols > 1 Then
grid1.Width = cols * 698 + 8
grid1.cols = cols
Else
cols = 2
End If
If cols > 20 Then cols = 20
If cols > 6 Then
Form1.Width = cols * 685 + 200
Else
Form1.Width = 4968
End If
If Not changing Then vscroll1.Value = cols
End Sub
Private Sub Text2_Change()
rows = Int(Val(text2.Text))
If rows > 1 Then
grid1.Height = rows * 290 + 14
grid1.rows = rows
Else
rows = 2
End If
If rows > 20 Then rows = 20
If rows > 8 Then
Form1.Height = (rows + 14) * 285
Else
Form1.Height = 6264
End If
command1.Top = Form1.Height - spacing1
command2.Top = Form1.Height - spacing1
command3.Top = Form1.Height - spacing1
command4.Top = Form1.Height - spacing1
frame2.Top = Form1.Height - spacing2
text5.Top = Form1.Height - spacing4
label7.Top = Form1.Height - spacing4
If Not changing Then vscroll2.Value = rows
End Sub
Private Sub VScroll1_Change()
changing = True
Text1.Text = Str$(vscroll1.Value)
changing = False
End Sub
Private Sub VScroll2_Change()
changing = True
text2.Text = Str$(vscroll2.Value)
changing = False
End Sub