PERMUTATION TEST IN R
Complete Exam Notes Sheet
Assembly Time Example | All Key Concepts | Ready-to-Use Code
1. The Problem Setup
Two groups of 9 employees each assembled a device. We want to test if the NEW procedure takes LESS time on
average than the STANDARD procedure.
Hypothesis
H₀: μ_standard = μ_new (no difference in average times) H₁: μ_new < μ_standard (new procedure is
faster = lower mean) This is a ONE-TAILED test — we only care if new is LOWER.
Step 1: Store Data and Calculate d
standard <- c(32, 37, 35, 28, 41, 44, 35, 31, 34)
new <- c(35, 31, 29, 25, 34, 40, 27, 32, 31)
d <- mean(standard) - mean(new) # observed difference of means
d # should be positive if standard takes longer
WHY d = mean(standard) - mean(new)? If the new procedure is faster, the standard mean will be HIGHER,
making d POSITIVE. We test if d is significantly large (right tail), i.e., is standard significantly slower?
2. Building the Permutation Test
Step 2a: Combine into one vector r
r <- c(standard, new)
# r has 18 values total (9 + 9)
# Under H0, all 18 values are exchangeable between groups
KEY IDEA: Under H₀, group labels don’t matter. So we pool all 18 values and ask: if we randomly split into
two groups of 9, how often do we see a difference as large as d?
Step 2b: Create index
index <- 1:length(r) # = 1:18
# index = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18)
# These are POSITIONS in r, not actual values
We use index (positions) rather than values because combn() works on positions to select subsets. This
way we can always extract r[selected positions] to get actual values.
Step 2c: All Combinations using combn()
p <- combn(index, 9)
# combn(n, k) = choose all combinations of k items from n items
# combn(18, 9) = C(18,9) = 48,620 combinations
# p is a MATRIX:
# - 9 ROWS (each row is one position in the chosen group)
# - 48,620 COLUMNS (each column is one combination)
dim(p) # [1] 9 48620
ncol(p) # 48620 <- number of combinations
nrow(p) # 9 <- group size
length(p) # 437580 <- total elements (9 x 48620) <-- NOT what we loop over!
CRITICAL: ncol(p) vs length(p) vs nrow(p)
ncol(p) = 48,620 → Number of COMBINATIONS (one per column) ✔ USE THIS in loop
nrow(p) = 9 → Group size (how many positions per combination)
length(p) = 437,580 → Total elements in matrix ✘ NEVER use in loop
We loop ncol(p) times because each COLUMN = one combination. Using length(p) would mean 437,580
iterations and p[,i] would crash after column 48,620.
Step 2d: The Loop to Calculate All Differences
dif <- numeric(ncol(p)) # pre-allocate vector of 48,620 zeros
for (i in 1:ncol(p)) {
group1 <- r[p[, i]] # extract 9 values at positions in column i
group2 <- r[-p[, i]] # remaining 9 values (negative indexing)
dif[i] <- mean(group1) - mean(group2)
}
# dif now has 48,620 differences, one per possible split
Code Part What it does
r[p[, i]] Extract values from r at positions stored in column i of p
r[-p[, i]] Extract all values from r EXCEPT those positions (the other group)
numeric(n) Create a vector of n zeros — pre-allocating is faster than growing
mean(group1) - Difference in means for this particular split
mean(group2)
3. Histogram and Visual Analysis
Step 3a: Plot the Histogram
hist(dif,
main = 'Permutation Distribution of Difference in Means',
xlab = 'Difference in Mean Assembly Times (Standard - New)',
ylab = 'Frequency',
col = 'lightgray',
border = 'white')
Step 3b: Add Critical Value Line (Green Dotted)
critical_value <- quantile(dif, 0.95) # 95th percentile for one-tailed 5% test
abline(v = critical_value,
col = 'green',
lty = 3, # lty=3 means DOTTED
lwd = 2)
Step 3c: Add Observed Statistic Line (Blue Dashed)
abline(v = d,
col = 'blue',
lty = 2, # lty=2 means DASHED
lwd = 2)
LINE TYPES (lty) in R: lty = 1 → Solid | lty = 2 → Dashed | lty = 3 → Dotted | lty = 4 → Dot-dash For
this question: Critical value = GREEN DOTTED (lty=3), Observed d = BLUE DASHED (lty=2)
4. Calculating the P-value
Exact Permutation P-value
p_value_exact <- mean(dif >= d)
# This counts: out of 48,620 combinations, what PROPORTION gave
# a difference >= our observed d ?
cat('Exact Permutation p-value:', p_value_exact)
HOW P-VALUE WORKS HERE: mean(dif >= d) gives a TRUE/FALSE vector, then averages it. TRUE = 1,
FALSE = 0, so mean = proportion of cases as extreme or more extreme than observed. This IS the p-value
for a one-tailed (right-tail) test.
Comparison with Normality Assumption
Permutation p-value ≈ 13-14% Normality-based p-value = 13.19% Both are well above 5% → FAIL TO
REJECT H₀ Conclusion: No sufficient evidence that the new procedure is faster. The close agreement
between methods suggests normality assumption is REASONABLE here.
5. Resampling Approach (10,000 Resamples)
Instead of ALL combinations (48,620), we take 10,000 RANDOM shuffles. This is faster and gives a good
approximation.
[Link](77) # ensures reproducibility - same results every run
n_resamples <- 10000
dif_resample <- numeric(n_resamples) # pre-allocate
for (i in 1:n_resamples) {
shuffled <- sample(r) # randomly shuffle all 18 values
group1 <- shuffled[1:9] # first 9 = 'standard'
group2 <- shuffled[10:18] # last 9 = 'new'
dif_resample[i] <- mean(group1) - mean(group2)
}
p_value_resample <- mean(dif_resample >= d)
cat('Resampling p-value:', p_value_resample)
Code Purpose
[Link](77) Fixes random number generator so results are reproducible
sample(r) Randomly shuffles all 18 values in r
shuffled[1:9] Take first 9 shuffled values as 'group 1'
shuffled[10:18] Take remaining 9 shuffled values as 'group 2'
mean(dif_resample >= d) Proportion of resamples with difference >= observed d = p-value
EXACT vs RESAMPLING: Exact: uses ALL C(18,9) = 48,620 combinations → precise but slow
Resampling: uses 10,000 random shuffles → approximate but fast Both should give similar p-values.
Resampling p-value ≈ Exact p-value ≈ 13%
6. Quick Reference Cheat Sheet
Function/Code What it does
combn(index, 9) All ways to choose 9 positions from 18 → C(18,9) = 48,620 columns
ncol(p) Number of combinations = number of loop iterations (48,620)
length(p) Total elements in matrix = 9 × 48,620 = 437,580 (NEVER loop over this)
r[p[, i]] Values at positions in column i of p (one group)
r[-p[, i]] All other values not in column i (other group)
numeric(n) Create vector of n zeros for pre-allocation
mean(dif >= d) P-value: proportion of permutations as extreme as observed
quantile(dif, 0.95) Critical value at 5% significance (one-tailed right)
abline(v=x, lty=2) Vertical dashed line at x on current plot
sample(r) Randomly shuffle vector r
[Link](77) Fix random seed for reproducibility
mean(group1) - Test statistic: difference in group means
mean(group2)
7. Common Exam Mistakes to Avoid
WRONG CORRECT
for (i in 1:length(p)) for (i in 1:ncol(p)) — loop per column/combination
lty=2 for critical value (green) lty=3 (dotted) for critical value, lty=2 (dashed) for observed d
mean(dif > d) for p-value mean(dif >= d) — include equal cases
sample(r, 9) twice for resampling sample(r) then split — avoids overlap
Forgetting [Link]() Always [Link]() before resampling for reproducibility
length(r) gives 18... using it as k in k in combn is 9 (group size), not 18 (total)
combn
8. Complete Code (All in One)
# ============================================================
# PART (i): Data and observed difference
# ============================================================
standard <- c(32, 37, 35, 28, 41, 44, 35, 31, 34)
new <- c(35, 31, 29, 25, 34, 40, 27, 32, 31)
d <- mean(standard) - mean(new)
# ============================================================
# PART (ii): Exact Permutation Test
# ============================================================
r <- c(standard, new) # (a) combined vector
index <- 1:length(r) # (b) positions 1 to 18
p <- combn(index, 9) # (c) all C(18,9) = 48,620 combos
dif <- numeric(ncol(p)) # (d) pre-allocate
for (i in 1:ncol(p)) {
group1 <- r[p[, i]]
group2 <- r[-p[, i]]
dif[i] <- mean(group1) - mean(group2)
}
# ============================================================
# PART (iii): Histogram
# ============================================================
hist(dif, main='Permutation Distribution',
xlab='Difference in Means', col='lightgray')
abline(v = quantile(dif, 0.95), col='green', lty=3, lwd=2) # critical
abline(v = d, col='blue', lty=2, lwd=2) # observed
# ============================================================
# PART (iv): P-value
# ============================================================
p_value_exact <- mean(dif >= d)
cat('Exact p-value:', p_value_exact)
# ============================================================
# PART (v-vi): Resampling with 10,000 samples
# ============================================================
[Link](77)
dif_resample <- numeric(10000)
for (i in 1:10000) {
shuffled <- sample(r)
dif_resample[i] <- mean(shuffled[1:9]) - mean(shuffled[10:18])
}
p_value_resample <- mean(dif_resample >= d)
cat('Resampling p-value:', p_value_resample)
FINAL CONCLUSION: Both p-values ≈ 13% > 5% significance level Fail to reject H₀ — insufficient
evidence that new procedure is faster Normality assumption is reasonable (parametric and non-parametric
agree)