6
$\begingroup$

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.

$\endgroup$
6
  • 2
    $\begingroup$ What specifically is your alternative hypothesis and what connection does it have with your data? $\endgroup$
    – whuber
    Jan 14 '15 at 19:40
  • $\begingroup$ I want to test whether variable T (column) is relevant to variable G (row). Therefore, H0: T is independent of G and H1: T is dependent to G. If p < alpha (0.05 in my case) I can conclude that T is relevant to G. Then, if p >= alpha and (1 - beta) is high I can conclude that T is irrelevant to G. Finally, if p >= alpha and (1 - beta) is low, I can't conclude that T is irrelevant. $\endgroup$
    – Matthias
    Jan 14 '15 at 21:19
  • 3
    $\begingroup$ That's standard--but it does not give us a quantifiable alternative to work with. When you compute power, you must specify the alternative hypotheses explicitly enough to be able to determine exactly which probability distributions they correspond to. Another way to put it is that you must xpress precisely how the alternative might deviate from independence. BTW, post-hoc power analysis is irrelevant for drawing the conclusions you mention: all that matters is whether $p$ is low enough for you. $\endgroup$
    – whuber
    Jan 14 '15 at 23:39
  • 4
    $\begingroup$ Put another way, if you want to use the observed data to set your alternative hypothesis for the power calculation, then the power in your example is equal to 1, since power is P(reject H0). Post hoc power is a silly exercise; use power only in planning future studies so that. They have ASA good chance of being able to detect an effect of clinical significance. $\endgroup$
    – Russ Lenth
    Jan 15 '15 at 0:05
  • 3
    $\begingroup$ I think you misunderstand the comments. In pwr.chisq.test, you must specify an effect size (w) to obtain the power; where w itself is computed from the cell probabilities under the alternative. Similarly, when computing power for a Fisher test, you must specify the precise alternative you mean. $\endgroup$
    – Glen_b
    Jan 15 '15 at 11:44
8
$\begingroup$

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.

$\endgroup$
9
  • 2
    $\begingroup$ "We take your existing dataset as the alternative hypothesis / model of the true data generating process" - I feel like there should be some kind of red flashing warning light around this! Although without further information (for example, on what effect size would be of practical significance and therefore we would care about our power to detect) it is the obvious choice. I generally query whether most people who seek an answer to this question actually understand what, philosophically, is the meaning of the number they're trying to obtain . $\endgroup$
    – Silverfish
    Jan 16 '15 at 22:56
  • 1
    $\begingroup$ Although I concur with the warnings given by @Silverfish, I am upvoting this answer because the last line puts a useful, positive construction on the whole exercise: we do power analyses to design future experiments (and maybe muse a little bit about existing data, but never in a formal way). $\endgroup$
    – whuber
    Jan 16 '15 at 23:16
  • 1
    $\begingroup$ @whuber I thought this was a very useful answer and I gave it my upvote when posting my comment. I'd also happily have upvoted an answer with no code that simply expounded what a numerical result to the OP's question might mean - in particular how it relates to P(Type II error) since that seems to be a common and frequently misunderstood concern. I suspect what many people hope to find is a simple answer to P(type II error) that no function or simulation could provide. $\endgroup$
    – Silverfish
    Jan 16 '15 at 23:38
  • $\begingroup$ Thanks for the comments. @Silverfish, note that the definition of post-hoc power takes the existing dataset as H1; anything else is a different type of power analysis. As I read the question, the OP is asking for post-hoc power. I also try to devote a little discussion to whether post-hoc power is ultimately meaningful, but it is what was asked, as best I can divine. $\endgroup$ Jan 16 '15 at 23:53
  • 1
    $\begingroup$ @gung Yes, my comment was meant to be read as supportive of your choice to infer an intention of post hoc analysis even though the OP didn't explicitly state it. As the comments on the question suggest, not all readers interpreted the question the same way and I think answers taking a different angle would have been valid too. But this was a fine answer. $\endgroup$
    – Silverfish
    Jan 17 '15 at 0:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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