Chapter 13 Randomisation Tests

When we observe a difference between two groups, we need to ask whether it provides evidence of a real effect or whether it could plausibly have arisen through random variation. A randomisation test approaches this question by generating results that we might obtain if the null hypothesis were true.

We repeatedly rearrange observations or labels in ways that are consistent with the null hypothesis and calculate a test statistic each time. The resulting values form a null distribution. We then compare the observed test statistic with this distribution. In this chapter, the test statistic is a difference between means, although other statistics can also be used.

Randomisation tests have a long history and were initially proposed by R. A. Fisher in the 1930s. Modern computers make it straightforward to repeat the required rearrangements many times, but computation is not what makes the test valid. The randomisation procedure must match the null hypothesis and respect the study design. If individuals were randomly allocated to control and treatment groups, we can examine other possible allocations. Paired or otherwise dependent observations cannot necessarily be shuffled independently.

The technical term is exchangeability: under the null hypothesis, the observations or labels being shuffled must be interchangeable in the way required by the test. We will apply this idea without developing its mathematical details. The examples also provide further practice with dplyr and ggplot2.

13.1 Randomisation test in R

A new drug has been developed that is supposed to reduce cholesterol levels in men. An experiment has been carried out where 12 human test subjects have been assigned randomly to two groups: “Control” and “Drug”. The pharmaceutical company is hoping that the “Drug” group will have lower cholesterol than the “Control” group. The aim here is to do a randomisation test to check that.

In this example, you will:

  • calculate the observed difference between two group means;
  • generate a null distribution by shuffling treatment labels;
  • use the null distribution to calculate a p-value; and
  • interpret the result in the context of the experiment.

Remember to load the dplyr, tidyr, and ggplot2 packages, and open the book’s RStudio project so that the relative data paths work.

Import the data, called cholesterol.csv.

ch <- read.csv("CourseData/cholesterol.csv")
ch$Treatment <- factor(ch$Treatment, levels = c("Control", "Drug"))

Let’s first take a look at the data by plotting it. I will plot a boxplot and add jittered points for clarity.

ggplot(ch, aes(x = Treatment, y = Cholesterol)) +
  geom_boxplot() +
  geom_jitter(colour = "grey", width = .1)

It looks like there might be a difference between the groups. Now let’s consider our test statistic and hypotheses. The test statistic is the mean cholesterol level in the drug group minus the mean in the control group. Negative differences therefore support the directional alternative.

The null hypothesis is that treatment assignment does not affect an individual’s cholesterol level.

The alternative hypothesis is that the drug reduces cholesterol levels, so the mean of the drug group is less than the mean of the control group.

13.1.1 Calculate the observed difference

There are a few ways of doing this. In base R you can use the function tapply (“table apply”), followed by diff (“difference”).

cholesterolMeans <- tapply(ch$Cholesterol, ch$Treatment, mean)
cholesterolMeans
##  Control     Drug 
## 205.6667 185.5000
diff(cholesterolMeans)
##      Drug 
## -20.16667

Because we are focusing on learning dplyr, you can also calculate the means like this:

ch %>% # ch is the cholesterol data
  group_by(Treatment) %>% # group the data by treatment
  summarise(mean = mean(Cholesterol), .groups = "drop") # calculate means
## # A tibble: 2 × 2
##   Treatment  mean
##   <fct>     <dbl>
## 1 Control    206.
## 2 Drug       186.

Here the pipes (%>%) pass the result of each function on as input to the next. You can use further commands, pull to get the mean vector from the summary table, and then use diff to calculate the difference between the groups, before assigning that to observedDiff.

observedDiff <- ch %>%
  group_by(Treatment) %>% # group the data by treatment
  summarise(mean = mean(Cholesterol), .groups = "drop") %>% # calculate means
  pull(mean) %>% # extract the mean vector
  diff()
observedDiff
## [1] -20.16667

This is a complicated set of commands. To make sure you understand it, try running it bit by bit to see what is going on.

Because the factor levels are explicitly ordered as Control and then Drug, diff() calculates Drug minus Control. The control mean is 205.6667, the drug mean is 185.5000, and the observed difference is -20.1667. The individuals receiving the drug therefore had cholesterol levels approximately 20.2 units lower, on average, than individuals in the control group. This difference alone does not tell us how strong the statistical evidence is.

13.1.2 Null distribution

Now we ask what results we might obtain if the null hypothesis were true. If the drug has no effect, the cholesterol measurements would be unchanged under another allocation of individuals to the two groups. Because the individuals were randomly allocated in this hypothetical experiment, we can explore other possible allocations by shuffling the treatment labels while leaving the measurements unchanged.

The command below uses mutate() and sample() to shuffle the Treatment variable. The 12 cholesterol measurements remain fixed, each shuffled dataset still contains six Control and six Drug labels, and the association between treatment and cholesterol is broken. This represents another allocation consistent with the null hypothesis of no treatment effect.

Let’s try that one time:

ch %>%
  mutate(Treatment = sample(Treatment)) %>%
  # shuffle the Treatment data
  group_by(Treatment) %>%
  summarise(mean = mean(Cholesterol), .groups = "drop") %>%
  pull(mean) %>%
  diff()
## [1] -9.5

In this instance, the difference with the shuffled Treatment values of -9.5 is rather different from our observed difference of -20.167.

There are choose(12, 6) possible ways to select six of the 12 individuals for the Drug group:

choose(12, 6)
## [1] 924
choose(24, 12)
## [1] 2704156

For this small dataset, an exact randomisation test could examine all 924 allocations. With 12 individuals in each of two groups, however, there would already be 2,704,156 possible allocations. A Monte Carlo randomisation test instead examines a random sample of possible allocations. Its p-value is an estimate and therefore varies slightly among runs. We will use this approach to practise sample() and replicate(); the exact result will provide a useful comparison later.

Doing one randomisation is not much help, so we need to repeat the process many times.

We can do this easily in R using replicate(), which is a wrapper that tells R to repeat a command n times and store the results in a vector. Let’s try it first 10 times to see how it works:

replicate(
  10,
  ch %>%
    mutate(Treatment = sample(Treatment)) %>%
    group_by(Treatment) %>%
    summarise(mean = mean(Cholesterol), .groups = "drop") %>%
    pull(mean) %>%
    diff()
)
##  [1]  -9.500000   1.833333   0.500000  -5.166667  -3.166667  11.500000
##  [7] -10.500000  26.833333  -7.833333   6.833333

You can see that replicate() performs the shuffle and recalculates the difference 10 times.

For the main analysis, we use 10,000 randomisations. We store that number as B, so the simulation and p-value calculation use the same value. The data frame contains one shuffled difference from each randomisation.

B <- 10000

shuffledData <- data.frame(
  shuffledDiffs = replicate(
    B,
    ch %>%
      mutate(Treatment = sample(Treatment)) %>%
      group_by(Treatment) %>%
      summarise(mean = mean(Cholesterol), .groups = "drop") %>%
      pull(mean) %>%
      diff()
  )
)

The .groups = "drop" argument tells summarise() that the result no longer needs to remain grouped. This avoids repeated grouping messages without changing a global R option.

13.1.3 Testing significance

Before formally testing the hypothesis, it is useful to visualise the null distribution. The x-axis shows differences between Drug and Control means, and the y-axis shows how frequently those differences occurred among the randomised datasets. Because treatment labels are unrelated to cholesterol under the null randomisation procedure, the exact null distribution has a mean of zero. A finite Monte Carlo sample will usually be centred only approximately on zero.

I use ggplot() to create a plot called p1. Putting the command in brackets both creates the plot object and prints it to the screen.

(p1 <- ggplot(shuffledData, aes(x = shuffledDiffs)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  theme_minimal() +
  labs(
    x = "Difference between means (Drug - Control)",
    y = "Number of randomised datasets",
    title = "Null distribution of differences between means"
  ))

You can now add your observed difference (calculated above) to this plot like this:

p1 + geom_vline(
  xintercept = observedDiff,
  colour = "#B2182B", linewidth = 1.2
)

Where does our observed difference lie relative to the null distribution? Does it appear to be an unusual result?

13.1.4 Testing the hypothesis

Recall that the alternative hypothesis predicts a negative difference: Drug minus Control should be less than zero. This is a one-sided test, so we count randomised differences that are less than or equal to the observed difference.

The p-value tells us how often a result at least as extreme as the observed difference would arise under the null hypothesis, using this test statistic and randomisation procedure.

b <- sum(shuffledData$shuffledDiffs <= observedDiff)
p_value <- (b + 1) / (B + 1)

b
## [1] 325
p_value
## [1] 0.03259674

Here, b is the number of simulated statistics at least as extreme as the observed statistic, and B is the number of randomisations. The +1 correction includes the observed arrangement in the Monte Carlo calculation and prevents an estimated p-value of zero. With this seed, b is 325 and the corrected Monte Carlo estimate is 0.0326. Another seed would give a slightly different estimate.

For comparison, we can enumerate every possible allocation for this small dataset:

possibleDrugGroups <- combn(nrow(ch), 6)

exactDiffs <- apply(possibleDrugGroups, 2, function(drugRows) {
  mean(ch$Cholesterol[drugRows]) - mean(ch$Cholesterol[-drugRows])
})

exact_b <- sum(exactDiffs <= observedDiff)
exact_p_value <- exact_b / length(exactDiffs)

c(extreme_allocations = exact_b, exact_p_value = exact_p_value)
## extreme_allocations       exact_p_value 
##         31.00000000          0.03354978

The exact one-sided result is 31/924, or approximately 0.0335. This is close to the Monte Carlo estimate but is not itself a simulated result.

If the alternative hypothesis had instead been that the drug changes cholesterol in either direction, we would use a two-sided test. For this test statistic, we could count values satisfying:

abs(shuffledData$shuffledDiffs) >= abs(observedDiff)

The alternative hypothesis determines which part of the null distribution counts as at least as extreme.

Before continuing, check that you can answer these questions:

  1. What exactly did we shuffle?
  2. Why was it appropriate to shuffle treatment labels?
  3. What does the null distribution represent?
  4. What does the p-value tell us?
  5. Would a small p-value necessarily imply an important treatment effect?

The result provides evidence against the null hypothesis in the predicted direction. It does not tell us the probability that the null hypothesis is true, and it does not by itself establish that the reduction is biologically or medically important. Conversely, a large p-value would mean that the observed statistic is not particularly unusual under the specified null model; it would not prove that the null hypothesis is true.

13.1.5 Writing it up

We can report the Monte Carlo analysis concisely:

“Mean cholesterol was 205.7 in the control group and 185.5 in the drug group, giving an observed difference of -20.2 (Drug minus Control). A one-sided Monte Carlo randomisation test was used to evaluate whether the drug reduced cholesterol. Treatment labels were randomly reallocated 10,000 times while preserving group sizes, and the difference between group means was recalculated after each allocation. The result provided evidence that cholesterol was lower in the drug group (p = 0.033).”

The p-value describes statistical evidence relative to the null model; it does not measure the size or medical importance of the effect. The estimated reduction is approximately 20.2 units, but interpreting its importance requires relevant biological or medical context. Larger samples often provide greater power to detect effects, but statistical significance and biological importance remain different questions.

13.2 Paired Randomisation Tests

A paired randomisation test is used when the observations in your data are paired or matched in some way. In this type of test, each observation in one group (e.g. Group A) is paired or matched with a corresponding observation in the other group (e.g. Group B), based on some relevant characteristic or factor.

You would use a paired randomisation test when you want to compare the mean values of the two groups while taking into account the paired nature of the observations. This type of test is particularly useful when the paired observations are expected to be more similar to each other than to the observations in the other group.

In the cholesterol example, each individual belonged to one of two independent groups. Here, each individual contributes two linked measurements. We must preserve those links: simply shuffling measurements among subjects would not respect the study design.

Here are three examples of situations where a paired randomisation test would be appropriate:

  1. Pre- and post-treatment measurements: Suppose you are conducting a study to evaluate the effectiveness of a new drug. You measure a specific outcome variable in each participant before and after administering the treatment.

  2. Left-right comparisons: Imagine you are studying the effect of a certain intervention on the strength of participants’ dominant and non-dominant hands. Each participant’s hand strength is measured twice, once for the dominant hand and once for the non-dominant hand.

  3. Matched case-control studies: In a case-control study design, where individuals with a particular condition (cases) are compared to individuals without the condition (controls), it is common to match cases and controls based on certain characteristics such as age, gender, or socioeconomic status. The idea is that the pairs (case vs. control) are similar in every way EXCEPT the condition being investigated.

The paired randomisation test is a one-sample randomisation test in which one difference is calculated for every pair. Here, we define each difference as After minus Before, so positive values represent weight gain. Under an appropriate null hypothesis, reversing the direction of a within-pair difference should be equally plausible. We can therefore multiply each difference randomly by +1 or -1.

This sign-flipping procedure has an assumption of its own: the signs must be exchangeable under the null hypothesis. In a randomised paired experiment, that justification can follow from the design. For observational before-and-after data, it relies more strongly on assumptions such as symmetry of the difference distribution. A mean difference of zero alone does not automatically justify arbitrary sign-flipping.

I will illustrate this with an example from Everitt (1994) who looked at using cognitive behaviour therapy as a treatment for anorexia. Everitt collected data on weights of people before and after therapy. These data are in the file anorexiaCBT.csv

anWide <- read.csv("CourseData/anorexiaCBT.csv")
head(anWide)
##   Subject Week01 Week08
## 1       1   80.5   82.2
## 2       2   84.9   85.6
## 3       3   81.5   81.4
## 4       4   82.6   81.9
## 5       5   79.9   76.4
## 6       6   88.7  103.6

These data are arranged in a so-called “wide” format. For plotting and analysis, they need to be rearranged into a tidy “long” format so that each observation is on a row. We can do this using pivot_longer():

anLong <- anWide %>%
  pivot_longer(
    cols = starts_with("Week"), names_to = "time",
    values_to = "weight"
  ) %>%
  mutate(time = factor(time, levels = c("Week01", "Week08"))) %>%
  arrange(Subject, time)
head(anLong)
## # A tibble: 6 × 3
##   Subject time   weight
##     <int> <fct>   <dbl>
## 1       1 Week01   80.5
## 2       1 Week08   82.2
## 3       2 Week01   84.9
## 4       2 Week08   85.6
## 5       3 Week01   81.5
## 6       3 Week08   81.4

We should always plot the data. So here goes.

(p1 <- ggplot(anLong, aes(x = weight, fill = time)) +
  geom_histogram(
    position = "identity", alpha = .7, bins = 15,
    colour = "white"
  ) +
  labs(x = "Weight", y = "Number of observations", fill = "Time")
)

Another useful way to plot this data is to use an interaction plot. In these plots the matched pairs (grouped by Subject) are joined together with lines. You can plot one like this:

(p2 <- ggplot(anLong, aes(x = time, y = weight, group = Subject)) +
  geom_point() +
  geom_line(linewidth = 1, alpha = 0.5) +
  labs(x = "Time", y = "Weight")
)

The second plot is especially informative because every line joins the two measurements from one subject. What we are interested in is whether weight changed after CBT. The null model for this sign-flipping test treats positive and negative within-subject differences as exchangeable; the directional alternative is that weight increased.

The starting point for the analysis is to calculate the observed change in weight.

anChanges <- anLong %>%
  group_by(Subject) %>%
  summarise(
    change = weight[time == "Week08"] - weight[time == "Week01"],
    .groups = "drop"
  )

You have created a dataset that looks like this:

head(anChanges)
## # A tibble: 6 × 2
##   Subject  change
##     <int>   <dbl>
## 1       1  1.70  
## 2       2  0.700 
## 3       3 -0.1000
## 4       4 -0.700 
## 5       5 -3.5   
## 6       6 14.9

And you can calculate the observed change like this:

obsChange <- mean(anChanges$change)
obsChange
## [1] 3.006897

The observed mean change is approximately 3.0069. This is the average of the 29 within-individual changes, not a difference between two independent samples.

13.2.1 The randomisation test

The logic of this test is that, under the sign-flipping null model, a within-subject difference is equally likely to have its observed sign or the opposite sign.

We randomly multiply each observed difference by +1 or -1. Each individual’s two measurements remain linked, the magnitude of each difference is preserved, and differences are not exchanged between subjects. replace = TRUE is required because every individual’s sign is chosen independently and all combinations of positive and negative signs must be possible.

head(anChanges)
## # A tibble: 6 × 2
##   Subject  change
##     <int>   <dbl>
## 1       1  1.70  
## 2       2  0.700 
## 3       3 -0.1000
## 4       4 -0.700 
## 5       5 -3.5   
## 6       6 14.9
anShuffled <- anChanges %>%
  mutate(sign = sample(c(1, -1),
    size = nrow(anChanges),
    replace = TRUE
  )) %>%
  mutate(shuffledChange = change * sign)

Let’s take a look at this new shuffled dataset:

head(anShuffled)
## # A tibble: 6 × 4
##   Subject  change  sign shuffledChange
##     <int>   <dbl> <dbl>          <dbl>
## 1       1  1.70       1         1.70  
## 2       2  0.700      1         0.700 
## 3       3 -0.1000     1        -0.1000
## 4       4 -0.700     -1         0.700 
## 5       5 -3.5        1        -3.5   
## 6       6 14.9       -1       -14.9

We need to calculate the mean of this shuffled vector. We can do this by pull to get the vector, and then mean.

anChanges %>%
  mutate(sign = sample(c(1, -1),
    size = nrow(anChanges),
    replace = TRUE
  )) %>%
  mutate(shuffledChange = change * sign) %>%
  pull(shuffledChange) %>%
  mean()
## [1] -3.296552

Now we build a null distribution of mean changes by repeating these three steps 10,000 times: choose signs, multiply the differences by those signs, and calculate the mean shuffled difference. We again store the number of randomisations in B.

B <- 10000

nullDist <- data.frame(
  change = replicate(
    B,
    anChanges %>%
      mutate(sign = sample(c(1, -1),
        size = nrow(anChanges),
        replace = TRUE
      )) %>%
      mutate(shuffledChange = change * sign) %>%
      pull(shuffledChange) %>%
      mean()
  )
)

13.2.2 Null distribution

(nullDistPlot <- ggplot(nullDist, aes(x = change)) +
  geom_histogram(bins = 30, fill = "steelblue", colour = "white") +
  theme_minimal() +
  labs(
    x = "Mean change in weight",
    y = "Number of sign randomisations",
    title = "Null distribution of paired mean changes"
  ))

This distribution represents mean changes expected under the sign-flipping null model and is centred approximately on zero. We can add the observed mean change as a line:

nullDistPlot + geom_vline(
  xintercept = obsChange,
  colour = "#B2182B", linewidth = 1.2
)

Where does the observed mean change lie relative to this null distribution?

13.2.3 The formal hypothesis test

The alternative hypothesis predicts weight gain, and After minus Before is positive for weight gain. We therefore use the upper tail, unlike the lower-tail comparison in the cholesterol example.

b_paired <- sum(nullDist$change >= obsChange)
p_value_paired <- (b_paired + 1) / (B + 1)

b_paired
## [1] 152
p_value_paired
## [1] 0.01529847

With this seed, 152 of the 10,000 randomised mean changes were greater than or equal to the observed change. The corrected Monte Carlo p-value is 0.0153, providing evidence that weight increased under the specified null model.

This analysis identifies statistical evidence of a before-and-after change. It does not, by itself, establish that CBT caused the change. For observational data, changes over time, concurrent interventions, and other explanations may also be relevant.

13.3 Exercise: Sexual selection in Hercules beetles

A Hercules beetle is a large rainforest species from South America. Researchers suspect that sexual selection has been operating on the species so that the males are significantly larger than the females. You are given data6 on width measurements in cm from 20 individuals: 10 females and 10 males. Can you use your skills to report whether males are significantly larger than females?

The data are called herculesBeetle.csv and are in the CourseData folder.

Use the preceding examples as a guide, but make the analytical decisions yourself.

Task 1: Explore the data. Import the dataset, examine its structure, visualise body size for males and females, and calculate descriptive statistics with units.

Task 2: Formulate hypotheses. State the null and alternative hypotheses. Decide whether a one-sided or two-sided test is appropriate and explain your decision.

Task 3: Calculate the observed statistic. Calculate the difference between group means and state clearly which mean is subtracted from which.

Task 4: Construct a null distribution. Use randomisation to generate a null distribution with a sufficiently large number of randomisations. Explain what you shuffle and why that shuffle is consistent with the null hypothesis and study design.

Task 5: Evaluate the evidence. Plot the null distribution, indicate the observed statistic, calculate a corrected Monte Carlo p-value, and interpret it in context.

Task 6: Write a short results paragraph. Report the observed group means, difference between means, randomisation method, p-value, and a restrained biological interpretation.

Optional hints

  • Use group_by() and summarise() to calculate group means.
  • Set the factor levels so the direction of diff() is explicit.
  • Use sample() to shuffle the sex labels while preserving the group sizes.
  • Use replicate() to repeat the randomisation, first a few times to test your code and then at least 10,000 times.
  • Use ggplot2 to visualise the null distribution and add the observed statistic with geom_vline().
  • Count simulated statistics in the tail specified by your alternative hypothesis and use the +1 Monte Carlo correction.

Finally, would a statistically significant difference necessarily mean that the difference in body size is biologically important? What additional information might help you decide?

  1. This example is from: https://uoftcoders.github.io/rcourse/lec09-Randomization-tests.html↩︎