Mathematical Statistics With Resampling And R

This gives you an approximate p-value based on 9999 random permutations. The important thing to notice is that permutation tests control the type I error rate exactly (up to Monte Carlo error) regardless of the underlying distribution. This is their main advantage over parametric tests. The disadvantage is that they can be computationally expensive with large samples and complex test statistics. Working with bootstrap confidence intervals on a ratio of medians from heavily right-skewed data revealed a subtle issue. The BCa interval returned NA for the acceleration parameter. The problem was that my jackknife replicate variances were all identical because the median statistic was constant across most leave-one-out samples when the sample size was small relative to the discreteness of the data. The acceleration factor became undefined. The workaround was to switch to the simple percentile interval and increase the number of bootstrap replicates from 2000 to 10000. With more replicates, the percentile interval stabilized. I also verified the result by comparing it against a normal approximation on the log-transformed data, which gave a very similar interval. The log transformation is the standard approach for ratio data, and it avoids the acceleration estimation problem entirely while still providing a reasonable interval on the original scale after back-transformation.

Common Pitfalls

There are several ways to misapply resampling methods without realizing it. The most common mistake is applying the bootstrap to statistics that are not smooth functions of the empirical distribution. The sample median, for example, is a discontinuous functional. The bootstrap can still work for medians, but convergence is slower, and the resulting intervals can be unreliable with small samples. Smooth statistics like the mean, regression coefficients, and correlation coefficients converge at the standard n^(-1/2) rate. Median-type statistics may converge more slowly. Another frequent error is using the bootstrap on dependent data without accounting for the dependence structure. The standard bootstrap assumes independent observations. Time series data, spatial data, and clustered data violate this assumption. For time series, the block bootstrap or stationary bootstrap preserves local dependence. For clustered data, you resample entire clusters rather than individual observations. A third pitfall is interpreting bootstrap confidence intervals as posterior Bayesian intervals without justification. They are not the same thing. A bootstrap interval is a frequentist object based on repeated sampling from the empirical distribution. A Bayesian credible interval conditions on the observed data and integrates over the parameter space with respect to a prior. The numerical values can be similar in some cases, but the interpretation is different, and conflating them leads to incorrect statements about what the interval means.

Cross-Validation as Resampling

Cross-validation is a resampling method that is often overlooked in statistics courses but is essential for model assessment. K-fold cross-validation partitions your data into K equal subsets, trains on K-1 subsets, and evaluates on the held-out subset. You repeat this K times and average the performance metric. Leave-one-out cross-validation is the special case where K equals the sample size. The cv package and caret make this straightforward:

library(caret)

train_control <- trainControl(method = "repeatedcv",
                              number = 10,
                              repeats = 3)

model - train(y ~ ., 
               data = my_data,
               method = "lm",
               trControl = train_control)

This runs 10-fold cross-validation repeated 3 times, giving 30 total fits. The repeated cross-validation reduces the variance of the performance estimate compared to a single run of 10-fold CV. The trade-off is computational cost. With moderate-sized datasets and simple models, this is usually acceptable. With large datasets or complex models like random forests, you need to balance the number of repeats against runtime. The jackknife is the older sibling of the bootstrap. It works by systematically leaving out one observation at a time and recomputing the statistic. The jackknife estimate of bias is (n-1) times the difference between the mean of the leave-one-out estimates and the original estimate. The jackknife estimate of standard error uses the standard deviation of the leave-one-out estimates scaled by sqrt((n-1)/n). The jackknife is deterministic—you get the same result every time. This can be an advantage for reproducibility. It is also computationally cheaper than the bootstrap for large datasets because it requires exactly n replicates instead of thousands. However, the jackknife is less accurate for statistics that are not smooth, and it does not generalize as easily to more complex resampling schemes.

Here is a jackknife variance estimate for a regression coefficient:

jackknife_se <- function(data, formula) {
  n <- nrow(data)
  coeffs <- matrix(NA, nrow = n, ncol = length(coef(lm(formula, data))))
  colnames(coeffs) <- names(coef(lm(formula, data)))
  
  for (i in 1:n) {
    fit <- lm(formula, data = data[-i, ])
    coeffs[i, ] <- coef(fit)
  }
  
  se - apply(coeffs, 2, sd)
  return(se * sqrt((n - 1) / n))
}

This function computes leave-one-out regression coefficients and returns the jackknife standard error for each coefficient. It is a direct implementation of the basic jackknife procedure. More sophisticated variants like the delete-d jackknife leave out d observations at a time, but these are less commonly used in practice. Resampling methods are not a universal solution. They fail in several scenarios. First, when your sample is too small relative to the complexity of the statistic. With fewer than 20 observations, bootstrap confidence intervals can be very inaccurate. The empirical distribution is a poor approximation of the true sampling distribution with so little data. Second, resampling cannot create information that is not already in your data. If your sample is biased or unrepresentative, no amount of resampling will fix that. Third, for heavy-tailed distributions where the variance is infinite, the bootstrap can converge very slowly or not at all in the standard form. For small samples, consider exact methods when available. For non-representative samples, the problem is design-based, not computational. For heavy-tailed data, consider robust statistics or transformations before applying resampling. The key insight is that resampling amplifies what is already there. It does not correct fundamental problems with the data collection process.

Practical Recommendations

Start with the percentile bootstrap for simple statistics and moderate sample sizes. It is fast to implement and easy to interpret. Move to BCa when your statistic is biased or your distribution is skewed. Use permutation tests when you need exact type I error control and your sample size is manageable. Apply block bootstraps to time series data. Use cross-validation for model assessment rather than trying to extract confidence intervals from it. Always check your results with multiple methods when possible. If the percentile, basic, and BCa intervals give very different answers, something is wrong with your data or your assumptions. Inspect the bootstrap distribution plot. Look for multimodality, heavy tails, or extreme outliers. A visual check takes seconds and can prevent hours of wasted effort on invalid inference. The R ecosystem for resampling methods is mature and well-documented. The boot package documentation alone contains references to the original research papers and practical examples. Spending time reading those examples before writing your own code will save considerable debugging time later. Most importantly, understand the assumptions behind each method and verify that your data satisfies them. Resampling is powerful, but it is not magic.