16
$\begingroup$

I have two datasets and I would like to know if they are significantly different or not (this comes from "Two groups are significantly different? Test to use").

I decided to use a permutation test, doing the following in R:

permutation.test <- function(coding, lncrna) {
    coding <- coding[,1] # dataset1
    lncrna <- lncrna[,1] # dataset2

    ### Under null hyphotesis, both datasets would be the same. So:
    d <- c(coding, lncrna)

    # Observed difference
    diff.observed = mean(coding) - mean(lncrna)
    number_of_permutations = 5000
    diff.random = NULL

    for (i in 1:number_of_permutations) {
        # Sample from the combined dataset
        a.random = sample (d, length(coding), TRUE)
        b.random = sample (d, length(lncrna), TRUE)
        # Null (permuated) difference
        diff.random[i] = mean(b.random) - mean(a.random)
    }

    # P-value is the fraction of how many times the permuted difference is equal or more extreme than the observed difference
    pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations
    pvalue
}

Nevertheless, p-values should not be 0 according to this paper: http://www.statsci.org/smyth/pubs/permp.pdf

What do you recommend me to do? Is this way to calculate the p-value:

pvalue = sum(abs(diff.random) >= abs(diff.observed)) / number_of_permutations

a good way? Or is it better doing the following?

pvalue = sum(abs(diff.random) >= abs(diff.observed)) + 1 / number_of_permutations + 1
$\endgroup$
  • $\begingroup$ (1) The final line in the question is erroneous because it does not include the parentheses necessary for it to execute the intended calculation. (It is guaranteed to produce results greater than $1$, which is impossible for any p-value.) (2) You are not actually conducting a permutation test: the two samples a.random and b.random will rarely comprise a random partition of the data but typically will overlap substantially. Instead, compute b.random as the complement of a.random within the union of coding and lncrna. $\endgroup$ – whuber Aug 18 '14 at 16:24
  • 1
    $\begingroup$ Because the p-value is the set of values at least as extreme as the observed, if one evaluates the permutation distribution, the observed statistic is in the counted "permutations". When doing randomization, it's common to count the observed statistic among the considered permutation statistics (for similar reasons). $\endgroup$ – Glen_b Aug 23 '14 at 9:27
17
$\begingroup$

Discussion

A permutation test generates all relevant permutations of a dataset, computes a designated test statistic for each such permutation, and assesses the actual test statistic in the context of the resulting permutation distribution of the statistics. A common way to assess it is to report the proportion of statistics which are (in some sense) "as or more extreme" than actual statistic. This is often called a "p-value."

Because the actual dataset is one of those permutations, its statistic will necessarily be among those found within permutation distribution. Therefore, the p-value can never be zero.

Unless the dataset is very small (less than about 20-30 total numbers, typically) or the test statistic has a particularly nice mathematical form, is not practicable to generate all the permutations. (An example where all permutations are generated appears at Permutation Test in R.) Therefore computer implementations of permutation tests typically sample from the permutation distribution. They do so by generating some independent random permutations and hope that the results are a representative sample of all the permutations.

Therefore, any numbers (such as a "p-value") derived from such a sample are only estimators of the properties of the permutation distribution. It is quite possible--and often happens when effects are large--that the estimated p-value is zero. There is nothing wrong with that, but it immediately raises the heretofore neglected issue of how much could the estimated p-value differ from the correct one? Because the sampling distribution of a proportion (such as an estimated p-value) is Binomial, this uncertainty can be addressed with a Binomial confidence interval.


Architecture

A well-constructed implementation will follow the discussion closely in all respects. It would begin with a routine to compute the test statistic, as this one to compare the means of two groups:

diff.means <- function(control, treatment) mean(treatment) - mean(control)

Write another routine to generate a random permutation of the dataset and apply the test statistic. The interface to this one allows the caller to supply the test statistic as an argument. It will compare the first m elements of an array (presumed to be a reference group) to the remaining elements (the "treatment" group).

f <- function(..., sample, m, statistic) {
  s <- sample(sample)
  statistic(s[1:m], s[-(1:m)])
}

The permutation test is carried out first by finding the statistic for the actual data (assumed here to be stored in two arrays control and treatment) and then finding statistics for many independent random permutations thereof:

z <- stat(control, treatment) # Test statistic for the observed data
sim<- sapply(1:1e4, f, sample=c(control,treatment), m=length(control), statistic=diff.means)

Now compute the binomial estimate of the p-value and a confidence interval for it. One method uses the built-in binconf procedure in the HMisc package:

require(Hmisc)                                    # Exports `binconf`
k <- sum(abs(sim) >= abs(z))                      # Two-tailed test
zapsmall(binconf(k, length(sim), method='exact')) # 95% CI by default

It's not a bad idea to compare the result to another test, even if that is known not to be quite applicable: at least you might get an order of magnitude sense of where the result ought to lie. In this example (of comparing means), a Student t-test usually gives a good result anyway:

t.test(treatment, control)

This architecture is illustrated in a more complex situation, with working R code, at Test Whether Variables Follow the Same Distribution.


Example

As a test, I generated $10$ normally distributed "control" values from a distribution with mean $0$ and $20$ normally distributed "treatment" values from a distribution with mean $1.5$.

set.seed(17)
control <- rnorm(10)
treatment <- rnorm(20, 1.5)

After using the preceding code to run a permutation test I plotted the sample of the permutation distribution along with a vertical red line to mark the actual statistic:

h <- hist(c(z, sim), plot=FALSE)
hist(sim, breaks=h$breaks)
abline(v = stat(control, treatment), col="Red")

Figure

The Binomial confidence limit calculation resulted in

 PointEst Lower        Upper
        0     0 0.0003688199

In other words, the estimated p-value was exactly zero with a (default 95%) confidence interval from $0$ to $0.00037$. The Student t-test reports a p-value of 3.16e-05, which is consistent with this. This supports our more nuanced understanding that an estimated p-value of zero in this case corresponds to a very small p-value that we can legitimately take to be less than $0.00037$. That information, although uncertain, usually suffices to make a definite conclusion about the hypothesis test (because $0.00037$ is far below common thresholds of $0.05$, $0.01$, or $0.001$).


Comments

When $k$ out of $N$ values in the sample of the permutation distribution are considered "extreme," then both $k/N$ and $(k+1)/(N+1)$ are reasonable estimates of the true p-value. (Other estimates are reasonable, too.) Normally there is little reason to prefer one to the other. If they lead to different decisions, that means $N$ is too small. Take a larger sample of the permutation distribution instead of fudging the way in which the p-value is estimated.

If greater precision in the estimate is needed, just run the permutation test longer. Because confidence interval widths typically scale inversely proportional to the square root of the sample size, to improve the confidence interval by a factor of $10$ I ran $10^2=100$ times as many permutations. This time the estimated p-value was $0.000005$ (five of the permutation results were at least as far from zero as the actual statistic) with a confidence interval from $1.6$ through $11.7$ parts per million: a little smaller than the Student t-test reported. Although the data were generated with normal random number generators, which would justify using the Student t-test, the permutation test results differ from the Student t-test results because the distributions within each group of observations are not perfectly normal.

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ The paper by Smyth & Phipson cited above clearly shows why k/N is a poor choice for a p-value estimator. In a nutshell, for relevant significance levels such as alpha=0.05, P((k/N)<alpha|H0) can be surprisingly greater than alpha. This means that a random permutation test using k/N as its p-value estimator and 0.05 as its rejection threshold will reject the null hypothesis more than 5% of the times! A zero p-value is an extreme case of this problem - with a criterion of alpha=0 we expect to never reject the null, yet b/m can equal zero under the null, leading to a false rejection. $\endgroup$ – Trisoloriansunscreen Aug 19 '14 at 12:49
  • 1
    $\begingroup$ @Tal "A poor choice" for a particular purpose. What distinguishes us as statisticians from others is our understanding of the role of variability in data analysis and decision making, along with our ability to quantify that variability appropriately. That is the approach exemplified (and implicitly advocated) in my answer here. When it is carried out there is no such problem as you describe, because the user of the permutation procedure is led to understand its limitations and its strengths and will have the freedom to act according to his or her objectives. $\endgroup$ – whuber Aug 19 '14 at 16:52
14
$\begingroup$

Since estimated p-values are used in order to decide whether to reject the null hypothesis, it is important to consider how the choice of estimator affects the probability of a false rejection. The cited paper by Smyth & Phipson's points out that the unbiased estimator ($\frac{B}{M}$) fails to control the type-I error rate correctly. In contrast, ($\frac{B+1}{M+1}$) is a valid (but conservative) p-value estimator - it doesn't lead to excess rejection of the null.

(B is the number of random permutations in which a statistic greater or equal than the observed one is obtained and M is the total number of random permutations sampled).

Smyth & Phipson also demonstrate that the invalidity of ($\frac{B}{M}$) becomes critical in multiple comparisons settings, where very small p-value estimates are derived and then corrected by multiplication with a factor. An estimation of a zero p-value under the null is especially disastrous in these settings, since it stays zero regardless of the corrections applied.

| cite | improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ +1 This is a good summary of the paper's main point. I especially appreciate your attention to the distinction between an estimated p-value and the true permutation p-value. $\endgroup$ – whuber Aug 19 '14 at 16:54

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.