Computing the power of Fisher's exact test in R

Question

Suppose that I have the following 3 x 3 contingency table:

      |    T1   |    T2   |    T3   |
------+---------+---------+---------+---
  G1  |   18    |   15    |   65    | 
------+---------+---------+---------+---
  G2  |   20    |   10    |   70    |
------+---------+---------+---------+---
  G3  |   15    |   55    |   30    |

Then I can run the Fisher's exact test (using the Monte Carlo simulation option) in R as follows:

table = matrix(c(18,20,15,15,10,55,65,70,30), 3, 3)
fisher.test(table, simulate.p.value=TRUE)

, yielding the following result (p-value):

Fisher's Exact Test for Count Data with simulated p-value (based on
    2000 replicates)

data:  table
p-value = 0.0004998
alternative hypothesis: two.sided

In my hypothesis testing I would also be interested in observing the probability of a type-II error (false-negative rate, beta) as well as the probability of correctly rejecting the null hypothesis (H0) when it is false (statistical power, 1- beta). Does anybody know whether this is possible with R and Fisher's exact test on m x n (with m > 2 and n > 2) contingency tables? Do you maybe know some other tools I can use to calculate this? Thanks!

P.S.: Please note that the above contingency table is just an example and the real obtained p-values can be > alpha (in my case 0.05) which is the more interesting case.

Answer 1

What you are asking for here is a post-hoc power analysis. (More specifically, "the probability of correctly rejecting the null hypothesis" is the power, and 1-power is beta, "the probability of a type-II error". You ask for both, but we only need one to know the other.) We take your existing dataset as the alternative hypothesis / model of the true data generating process. I don't know of a specialized, pre-existing function (e.g., in the pwr package) to do this, but, yes, this can be done in R. You will just have to simulate it. For (considerably) more information on power analyses, and simulating them in R, you should read my answer here: Simulation of logistic regression power analysis - designed experiments. In this case, I will just give a quick, adapted version for dealing with Fisher's exact test. (I usually write code as close to pseudocode as possible so that it may be more widely understood, but because this has the potential of taking so long to run, I try to move as much as possible out of the for loop, and use some of R's unique capacities.)

table = matrix(c(18,20,15,15,10,55,65,70,30), 3, 3)
table
#      [,1] [,2] [,3]
# [1,]   18   15   65
# [2,]   20   10   70
# [3,]   15   55   30
N = sum(table)  # this is the total number of observations
N
# [1] 298
probs = prop.table(table)       
      # these are the probabilities of an observation
probs                           #  being in any given cell
#            [,1]       [,2]      [,3]
# [1,] 0.06040268 0.05033557 0.2181208
# [2,] 0.06711409 0.03355705 0.2348993
# [3,] 0.05033557 0.18456376 0.1006711
probs.v = as.vector(probs)      
      # notice that the probabilities read column-wise
probs.v
# [1] 0.06040268 0.06711409 0.05033557 0.05033557 0.03355705 
    0.18456376 0.21812081
# [8] 0.23489933 0.10067114
cuts = c(0, cumsum(probs.v))    
        # notice that I add a 0 on the front
cuts
# [1] 0.00000000 0.06040268 0.12751678 0.17785235 0.22818792 
      0.26174497
# [7] 0.44630872 0.66442953 0.89932886 1.00000000

set.seed(4941)      # this makes it exactly reproducible
B      = 10000      # number of iterations in simulation
vals   = runif(N*B) # generate random values / probabilities
cats   = cut(vals, breaks=cuts, labels=c("11", "21", "31", 
             "12", "22", "32", "13", "23", "33"))
cats   = matrix(cats, nrow=N, ncol=B, byrow=F)
counts = apply(cats, 2, function(x){ as.vector(table(x)) })

rm(table, N, vals, probs, probs.v, cuts, cats) 
p.vals = vector(length=B) # this will store the outputs
ptm = proc.time()         # this lets me time the simulation
for(i in 1:B){
  mat       = matrix(counts[,i], nrow=3, ncol=3, byrow=T)
  p.vals[i] = fisher.test(mat, simulate.p.value=T)$p.value
}
proc.time() - ptm               # not too bad, really
#  user  system elapsed 
# 28.66    0.32   29.08 
#
mean(p.vals>=.05)               
     # the estimated probability of type II errors is 0
# [1] 0
c(0, 3/B)                       
     # using the rule of 3 to estimate the 95% CI
# [1] 0e+00 3e-04

Given how far your data diverge from the null hypothesis in Fisher's exact test, and the amount of data you have, this simulation does not turn up a single type II error in 10,000 iterations. Because each iteration can be understood as a draw from a binomial distribution with probability $p$ (which we are estimating as the proportion of type II errors observed), this simulation is actually an estimate with some stochastic variability. We can form a 95% confidence interval bounding the true probability of a type II error. To get around the fact that we didn't actually find any type II errors, we will use the rule of 3 ($3/N$) to estimate the upper limit of the CI. Thus, the 95% CI for true type II error rate is $[0,\ 0.0003]$.

---

On a different note, @rvl points out in the comments that "[p]ost hoc power is a silly exercise". That is largely true. I have seen people make the argument, in effect, 'my results are not significant, but I don't have any power, so there's no reason to believe my theory isn't right', which is fairly bizarre on any number of levels. On the other hand, since your results are significant, it isn't clear what difference knowing the post-hoc power for your study is either. I find that understanding post-hoc power can be useful pedagogically to help people begin to understand the topic. And we can also take this as a starting point for a-priori power analyses for planning future studies.