Computing Fisher's Exact Test for rxc tables

Question

I need to implement Fisher's Exact Test in R for rxc tables (c > 2). I know there already is fisher.test to compute such a thing but I need to get it done step by step so I can have custom requirements for p-value computing.

I have tried to first get the 'standard' results of the FET. I found in this article the generalization of a table probability computing using a multiple hypergeometric law :

$$ P = \frac{\prod_i{n_{i.}!}\prod_j{n_{.j}!}}{n!\prod_{ij}{n_{ij}!}} $$

As the tables I will be working on are 2x16 with values ranging from 0 to 20-30, $n!$ is too much for R to compute as is, so I used a natural log transformation so that I could compute the probability using lfactorial instead of factorial. I got the following equation :

$$ ln(P) = \sum_i\ln (n_{i.}!) + \sum_j\ln (n_{.j}!) - \ln (n!) - \sum_{ij}\ln (n_{ij}!) $$

I use r2dtable to get a set number of random table with given marginals (as simulations should be enough and a good compromise), compute the probability for all those tables and sum the probability lesser than my observed table's. I then end up with the following code.

set.seed(42)

# Dummy matrix at real dimensions
test1 <- matrix(floor(runif(32)*20), nrow = 2)

# Dummy 2x2 matrix
test2 <- matrix(floor(runif(4)*20), nrow = 2)

# Function to compute the probability of a table
get.tableP <- function(test){
  col.fact <- sum(lfactorial(colSums(test)))
  row.fact <- sum(lfactorial(rowSums(test)))
  cell.fact <- sum(lfactorial(as.vector(test)))
  tot.fact <- lfactorial(sum(test))

  lP <- col.fact + row.fact - cell.fact - tot.fact
}

# Computing the "cut-off" probability
ref.P <- get.tableP(test1)

# Computing a set number of random matrix with given marginals
p.values <- sapply(r2dtable(1000, rowSums(test1), colSums(test1)), get.tableP)

# Computing p-value
sum(exp(p.values[p.values < ref.P])) + exp(ref.P)

# Results given by R standard function
fisher.test(test1, simulate.p.value=T, B=1000)


# Same with 2x2 table

ref.P <- get.tableP(test2)
p.values <- sapply(r2dtable(1000, rowSums(test2), colSums(test2)), get.tableP)

sum(exp(p.values[p.values < ref.P])) + exp(ref.P)
fisher.test(test2, simulate.p.value=T, B=1000)

My problem is that when comparing my results to those returned by fisher.testI get huge differences when working on my 2x16 tables whereas the 2x2 tables are closer.

Is there anything I got wrong with the math or is it a problem with my implementation ?

Thanks !

Answer 1

Ok, the problem was that r2dtable generates random tables with replacement so when I was trying to test 1000 tables I was really testing about 60 unique tables. That's why such a gap was observed I think.

So the next question is how to make sure to have 1000 unique tables without having to generate 1e8 tables and filter out duplicates ?

Thanks for your help