1.1 The problem with parametric assumptions

Consider measuring dissolved oxygen at 8 stream sites. You want a 95% confidence interval for the median DO. The t-interval applies to means, not medians. There is no simple formula for the SE of a median. What do you do?

Resampling answers: use the data to simulate the sampling process.

1.2 Bootstrapping — resampling with replacement

The idea: your observed sample of \(n\) values is the best approximation you have of the population. Repeatedly draw new samples of size \(n\) with replacement from your observed data, compute your statistic each time, and use the distribution of those statistics as an approximation of the true sampling distribution.

## Original sample (n=8): 7.2 6.8 8.1 5.9 7.4 6.3 8.8 7

## One bootstrap resample: 7.2 7.4 7.2 7.2 6.8 5.9 6.8 6.8

## Another resample:       7.2 7 8.8 7 8.8 5.9 7.2 7.4

1.3 Permutation testing — resampling without replacement

The idea: if the null hypothesis is true (e.g., group A and group B have the same mean), then the group labels are arbitrary — any assignment of observations to groups is equally likely. Randomly shuffle the labels many times, computing the test statistic each time, to build the null distribution. Compare your observed statistic to this null distribution.

## Observed mean difference: 3.4 mg/L

## One permutation:

##   Shuffled A: 8.2 8.9 11.4   mean = 9.5

##   Shuffled B: 9.4 12.1 13.2   mean = 11.57

2.1 CI for the median — a case where formulas fail

We have 15 stream temperature readings (deg C) from a summer survey. The distribution is right-skewed, so the median is more appropriate than the mean. There is no standard formula for the SE of a median.

set.seed(7)
stream_temp <- round(c(
  rnorm(12, mean = 18, sd = 2),
  c(24.1, 26.3, 25.8)          # three hot-day outliers
), 1)

cat("n =", length(stream_temp), "\n")

## n = 15

cat("Mean:  ", round(mean(stream_temp), 2), "deg C\n")

## Mean:   20.03 deg C

cat("Median:", median(stream_temp), "deg C\n")

## Median: 18.7 deg C

set.seed(123)
B        <- 5000
boot_med <- replicate(B, median(sample(stream_temp, replace = TRUE)))

ci_pct <- quantile(boot_med, c(0.025, 0.975))
cat("Bootstrap 95% CI for median:", round(ci_pct, 2), "deg C\n")

## Bootstrap 95% CI for median: 16.6 22.6 deg C

cat("Observed median:", median(stream_temp), "deg C\n")

## Observed median: 18.7 deg C

hist(boot_med,
     breaks = 40, col = "#4a7c5e", border = "white",
     main   = "Bootstrap Distribution of the Median",
     xlab   = "Bootstrap median (deg C)",
     ylab   = "Frequency")
abline(v = ci_pct, col = "#c8832a", lty = 2, lwd = 2)
abline(v = median(stream_temp), col = "#1a1a18", lwd = 2)
legend("topright",
       legend = c("Observed median", "95% CI"),
       col    = c("#1a1a18", "#c8832a"),
       lty    = c(1, 2), lwd = 2, bty = "n", cex = 0.85)

Bootstrap distribution of the sample median (B = 5000 resamples). The histogram approximates the sampling distribution of the median; the dashed lines mark the 95% percentile CI.

2.2 The percentile method vs. the BCa method

The simplest CI is the percentile method: take the 2.5th and 97.5th percentiles of the bootstrap distribution directly. It works well when the bootstrap distribution is symmetric and unbiased.

For skewed statistics or small samples, the bias-corrected and accelerated (BCa) interval is more accurate. The boot package implements this automatically.

library(boot)

# boot() requires a function(data, indices) that returns the statistic
median_stat <- function(data, indices) median(data[indices])

boot_obj <- boot(stream_temp, statistic = median_stat, R = 5000)
boot_obj

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = stream_temp, statistic = median_stat, R = 5000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     18.7 0.60514    1.813746

# BCa confidence interval
boot.ci(boot_obj, type = "bca")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 5000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boot_obj, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   (16.6, 22.6 )  
## Calculations and Intervals on Original Scale

2.3 Bootstrap CI for a regression slope

The bootstrap is particularly valuable for regression when residuals are non-normal or samples are small. Here we use PM2.5 data to estimate the relationship between temperature and ozone, with bootstrap CIs on the slope.

data(airquality)
aq <- airquality[complete.cases(airquality[, c("Ozone", "Temp")]), ]

# Observed slope
lm_obs   <- lm(Ozone ~ Temp, data = aq)
slope_obs <- coef(lm_obs)["Temp"]

# Bootstrap: resample rows
set.seed(99)
boot_slope <- replicate(5000, {
  idx  <- sample(nrow(aq), replace = TRUE)
  coef(lm(Ozone ~ Temp, data = aq[idx, ]))["Temp"]
})

ci_slope <- quantile(boot_slope, c(0.025, 0.975))
ci_param <- confint(lm_obs)["Temp", ]

cat("Observed slope:          ", round(slope_obs, 3), "\n")

## Observed slope:           2.429

cat("Bootstrap 95% CI:        ", round(ci_slope, 3), "\n")

## Bootstrap 95% CI:         2.086 2.804

cat("Parametric 95% CI (lm):  ", round(ci_param, 3), "\n")

## Parametric 95% CI (lm):   1.967 2.891

hist(boot_slope,
     breaks = 50, col = "#4a7c5e", border = "white",
     main = "Bootstrap Distribution: Ozone ~ Temp Slope",
     xlab = "Slope (ppb per deg F)", ylab = "Frequency")
abline(v = ci_slope,    col = "#c8832a", lty = 2, lwd = 2)
abline(v = slope_obs,   col = "#1a1a18", lwd = 2)
abline(v = ci_param,    col = "#1565c0", lty = 3, lwd = 1.5)
legend("topright",
       legend = c("Observed", "Bootstrap 95% CI", "Parametric 95% CI"),
       col    = c("#1a1a18", "#c8832a", "#1565c0"),
       lty    = c(1, 2, 3), lwd = 2, bty = "n", cex = 0.8)

Bootstrap distribution of the OLS slope for Ozone ~ Temp in the airquality dataset. The dashed lines show the 95% percentile CI. Compare with the parametric CI from confint(lm(…)).

3.1 Bootstrap CI for Pearson r

Pearson correlation is sensitive to outliers and assumes bivariate normality. The bootstrap CI sidesteps both issues.

aq_cc <- airquality[complete.cases(airquality[, c("Wind", "Ozone")]), ]
r_obs <- cor(aq_cc$Wind, aq_cc$Ozone)

set.seed(55)
boot_r <- replicate(5000, {
  idx <- sample(nrow(aq_cc), replace = TRUE)
  cor(aq_cc$Wind[idx], aq_cc$Ozone[idx])
})

ci_r <- quantile(boot_r, c(0.025, 0.975))
cat("Observed r:", round(r_obs, 3), "\n")

## Observed r: -0.602

cat("Bootstrap 95% CI:", round(ci_r, 3), "\n")

## Bootstrap 95% CI: -0.695 -0.492

# Compare with Fisher z parametric CI
n_r   <- nrow(aq_cc)
z_r   <- atanh(r_obs)
se_z  <- 1 / sqrt(n_r - 3)
ci_z  <- tanh(z_r + c(-1, 1) * 1.96 * se_z)
cat("Fisher z 95% CI:", round(ci_z, 3), "\n")

## Fisher z 95% CI: -0.706 -0.471

hist(boot_r,
     breaks = 40, col = "#4a7c5e", border = "white",
     main = "Bootstrap Distribution: cor(Wind, Ozone)",
     xlab = "Pearson r", ylab = "Frequency")
abline(v = ci_r,  col = "#c8832a", lty = 2, lwd = 2)
abline(v = r_obs, col = "#1a1a18", lwd = 2)

Bootstrap distribution of Pearson r between Wind and Ozone in the airquality dataset. Skewness in the bootstrap distribution suggests the Fisher z parametric CI may differ.

3.2 Bootstrap CI for a ratio

Ratios are common in environmental science — N:P ratios, species diversity indices, per-capita emissions. Their sampling distributions are rarely normal. The bootstrap handles them without any transformation.

# Ratio of mean upstream to mean downstream nitrate (simulated)
set.seed(14)
upstream   <- rnorm(12, mean = 0.8, sd = 0.3)  # mg/L NO3
downstream <- rnorm(12, mean = 2.1, sd = 0.8)

ratio_obs <- mean(upstream) / mean(downstream)
cat("Observed ratio (up/down):", round(ratio_obs, 3), "\n")

## Observed ratio (up/down): 0.407

boot_ratio <- replicate(5000, {
  i <- sample(12, replace = TRUE)
  mean(upstream[i]) / mean(downstream[i])
})

ci_ratio <- quantile(boot_ratio, c(0.025, 0.975))
cat("Bootstrap 95% CI:", round(ci_ratio, 3), "\n")

## Bootstrap 95% CI: 0.331 0.496

4.1 Logic of the permutation test

A permutation test answers: if the null hypothesis were true and group membership were arbitrary, how often would we see a difference as large as the one we observed just by chance?

The p-value is simply:

\[p = \frac{\text{number of permutations with } |T^*| \geq |T_{\text{obs}}|}{B}\]

No distributional assumptions. No minimum sample size. Valid for any test statistic you can compute.

4.2 Two-sample permutation test

Stream nitrate concentrations upstream and downstream of an agricultural field. The t-test requires normality; with n = 7 and n = 8, that is hard to verify. The permutation test makes no such requirement.

set.seed(42)
upstream_no3   <- c(0.42, 0.38, 0.51, 0.44, 0.47, 0.39, 0.45)
downstream_no3 <- c(1.82, 2.10, 1.95, 2.24, 1.88, 2.01, 1.77, 2.15)

obs_diff <- mean(downstream_no3) - mean(upstream_no3)
cat("Observed mean difference:", round(obs_diff, 3), "mg/L\n")

## Observed mean difference: 1.553 mg/L

# Pool all observations
all_no3 <- c(upstream_no3, downstream_no3)
n_up    <- length(upstream_no3)
n_dn    <- length(downstream_no3)
n_tot   <- n_up + n_dn

# Permutation distribution
B <- 9999
perm_diffs <- replicate(B, {
  perm   <- sample(all_no3)
  mean(perm[(n_up + 1):n_tot]) - mean(perm[1:n_up])
})

p_perm <- (sum(abs(perm_diffs) >= abs(obs_diff)) + 1) / (B + 1)
cat("Permutation p-value (two-tailed):", round(p_perm, 4), "\n")

## Permutation p-value (two-tailed): 1e-04

cat("t-test p-value (reference):      ",
    round(t.test(downstream_no3, upstream_no3)$p.value, 4), "\n")

## t-test p-value (reference):       0

hist(perm_diffs,
     breaks = 60, col = "#8aad96", border = "white",
     main = "Permutation Null Distribution: Downstream - Upstream Nitrate",
     xlab = "Permuted mean difference (mg/L)", ylab = "Frequency")
abline(v = c(-obs_diff, obs_diff), col = "#c8832a", lty = 2, lwd = 2)
legend("topright",
       legend = c(paste0("Observed |diff| = ", round(obs_diff, 2)),
                  paste0("p = ", round(p_perm, 4))),
       col    = c("#c8832a", NA), lty = c(2, NA),
       bty    = "n", cex = 0.85)

Permutation null distribution for the difference in mean nitrate between downstream and upstream sites (B = 9999 permutations). The observed difference (dashed amber line) falls far into the tail, giving a small p-value.

4.3 Permutation test for correlation

Does PM2.5 correlate with temperature? A permutation test on r directly answers the question without the bivariate normality assumption.

data(airquality)
aq_c <- airquality[complete.cases(airquality[, c("Temp", "Ozone")]), ]

obs_r <- cor(aq_c$Temp, aq_c$Ozone)
cat("Observed r:", round(obs_r, 3), "\n")

## Observed r: 0.698

set.seed(77)
perm_r <- replicate(9999,
  cor(sample(aq_c$Temp), aq_c$Ozone)   # shuffle Temp, keep Ozone fixed
)

p_r <- (sum(abs(perm_r) >= abs(obs_r)) + 1) / (9999 + 1)
cat("Permutation p-value:", round(p_r, 4), "\n")

## Permutation p-value: 1e-04

cat("cor.test p-value:   ", round(cor.test(aq_c$Temp, aq_c$Ozone)$p.value, 4), "\n")

## cor.test p-value:    0

4.4 Randomisation test for ANOVA (multiple groups)

For three or more groups, permute the group labels and compute the F-statistic each time to build a null distribution.

# Divide airquality into three seasonal blocks
aq_c$season <- cut(aq_c$Month,
                   breaks = c(4, 6, 8, 9),
                   labels = c("May-Jun", "Jul-Aug", "Sep"),
                   include.lowest = TRUE)
aq_s <- aq_c[!is.na(aq_c$season), ]

obs_F <- summary(aov(Ozone ~ season, data = aq_s))[[1]]["season", "F value"]
cat("Observed F:", round(obs_F, 3), "\n")

## Observed F: 17.197

set.seed(33)
perm_F <- replicate(4999, {
  aq_s$season_perm <- sample(aq_s$season)
  summary(aov(Ozone ~ season_perm, data = aq_s))[[1]]["season_perm", "F value"]
})

p_F <- (sum(perm_F >= obs_F) + 1) / (4999 + 1)
cat("Permutation p-value:", round(p_F, 4), "\n")

## Permutation p-value: 2e-04

cat("Parametric ANOVA p: ",
    round(summary(aov(Ozone ~ season, data = aq_s))[[1]]["season", "Pr(>F)"], 4), "\n")

## Parametric ANOVA p:  0

hist(perm_F,
     breaks = 50, col = "#8aad96", border = "white",
     main = "Permutation F-test: Ozone by Season",
     xlab = "Permuted F-statistic", ylab = "Frequency")
abline(v = obs_F, col = "#c8832a", lty = 2, lwd = 2)

Permutation null distribution of the F-statistic for comparing ozone across three seasonal groups. The observed F (amber dashed) is well into the tail.