[R] power 2x3 exact test

[(Ted Harding)]

On 09-May-07 22:00:27, Duncan Murdoch wrote:
> On 09/05/2007 5:11 PM, Bingshan Li wrote:
>  > Hi, all,
>  >
>  > I am wondering if there is an algorithm for calculating power of 2x3
>  > table using Fisher exact test. There is one (algorithm 280) for 2x2
>  > Fisher exact test but I couldn't find one for 2x3 table. If we are
>  > not lucky enough to have one, is there any other way to calculate
>  > exact power of 2x3 table? The reason why I want exact power is
>  > because some cells are assumed to be very small and chi square
>  > approximation is not valid.
> 
> I think there are lots of possible alternatives to the null in a 2x3 
> table, so you may have trouble finding a single answer to this
> question. 
>   But assuming you have one in mind, I'd suggest doing a Monte Carlo 
> power calculation:  simulate a few thousand tables from the alternative
> distribution, and see what the distribution of p-values looks like.
> 
> Duncan Murdoch

I'd back Duncan on that point!

More specifically, for the 2x2 table, the table, conditional on the
marginals, is a function of one element (say top left-hand corner),
and the probability of any table depends on the single parameter
which is the odds-ratio of the 4 cell probabilities.

so this case is relatively easy and straightforward and, indeed,
for the 2x2 table R'a fisher.test() allows you to specify the
odds-ratio as a "null" parameter.

This is not the case with fisher.test() for a larger (say 2x3
table), so to investigate that case you cannot use fisher.test().

In all cases, however (according to the FORTRAN code on which
itis based -- see the reference in "?fisher.table"), the rejection
region for the exact fisher.test() consists of those table with
the smallest probabilities.

For the 2x3 table, say (cell counts with margins, and probabilities):


   a1  b1  c1 | d1          p1  q1  r1
              |
   a2  b2  c2 | d2          p2  q2  r2
   ---------------
   a   b   c  |  n

so that

   a1+b1+c1 = d1, a2+b2+c2 = d2,
   a1+a2 = a, b1+b2 = b, c1+c2 = c

the table is a function of any two functionally independent cells
(say a1 and b1), and its probability is a function of some two
odds-ratios, say

   (p1*r2)/(r1*p2)

   (q1*r2)/(r1*q2)

which, for the standard null hypothesis, are both equal to 1.
However, as Duncan says, alternatives are 2-dimensional and
so there is not a unique natural form for an alternative (as
opposed  to the 2x2 case, where it boils down to (p1*q2)/(p2*q1)
being not equal to 1, therefore greater than 1, or less than 1,
or 2-sidedly either >1 or <1).

The probability of the 2x3 table is proportional to

  ((p1*r2)/(r1*p2))^a1 * ((q1*r2)/(r1*q2))^b1

(or equivalent), divided by the product of the factorials of
a1, b1, c1, a2, b2, c2, subject to summing to 1 over all
combinations of (a1,b1) giving rise to a table compatible
with the marginal constraints.

Given that you expect some cells to be small, it should not
be a severe task to draw up a list of (a1,b1) values which
correspond to rejection of the null hypothesis (that both
ORs equal 1), and then the simulation using different values
of the two odds-ratios will give you the power for each such
pair of odds-ratios.

The main technical difficulty will be simulation of random
tables, conditional on the marginals, with the probabilities
as given above.

I don't know of a good suggestion for this. The fisher.test()
function will not help (see above). In the case of the 2x2
table, is is a straightforward hypergeometric distribution,
but 2x3 is not straightforward. Maybe someone has written
a function for this kind of application, and can point it
out to us???

A quick R-site search did not help!

Best wishes,
ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at manchester.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 10-May-07                                       Time: 00:12:29
------------------------------ XFMail ------------------------------