A permutation test
Section of Probability and simulation in R, part of Using R
A permutation test
| Type | Section |
| Taught in | Probability and simulation in R |
| Part | Using R |
InsectSprays records insect counts on agricultural units, each sprayed with one of six treatments. Restrict to sprays A and B. The observed contrast is the difference of mean counts. Under a null that says the spray labels are exchangeable, any reassignment of those labels to the same counts is as plausible as the observed assignment. The permutation test compares the observed difference with the differences obtained by shuffling the labels.
sprays_ab <- datasets::InsectSprays |>
dplyr::filter(spray %in% c("A", "B")) |>
dplyr::mutate(spray = droplevels(spray))
observed <- with(sprays_ab, mean(count[spray == "A"]) - mean(count[spray == "B"]))
set.seed(1)
permuted <- purrr::map_dbl(seq_len(1999), function(i) {
label <- sample(sprays_ab$spray)
mean(sprays_ab$count[label == "A"]) - mean(sprays_ab$count[label == "B"])
})
p_mc <- (1 + sum(abs(permuted) >= abs(observed))) / (length(permuted) + 1)
c(observed_mean_difference = observed, monte_carlo_p = p_mc)droplevels drops the unused spray labels so that a later formula does not keep empty levels. sample without a probability argument is a random permutation. Each of the 1999 shuffles is one draw from the null that the labels are exchangeable. purrr::map_dbl insists that every shuffle return one number. A shuffle that accidentally returned a vector would fail here, which is the right failure.
The + 1 in the numerator and the denominator counts the observed arrangement among the reference arrangements. Without it, a contrast more extreme than every shuffle would be reported as probability zero, which this Monte Carlo sample cannot justify. With it, the smallest attainable value is 1/2000.
The simulation approximates the null distribution under exchangeability. It does not make exchangeability true. If the units were arranged in space so that neighbouring plots share a pest load, shuffling labels as if the counts were exchangeable answers a different question from the one the field layout poses. Re-running the shuffles moves the Monte Carlo probability slightly. The seed and the number of shuffles are part of the result. Two thousand shuffles resolve a probability to roughly the nearest half of one thousandth, which is enough to see whether the observed contrast is typical, and not enough to quote many decimal places.
This is a different object from t.test(count ~ spray, data = sprays_ab). The t procedure adds a model for the counts. The permutation procedure uses only the exchangeability of the labels. If those two disagree, the disagreement is the model, not a bug in either command.