A nonparametric approach to computing the p-value for any test statistic.


Overview

In almost all hypothesis testing scenarios, the null hypothesis can be interpreted as follows.

\(H_0\): Any pattern that has been witnessed in the sampled data is simply due to random chance.

Permutation Tests depend completely on this single idea. If all patterns in the data really are simply due to random chance, then random re-samples of the data should show similar patterns. Thus, the process of a permutation test is:

  1. Compute a test statistic for the data.
  2. Reuse the data thousands of times to create a sampling distribution of the test statistic.

The sampling distribution is created by permuting (randomly rearranging) the data thousands of times and calculating a test statistic on each permuted version of the data. A histogram of the test statistics then provides the sampling distribution of the test statistic needed to compute the p-value of the original test statistic.


R Instructions

Any permutation test can be performed in R with a for loop.

myTest <- …perform the initial test… This could be a t.test, wilcox.test, aov, kruskal.test, lm, glm, chisq.test, or any other R code that results in a test statistic. It could even simply be the mean or standard deviation of the data.
observedTestStat <- …get the test statistic… Save the test statistic of your test into the object called observedTestStat. For tests that always result in a single test statistic like a t.test, wilcox.test, kruskal.test, and chisq.test it is myTest$statistic. For an aov, lm, or glm try printing summary(myTest)[] to see what values you are interested in using.

N <- 2000       N is the number of times you will reuse the data to create the sampling distribution of the test statistic. A typical choice is 2000, but sometimes 10000, or 100000 reuses are needed before useful answers can be obtained.
permutedTestStats <-  This is a storage container that will be used to store the test statistics from each of the thousands of reuses of the data. rep(NA, N) The rep() function repeats a given value N times. This particular statement repeats NA’s or “missing values” N times. This gives us N “empty” storage spots inside of permutedTestStats that we can use to store the N test statistics from the N reuses of the data we will make in our for loop.
for The for loop is a programming tool that lets us tell R to run a certain code over and over again for a certain number of times.  (i in   In R, the for loop must be followed by a space, then an opening parenthesis, then a variable name (in this case the variable is called “i”), then the word “in” then a list of values. 1:N The 1:N gives R the list of values 1, 2, 3, … and so on all the way up to N. These values are passed into i one at a time and the code inside the for loop is performed first for i=1, then again for i=2, then again for i=3 and so on until finally i=N. At that point, the for loop ends. ) Required closing parenthesis on the for (i in 1:N) statement. { This bracket opens the code section of the for loop. Any code placed between the opening { and closing } brackets will be performed over and over again for each value of i=1, i=2, … up through i=N.
   Two spaces in front of every line inside of the opening { and closing } brackets helps keep your code organized. permutedData <- …randomly permute the data… This is the most important part of the permutation test and takes some thinking. The data must be randomly reorganized in a way consistent with the null hypothesis. What that means exactly is specific to each scenario. Read the Explanation tab for further details on the logic you should use here.
   Two spaces in front of every line inside of the opening { and closing } brackets helps keep your code organized. permutedTest <- …perform test with permutedData… The same test that was performed on the original data, should be performed again on the randomly permuted data.
   Two spaces in front of every line inside of the opening { and closing } brackets helps keep your code organized. permutedTestStats This is the storage container that was built prior to the for (i in 1:N) code. Inside the for loop, this container is filled value by value using the square brackets [i]. [i] The square brackets [i] allows us to access the “i”th position of permutedTestStats. Remember, since this code is inside of the for loop, i=1 the first time through the code, then i=2 the second time through the code, i=3 the third time through, and so on up until i=N the Nth time through the code.  <- …get test statistic… The test statistic from permutedTest is accessed here and stored into permutedTestStats[i].
} The closing } bracket ends the code that is repeated over and over again inside the for loop.
hist(permutedTestStats) Creating a histogram of the sampling distribution of the test statistics obtained from the reused and permuted data allows us to visually compare the observedTestStat to the distribution of test statistics to visually see the percentage of test statistics that are as extreme or more extreme than the observed test statistic value. This is the p-value.
abline(v=observedTestStat) This adds the observedTestStat to the distribution of test statistics to visually see the percentage of test statistics that are as extreme or more extreme than the observed test statistic value. This is the p-value.
sum(permutedTestStats >= observedTestStat)/N This computes a “greater than” p-value. A two-sided p-value could be obtained by multiplying this value by 2 if the observed test statistic was on the right hand side of the histogram.
sum(permutedTestStats <= observedTestStat)/N This computes a “less than” p-value. A two-sided p-value could be obtained by multiplying this value by 2 if the observed test statistic was on the left hand side of the histogram.


Explanation

The most difficult part of a permutation test is in the random permuting of the data. How the permuting is performed depends on the type of hypothesis test being performed.

Paired Data Example

See the Sleep Paired t Test example for the background and context of the study. Here is how to perform the test as a permutation test instead of a t test.

The question that this sleep data can answer concerns which drug is more effective at increasing the amount of extra sleep an individual receives. The associated hypotheses would be \[ H_0: \mu_d = 0 \] \[ H_a: \mu_d \neq 0 \] where \(\mu_d\) denotes the true mean of the differences between the observations for each drug obtained from each individual. Differences would be obtained by \(d_i = \text{extra}_{1i} - \text{extra}_{2i}\).

To perform a permutation test of the hypothesis that the drugs are equally effective, we use the following code.

# Perform the initial test:
myTest <- with(sleep, t.test(extra[group==1], extra[group==2], paired = TRUE, mu = 0))
# Get the test statistic from the test:
observedTestStat <- myTest$statistic


# Obtain the permutation sampling distribution 
N <- 2000
permutedTestStats <- rep(NA, N)
for (i in 1:N){
  permuteData <- sample(x=c(-1,1), size=10, replace=TRUE) 
  permutedTest <- with(sleep, t.test(permuteData*(extra[group==1] - extra[group==2]), mu = 0))
  #Note, t.test(group1 - group2) is the same as t.test(group1, group2, paired=TRUE).
  permutedTestStats[i] <- permutedTest$statistic
}
hist(permutedTestStats)
abline(v=observedTestStat, col='skyblue', lwd=3)

# Greater than p-value: (not what we want here)
sum(permutedTestStats >= observedTestStat)/N
## [1] 1
# Less than p-value:
sum(permutedTestStats <= observedTestStat)/N
## [1] 0.0025
# Correct two sided p-value for this study:
2*sum(permutedTestStats <= observedTestStat)/N
## [1] 0.005