4
$\begingroup$

How are p-values calculated when using a Monte-Carlo approximation of the Fisher-Pittman test?


I was under the impression[1] that $p$-values generated by randomization tests should always be of the form $$ p = \frac{k+1}{N_\textrm{rounds}+1}$$ where $N_\textrm{rounds}$ is the total number of times that the data was shuffled into surrogate groups and $k$ is the number of times that the test statistic (e.g., the difference in group means) was at least as large in the surrogate groups as in the actual observations. Alternately, I think I've seen $p = \frac{\max(1, k)}{N_\textrm{rounds}}$ as well, which may just reflect a difference in whether the authors define $N_\textrm{rounds}$ to include the observed data or just the shuffled versions.

However, R's coin package doesn't seem to do this. Instead, it implausibly-small p-values, as in this example:

require('coin')
set.seed(24601)
group <- rep(c(0, 1), each=100)
y <- rnorm(length(group)) + (0.5 * group)
group <- as.factor(group)

oneway_test(y~group, distribution=approximate(B=1000)) # B = # of monte carlo replicates
t.test(y~group)
wilcox_test(y~group)

The Fisher-Pittman test returns a p-value of < 2.2e-16 (i.e., near machine epsilon, and not a multiple of 1/1000), while the other tests return p-values on the order of 1e-5.

It seems like the reported p-value should be something like $p<0.001$ in this case, unless coin/the Fisher-Pittman test has some way of calculating a tighter bound. Does it--or other Monte Carlo approximations--have a way of doing so, or does this reflect a (slightly) broken implementation? For example, I could imagine an algorithm that partitions the data in specific ways, but I have never actually heard of such a thing.


1 Davison AC, Hinkley DV (1997) Bootstrap methods and their application. Cambridge University Press, Cambridge, United Kingdom.

$\endgroup$
2
$\begingroup$

The number of replicates is too low. The calculated p-value is exactly zero. Try increasing the number of replicates:

set.seed(42)
res <- oneway_test(y~group, distribution=approximate(B=1e3))
#   Approximative Two-Sample Fisher-Pitman Permutation Test
#
#data:  y by group (0, 1)
#Z = -4.1669, p-value < 2.2e-16
#alternative hypothesis: true mu is not equal to 0

print(res)

res@distribution@pvalue(-4.1669)
#[1] 0
#99 percent confidence interval:
# 0.000000000 0.005284306 

Note the confidence interval.

set.seed(42)
res <- oneway_test(y~group, distribution=approximate(B=1e6))

print(res)
#   Approximative Two-Sample Fisher-Pitman Permutation Test
#
#data:  y by group (0, 1)
#Z = -4.1669, p-value = 2.1e-05
#alternative hypothesis: true mu is not equal to 0

res@distribution@pvalue(-4.1669)
#[1] 2.1e-05
#99 percent confidence interval:
# 1.106928e-05 3.594601e-05 

My rule of thumb (without taking into account the specific simulation): The number of MC replicates should be at least about one order of magnitude larger than 1/p. Otherwise rare "successes" might not be represented correctly (or even not occur) in the simulation.


To recap:

  1. No "successes" occur during the MC replicates

  2. The p-value is calculated as $p=n_{\text{success}}/n_{\text{replicates}}$, so the calculated p-value is exactly zero.

  3. When printed, p-values that are at or very close to zero are rendered as "p < 2.2e-16." That particular value was chosen because it is machine epsilon (the upper bound on rounding errors).

$\endgroup$
  • $\begingroup$ In a situation where none of the replicates are as extreme as the true value, should it actually be returning zero here? As I wrote at the beginning of my question, I would have expected it to return a slightly-too-conservative $p$-value of 1/(B+1) instead of a wildly-incorrect zero. A warning or diagnostic message wouldn't seem out of place here either. Am I totally off base here? $\endgroup$ – Matt Krause Apr 5 '16 at 15:11
  • $\begingroup$ It's not obvious to me what appropriate conditions for a warning might be. Instead I would prefer if the confidence interval of the p-value (or some other measure of uncertainty) was reported in the print output. $\endgroup$ – Roland Apr 5 '16 at 15:27
  • $\begingroup$ I was thinking something like "No 'successes' in any MC replicate. Consider increasing the number of replications." (which showing the confidence interval would also imply, I guess). $\endgroup$ – Matt Krause Apr 6 '16 at 19:19
  • $\begingroup$ OK. And what about exactly one "success"? We shouldn't really trust such a result either. And I wouldn't feel comfortable with two "successes". Are three "successes" sufficient? Now I would ask for some kind of uncertainty/power analysis. And I assume that the confidence interval is based on such an analysis. $\endgroup$ – Roland Apr 7 '16 at 6:52
  • $\begingroup$ The zero case seems especially pathological. Suppose you have 1/1000 successes, and thus p=1/1000. The binomial CDF (which generates the confidence interval) says that there's only a ~5% chance of seeing 15+ successes out of 10,000 iterations (and likewise for <5), so there is a fairly small chance that the (rounded) p-value will change with more iterations. This is not at all the case when observing 0/1000 successes, as your answer shows. $\endgroup$ – Matt Krause Apr 7 '16 at 13:09

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.