Random Contingency Tables
Vector 1:3, page 160
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 160
A contingency table is a 2-dimensional array of positive numbers together with their row and column sums and the overall total, such as
| 3009 | 2832 | 3008 | 8849 |
| 3047 | 3051 | 2992 | 9095 |
| 2974 | 3038 | 3018 | 9030 |
| 9030 | 8921 | 9023 | 26974 |
An obvious example is a spread-sheet of sales figures. A standard problem in simulation is that of finding random figures which have the same row and column totals as the originals; for instance, for finding possible sales figures that would have given the same overall results. From this one can tell whether an outstandingly high, or low, figure could reasonably be due to chance. In APL terms, this involves finding, for an RxC dimensional original, an NxRxC array of integers with the right marginal totals.
We tackle this problem according to two sound APL principles:
i) Work with arrays, not loops, where possible.
ii) Use existing idioms where possible.
The existing idiom for n-dimensional arrays of independent random normal deviates is well developed, e.g.
∇ K←NRAN R;N;Z
[1] K←(¯0.5+?(Z←2,⌈(N←×/,R)÷2)⍴10000000)÷10000000 ⍝ BOX-MULLER SIM.
[2] K←R⍴N↑((¯2×⍟(×/Z)⍴K[1;])*0.5)×, 1 2 ∘.○○2×K[2;] ⍝ R-SHAPED OUTPUT
∇
so we use a downright ancient technique to relate random contingency tables to random normals, drawing on the linked papers of Irwin (1949) and Lancaster (1949).
The condition of having fixed row and column totals can be expressed by using fixed row (or column) probabilities \(p_1,\ldots,p_n\) where \(p_i\) will be the ith row total divided by the grand total, and so on. The expected frequencies in the ith row and jth column will be \(e_{ij} = \text{grand total} \times \text{row probability} \times \text{column probability}\), and if the observed count in that cell is \(n_{ij}\) then \((n_{ij} - e_{ij})/\sqrt{e_{ij}} = q_{ij}\) gives a matrix \(Q\) of standardised deviations from the expected frequencies.
In the numerical example, this is
The known row or column probabilities can be used to construct a pair of matrices using the formula
for the row probabilities to get
These matrices are orthogonal, that is, with their inverse equal to their transpose. They are calculated by the APL function
∇ R←ORTHOG P;I;J;N;CP
[1] CP←+\P ⍝ HELMERT TRANSFORMATION MATRICES NEEDED FOR <RANTAB>.
[2] N←⍴P ⍝ COPIED DIRECTLY FROM: IRWIN, 'A NOTE ON THE
[3] R←(N,N)⍴0 ⍝ SUB-DIVISION OF CHI-SQUARED INTO COMPONENTS',
[4] R[1;]←P*0.5 ⍝ BIOMETRIKA, 1949, PP. 130-4.
[5] I←2
[6] K:J←1 ⍝ <P> IS A VECTOR OF ROW OR COLUMN PROBABILITIES.
[7] R[I;I]←-(CP[I-1]÷CP[I])*0.5
[8] L:R[I;J]←(P[I]×P[J]÷(CP[I]×CP[I-1]))*0.5
[9] →((J←J+1)<I)/L
[10] →((I←I+1)≤N)/K
∇
which is an exact transcription.
We then find that
has zeros along the first row and column and independent random normal deviates elsewhere, e.g.
If we generate say 1,000 such sets of random normal deviates and reverse the transformation, we will get 1,000 sets of numbers with the right row totals, column totals, and statistical independence properties. If we round them to whole numbers, our task is done. Hence
∇ Z←NS RANTAB RC;ROW;COL;TOT;R;C;E;N;A;B
[1] ⍝ SIMULATION OF RANDOM R×C CONTINGENCY TABLES.
[2] ⍝ FOR METHOD SEE LANCASTER, 'THE DERIVATION AND PARTITION OF
[3] ⍝ CHI-SQUARE IN CERTAIN DISCRETE DITRIBUTIONS',
[4] ⍝ BIOMETRIKA, 1949, PP117-29.
[5] N←1↑NS ⍝ NUMBER OF SIMULATIONS REQUIRED
[6] A←NS[2] ⍝ NUMBER OF ROWS
[7] B←NS[3] ⍝ NUMBER OF COLUMNS
[8] ROW←A↑RC ⍝ ROW TOTALS
[9] COL←(-B)↑RC ⍝ COLUMN TOTALS
[10] TOT←+/ROW
[11] →(TOT=+/COL)/C ⍝ CHECK CONSISTENCY OF SPECIFICATIONS
[12] 'OVERALL TOTALS INCONSISTENT'
[13] →0
[14] C:R←ORTHOG ROW÷TOT
[15] C←ORTHOG COL÷TOT
[16] E←ROW∘.×COL÷TOT
[17] Z←(2 3 1 ⍉(⍉R)+.×0,[2]0,[1]NRAN(A-1),(B-1),N)+.×C ⍝ SIMULATION
[18] Z←⌊0.5+E+((E←((N,(A-1),(B-1))⍴ ¯1 ¯1 ↓E))*0.5)× 0 ¯1 ¯1 ↓Z ⍝ ROUND
[19] Z←Z,[2]((N,B)⍴COL)-+/[2]Z←Z,[3]((N,(A-1))⍴¯1↓ROW)-+/[3]Z ⍝ FILL
[20] Z←(E←∧/[2]∧/[3]Z≥0)/[1]Z ⍝ CHECK FOR NEGATIVE FREQUENCIES
[21] →(N>+/E)×E ⍝ STOP IF FREQUENCIES ALL POSITIVE, OTHERWISE
[22] E:Z←Z,[1]((N-+/E),A,B)RANTAB RC ⍝ GO BACK FOR MORE
∇
Breaking this down, lines 5 to 13 split the first parameter up into the dimensions of the result NxAxB and the second parameter up into the required row and column totals.
[5] N←1↑NS
[6] A←NS[2]
[7] B←NS[3]
[8] ROW←A↑RC
[9] COL←(-B)↑RC
[10] TOT←+/ROW
[11] →(TOT=+/COL)/C
[12] 'OVERALL TOTALS INCONSISTENT'
[13] →0
The next three lines
generate the non-random matrices \(E\), \(R\), and \(C\) Notice that ORTHOG only has to be used twice, however many simulations are required, so that the presence of loops in this function does not much matter.
The generation of Nx(A − 1) × (B − 1) independent random normals is, of course, a single call on the function NRAN.
The ‘3-dimensional matrix’ Z that results when zeros are glued on in the first row and column.
is of dimension AxBxN in order that it may be pre-multiplied by the AxA matrix \(R^T\). It is then converted to the intended form by
and post-multiplied by the BxB matrix \(C\) , to get a 3-dimensional form of the equation
If we did not have to round the result to whole numbers, we could then take
As it is, we drop the last row and column and round everything that is left to the nearest whole number, giving
The compound of all this is lines 17 and 18
[17] Z←(2 3 1 ⍉(⍉R)+.×0,[2]0,[1]NRAN(A-1),(B-1),N)+.×C
[18] Z←⌊0.5+E+((E←((N,(A-1),(B-1))⍴ ¯1 ¯1 ↓E))*0.5)× 0 ¯1 ¯1 ↓Z
It is now a simple matter to fix the row and column totals once and for all in line 19. Occasionally, this may give rise to a negative entry, which is detected in line 20 and edited out of the list of random tables. If this happens, line 21 sends the function back for more in line 22.
[19] Z←Z,[2]((N,B)⍴COL)-+/[2]Z←Z,[3]((N,(A-1))⍴¯1↓ROW)-+/[3]Z
[20] Z←(E←∧/[2]∧/[3]Z≥0)/[1]Z
[21] →(N>+/E)×E
[22] E:Z←Z,[1]((N-+/E),A,B)RANTAB RC
∇
Unless the expected frequencies are very small, it is unlikely that line 22 will be activated. Thus the more elegant recursive approach has been preferred to a loop, despite the need to repeat lines 5-16.
For any heretics who dislike the idea of ‘3-dimensional matrices’ and would prefer a loop, an alternative version is offered, with the rounding procedure hived off to a separate function.
∇ O←ROW RTABLE1 COL;TOT;R;C;EXP;N;A
[1] 'INPUT NUMBER OF CYCLES'
[2] N←⎕
[3] O←(0,(⍴ROW),(⍴COL))⍴0
[4] TOT←+/ROW
[5] →(TOT=+/COL)/C
[6] 'OVERALL TOTALS INCONSISTENT'
[7] →0
[8] C:R←ORTHOG ROW÷TOT
[9] C←ORTHOG COL÷TOT
[10] EXP←ROW∘.×COL÷TOT
[11] S:A←NRAN((⍴ROW)-1),((⍴COL)-1)
[12] A←0,[1]0,A
[13] A←(⍉R)+.×A+.×C
[14] A←EXP+A×EXP*0.5
[15] A←ROUNDB A
[16] →(⎕←0<(∨/∨/(A<0)))/S
[17] O←O,[1]A
[18] →((N←N-1)>0)/S
∇
∇ R←ROUNDB A
[1] R← ¯1 ¯1 ↓⌊A+0.5
[2] R←R,(¯1↓+/A)-+/R
[3] R←R,[1](+⌿A)-+⌿R
∇
A quad in line 16 is to keep a terminal awake (this version is rather slow) and is entirely optional.
The nearest equivalent program in FORTRAN is explained and listed in Patefield (1981). It would, however, be very difficult to use it as in the following example.
Example¶
⍝ RESULTS OF SALES FROM TOWNS A THROUGH H ARE AS FOLLOWS:
TOWN: A B C D E F G H
COUNCIL SALES: 70 380 201 157 202 208 262 370 | 1850
DOMESTIC SALES: 80 114 215 130 163 186 310 201 | 1399
----------------------------------------|-------
150 494 416 287 365 394 572 571 | 3249
⍝ QUESTION: IS TOWN B RELYING TOO MUCH ON COUNCIL BUYING?
Y←(100 2 8 RANTAB ROWS,COLS)[;1;2]
⌈/Y←Y[⍋Y]
306
⍝ THE ANSWER WOULD APPEAR TO BE 'YES'
⍝ TO SEE WHAT THE UPPER STATISTICAL LIMITS ARE:
Y[90 95 99 100]
293 295 304 306
⍝ JUST TO MAKE ABSOLUTELY SURE:
⌈/Y←(100 2 8 RANTAB ROWS,COLS)[;1;2]
305
⌈/Y←(1000 2 8 RANTAB ROWS,COLS)[;1;2]
320
⍝ A SECOND RUN PRODUCES THE SAME ANSWER AS THE FIRST.
⍝ A RUN OF 1,000 TRIALS PUSHES THE LIMIT UP A BIT,
⍝ BUT TOWN B IS STILL OVER THE TOP.
References¶
Irwin, J.O. (1949), A Note on the Subdivision of \(\chi^2\) into Components. Biometrika 36, 130.
Lancaster, H.O. (1949), The Derivation and Partition of \(\chi^2\) in Certain Discrete Distributions. Biometrika 36 117.
Patefield W. M. (1981) AS 159: An efficient method of generating random RxC tables with given row and column totals. Applied Statistics 30 pp. 91-97.
Dr. A.G. Prys Williams, Management Science of Statistics, University College of Swansea, SWANSEA SA2 8PP