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

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.

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

• 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? – Matt Krause Apr 5 '16 at 15:11
• 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. – Roland Apr 5 '16 at 15:27
• 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). – Matt Krause Apr 6 '16 at 19:19
• 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. – Roland Apr 7 '16 at 6:52
• 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. – Matt Krause Apr 7 '16 at 13:09