Permutation Tests

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 the null hypothesis is true. Further, random re-samples of the data should show similar lack of patterns. However, if the pattern in the data is real, then random re-samples of the data will show very different patterns from the original.

Consider the following image. In that image, the toy blocks on the left show a clear pattern or structure. They are nicely organized into colored piles. This suggests a real pattern that is not random. Someone certainly organized those blocks into that pattern. The blocks didn’t land that way by random chance. On the other hand, the pile of toy blocks shown on the right is certainly a random pattern. This is a pattern that would result if the toy blocks were put into a bag, shaken up, and dumped out. This is the idea of the permutation test. If there is structure in the data, then “mixing up the data and dumping it out again” will show very different patterns from the original. However, if the data was just random to begin with, then we would see a similar pattern by “mixing up the data and dumping it out again.”

The process of a permutation test is:

Compute a test statistic for the original data.
Re-sample the data (“shake it up and dump it out”) thousands of times, computing a new test statistic each time, to create a sampling distribution of the test statistic.
Compute the p-value of the permutation test as the percentage of test statistics that are as extreme or more extreme than the one originally observed.

In review, 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.

#Step 1 Compute a test statistic for the original data.
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.

#Step 2 Re-sample the data (“shake it up and dump it out”) thousands of times, computing a new test statistic each time, to create a sampling distribution of the test statistic.
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.

#Step 3 Compute the p-value of the permutation test as the percentage of test statistics that are as extreme or more extreme than the one originally observed.
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.

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 <- rep(NA, N) This is a storage container that will be used to store the test statistics from each of the thousands of reuses of the data. 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 (i in 1:N)\{ 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. 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. 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. There is a required closing parenthesis on the for (i in 1:N) statement. 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.
permutedTest <- ...perform test with permutedData... The same test that was performed on the original data, should be performed again on the randomly permuted data.
permutedTestStats[i] <- ...get test statistic... 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]. 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. 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. It is important to remember that the permutation test only changes the way the p-value is calculated. Everything else about the original test is unchanged when switching to a permutation test.

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

# Correct two sided p-value for this study:
2*sum(permutedTestStats <= observedTestStat)/N

## [1] 0.002

Other Examples from Students

Independent Samples t Test (and Wilcoxon) (click to show/hide)

ANOVA (click to show/hide)

Simple Linear Regression (click to show/hide)

Examples

Example 1

# Create the data:
set.seed(1140411)
sample1 <- rnorm(30, 69, 2.5)
sample2 <- rnorm(30, 69, 2.5)
theData <- data.frame(values = c(sample1,sample2), group = rep(c(1,2), each=30))
#View(theData)
boxplot(values ~ group, data = theData)

# Run the permutation test:

myTest <-  t.test(values~group,data = theData,mu=0)
observedTestStat <- myTest$statistic; observedTestStat

##          t 
## -0.7304735

N <- 2000      
permutedTestStats <-  rep(NA, N)
for  (i in 1:N ) {
  permutedData <- sample(x=theData$group)
  permutedTest <- t.test(values~permutedData ,data = theData,mu=0)
  permutedTestStats[i] <- permutedTest$statistic
}
hist(permutedTestStats)
abline(v=observedTestStat)

sum(permutedTestStats >= observedTestStat)/N

## [1] 0.7705

sum(permutedTestStats <= observedTestStat)/N

## [1] 0.2295

Example 2

set.seed(121)
sample1 <- rnorm(30, 185, 8)
sample2 <- sample1 - rnorm(30, 0, 3.5)
theData <- data.frame(values = c(sample1,sample2), group = rep(c(1,2), each=30), id = rep(c(1:30),times=2))
#View(theData)
with(theData, hist(values[group==1] - values[group==2]))

# Perform the permutation test:

myTest <- t.test(values ~ group, data = theData, paired = TRUE, mu = 0)
observedTestStat <- myTest$statistic

N <- 2000      
permutedTestStats <-  rep(NA, N)
for  (i in 1:N ) {
  permutedData <- sample(x=c(1,-1), size=30, replace = TRUE)
  permutedTest <- with(theData, t.test(permutedData*(values[group == 1] - values[group == 2]), mu = 0))
  permutedTestStats[i]  <-  permutedTest$statistic
}
hist(permutedTestStats)
abline(v=observedTestStat)

sum(permutedTestStats >= observedTestStat)/N

## [1] 0.6905

sum(permutedTestStats <= observedTestStat)/N

## [1] 0.3095