Suspiciously small p-values from randomization (Fisher-Pittman) test

Question

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.

Answer 1

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).