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.


# 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 !


1 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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.