The nonparametric bootstrap from scratch

resampling

statistics

ecology tutorial

Resample your data to approximate a sampling distribution, and build percentile, basic and BCa confidence intervals by hand for a skewed ecological response ratio.

Author

Tidy Ecology

Published

2026-08-16

Many quantities an ecologist reports have no textbook standard error. A response ratio, a coefficient of variation, a diversity index, a ratio of slopes: each is a smooth function of the data, but the formula for its uncertainty is either ugly or unknown. The bootstrap sidesteps that. It treats the sample as a stand-in for the population and asks a mechanical question: if I could redraw the data many times, how much would my estimate move? We answer it by redrawing from the data we already have.

This post builds the nonparametric bootstrap by hand for a response ratio, the ratio of two group means that turns up constantly in impact studies and meta-analysis. We compute a bootstrap standard error, then three flavours of confidence interval (percentile, basic and bias-corrected-and-accelerated), and finish by checking which one actually covers the truth. The recurring lesson of this series starts here: a resampled interval is only as honest as the assumptions behind the flavour you picked.

The idea in one paragraph

Write the estimator as a function of the sample, theta_hat = t(x). Its sampling distribution (how theta_hat would vary across repeated samples from the population) is what we cannot see, because we have one sample. The bootstrap replaces the unknown population with the empirical distribution: draw a resample of the same size with replacement from the observed data, recompute the statistic, and repeat many times. The spread of those resampled statistics estimates the spread of the real sampling distribution. Nothing about the estimator needs to be linear or normal.

A response ratio with a skewed sampling distribution

Abundance data are typically right-skewed, so we simulate control and impact counts from lognormal distributions and take the ratio of their means. The true ratio is fixed by the population means, so we can judge every interval against it later.

set.seed(34)
mlc <- log(26); slc <- 0.62          # control: meanlog, sdlog
mli <- log(52); sli <- 0.78          # impact: larger and more variable
rr_true <- exp(mli + sli^2/2) / exp(mlc + slc^2/2)   # ratio of population means

nc <- 20; ni <- 22
xc <- rlnorm(nc, mlc, slc)           # control abundances
xi <- rlnorm(ni, mli, sli)           # impact abundances
resp_ratio <- function(a, b) mean(b) / mean(a)
theta_hat <- resp_ratio(xc, xi)

The two samples give control and impact means of 29.5 and 67.3, so the observed response ratio is 2.286 against a true value of 2.237. A single sample can land some way from the truth; that is exactly why we want an interval rather than a point.

Resampling by hand

Each bootstrap replicate resamples within each group with replacement, keeping the two sample sizes fixed, and recomputes the ratio. The standard error is just the standard deviation of those replicates.

set.seed(101)
B <- 4000
boot <- replicate(B, mean(sample(xi, ni, replace = TRUE)) /
                     mean(sample(xc, nc, replace = TRUE)))
se_boot   <- sd(boot)
bias_boot <- mean(boot) - theta_hat
skew_boot <- mean((boot - mean(boot))^3) / sd(boot)^3

The bootstrap standard error is 0.599, and the bootstrap distribution is distinctly right-skewed (skewness 0.76). That skew is the whole reason the interval flavours will disagree: a ratio of positive quantities has a longer upper tail, and any method that forces symmetry will get the shape wrong.

Three intervals, three shapes

The normal interval takes theta_hat plus or minus 1.96 bootstrap standard errors, so it is symmetric by construction. The percentile interval reads the 2.5th and 97.5th percentiles straight off the bootstrap distribution, inheriting its skew. The basic (reverse-percentile) interval reflects those percentiles through the estimate. The BCa interval adjusts the percentiles for two things a bare percentile ignores: median bias (through z0, the fraction of replicates below the estimate, mapped to a normal score) and how fast the standard error changes with the parameter (the acceleration a, from the jackknife skewness).

ci_norm  <- theta_hat + c(-1, 1) * qnorm(0.975) * se_boot
ci_perc  <- as.numeric(quantile(boot, c(0.025, 0.975)))
ci_basic <- as.numeric(2 * theta_hat - quantile(boot, c(0.975, 0.025)))

# BCa: bias-correction z0 and acceleration a (jackknife)
z0   <- qnorm(mean(boot < theta_hat))
jack <- c(sapply(seq_len(nc), function(i) resp_ratio(xc[-i], xi)),
          sapply(seq_len(ni), function(i) resp_ratio(xc, xi[-i])))
jbar <- mean(jack)
acc  <- sum((jbar - jack)^3) / (6 * (sum((jbar - jack)^2))^1.5)
zc   <- qnorm(c(0.025, 0.975))
adj  <- pnorm(z0 + (z0 + zc) / (1 - acc * (z0 + zc)))
ci_bca <- as.numeric(quantile(boot, adj))

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_bar()`).

A histogram of resampled response ratios with a long right tail. Below it, four horizontal interval bars are stacked: the normal interval is centred and symmetric, the percentile and basic intervals are mildly asymmetric, and the BCa interval is shifted furthest to the right. — Figure 1: The bootstrap distribution of the response ratio (4000 resamples) with the point estimate (solid) and true value (dashed). The four intervals beneath it show how the normal interval stays symmetric while the percentile and BCa intervals lean into the upper tail.

The normal interval runs 1.11 to 3.46, with equal arms of 1.17 on each side. The BCa interval runs 1.51 to 4.17, with a short lower arm of 0.77 and a long upper arm of 1.88. The bias correction here is mild (z0 is 0.071, acceleration 0.082), yet it still pushes the interval into the tail where the real uncertainty lives. The basic interval reflects the skew the other way, which is a warning sign: reflecting a skewed distribution can send the interval away from where the mass actually is.

Which interval covers the truth?

Shapes are suggestive; coverage is the test. We know the true ratio, so we can redraw the whole study many times, build each interval on each dataset, and count how often it traps the truth. A 95% interval should succeed 95% of the time.

covsim <- function(R = 1500, Bb = 800) {
  hit <- matrix(0, R, 4); mlo <- matrix(0, R, 4); mhi <- matrix(0, R, 4)
  for (r in seq_len(R)) {
    a <- rlnorm(nc, mlc, slc); b <- rlnorm(ni, mli, sli)
    th <- resp_ratio(a, b)
    bt <- replicate(Bb, mean(sample(b, ni, TRUE)) / mean(sample(a, nc, TRUE)))
    se <- sd(bt)
    c_n <- th + c(-1, 1) * qnorm(0.975) * se
    c_p <- quantile(bt, c(0.025, 0.975))
    c_b <- 2 * th - quantile(bt, c(0.975, 0.025))
    z0b <- qnorm(mean(bt < th))
    jk  <- c(sapply(seq_len(nc), function(i) resp_ratio(a[-i], b)),
             sapply(seq_len(ni), function(i) resp_ratio(a, b[-i])))
    jbb <- mean(jk); ac <- sum((jbb - jk)^3) / (6 * (sum((jbb - jk)^2))^1.5)
    ad  <- pnorm(z0b + (z0b + zc) / (1 - ac * (z0b + zc)))
    c_c <- quantile(bt, ad)
    L <- list(c_n, as.numeric(c_p), as.numeric(c_b), as.numeric(c_c))
    for (k in 1:4) {
      hit[r, k] <- (rr_true >= L[[k]][1] & rr_true <= L[[k]][2])
      mlo[r, k] <- (rr_true <  L[[k]][1])          # missed below
      mhi[r, k] <- (rr_true >  L[[k]][2])          # missed above
    }
  }
  list(cov = colMeans(hit), miss_lo = colMeans(mlo), miss_hi = colMeans(mhi))
}
set.seed(202)
cvr <- covsim(R = 1500, Bb = 800)
cov_norm <- cvr$cov[1]; cov_perc <- cvr$cov[2]; cov_basic <- cvr$cov[3]; cov_bca <- cvr$cov[4]
hi_norm <- cvr$miss_hi[1]; lo_norm <- cvr$miss_lo[1]
hi_basic <- cvr$miss_hi[3]; lo_basic <- cvr$miss_lo[3]
lo_bca <- cvr$miss_lo[4]; hi_bca <- cvr$miss_hi[4]

Two panels. Left: a bar chart of overall coverage by method, all bars below a dashed line at 0.95, basic lowest. Right: grouped bars of miss-below and miss-above for each method; BCa has near-equal small bars, while normal and especially basic have a tiny miss-below bar and a large miss-above bar. — Figure 2: Left: empirical coverage of each 95% interval over 1500 simulated studies, all short of the nominal level (dashed). Right: where each interval misses, split into failures below and above the true ratio. A well-calibrated interval should miss about 2.5% on each side. BCa is nearly balanced; the normal and basic intervals miss almost entirely on the high side because they ignore the right skew.

Every method falls short of 95%: percentile reaches 93.7%, normal 93.4%, BCa 92.0%, and the reflected basic interval is clearly worst at 88.7%. But the overall number hides the more telling failure, which is where each interval misses. A calibrated 95% interval should miss about 2.5% above the truth and 2.5% below. BCa comes closest to that split, missing 3.7% low and 4.3% high. The symmetric normal interval misses high 5.8% of the time against only 0.8% low, and the basic interval is worse still, with almost all of its 11.3% failures (11.1% high versus 0.1% low) on the high side. Both sit too far left because they refuse to lean into the right skew. That imbalance, not the headline coverage, is the real reason to prefer BCa: it splits the risk correctly even when no method reaches the promised 95%.

What to take away

The nonparametric bootstrap is close to free: a few lines of resampling give you a standard error and an interval for an estimator that has no clean formula. The price is that the interval is an approximation whose accuracy depends on the flavour. For a skewed statistic, a symmetric normal interval mislabels the shape, the basic interval can reflect the skew the wrong way, and BCa is worth its small extra width because it corrects for bias and for how the spread changes with the parameter. Whichever you use, remember that reported coverage is optimistic, and always say which interval you computed and how many resamples you used. The next post keeps the resampling idea but changes what we resample from: instead of the data, we resample from a fitted model.

References

Efron 1979 Annals of Statistics 7(1):1-26 (10.1214/aos/1176344552)

Efron 1987 Journal of the American Statistical Association 82(397):171-185 (10.1080/01621459.1987.10478410)

DiCiccio and Efron 1996 Statistical Science 11(3):189-228 (10.1214/ss/1032280214)

Puth, Neuhauser and Ruxton 2015 Journal of Animal Ecology 84(4):892-897 (10.1111/1365-2656.12382)

Efron and Tibshirani 1993 An Introduction to the Bootstrap. ISBN 978-0-412-04231-7

--- title: "The nonparametric bootstrap from scratch" description: "Resample your data to approximate a sampling distribution, and build percentile, basic and BCa confidence intervals by hand for a skewed ecological response ratio." date: "2026-08-16 09:00" categories: [R, resampling, statistics, ecology tutorial] image: thumbnail.png image-alt: "A right-skewed bootstrap distribution of a response ratio with three differently shaped confidence intervals drawn beneath it." --- Many quantities an ecologist reports have no textbook standard error. A response ratio, a coefficient of variation, a diversity index, a ratio of slopes: each is a smooth function of the data, but the formula for its uncertainty is either ugly or unknown. The bootstrap sidesteps that. It treats the sample as a stand-in for the population and asks a mechanical question: if I could redraw the data many times, how much would my estimate move? We answer it by redrawing from the data we already have. This post builds the nonparametric bootstrap by hand for a response ratio, the ratio of two group means that turns up constantly in impact studies and meta-analysis. We compute a bootstrap standard error, then three flavours of confidence interval (percentile, basic and bias-corrected-and-accelerated), and finish by checking which one actually covers the truth. The recurring lesson of this series starts here: a resampled interval is only as honest as the assumptions behind the flavour you picked. ## The idea in one paragraph Write the estimator as a function of the sample, `theta_hat = t(x)`. Its sampling distribution (how `theta_hat` would vary across repeated samples from the population) is what we cannot see, because we have one sample. The bootstrap replaces the unknown population with the empirical distribution: draw a resample of the same size *with replacement* from the observed data, recompute the statistic, and repeat many times. The spread of those resampled statistics estimates the spread of the real sampling distribution. Nothing about the estimator needs to be linear or normal. ## A response ratio with a skewed sampling distribution Abundance data are typically right-skewed, so we simulate control and impact counts from lognormal distributions and take the ratio of their means. The true ratio is fixed by the population means, so we can judge every interval against it later. ```{r} #| label: data set.seed(34) mlc <- log(26); slc <- 0.62 # control: meanlog, sdlog mli <- log(52); sli <- 0.78 # impact: larger and more variable rr_true <- exp(mli + sli^2/2) / exp(mlc + slc^2/2) # ratio of population means nc <- 20; ni <- 22 xc <- rlnorm(nc, mlc, slc) # control abundances xi <- rlnorm(ni, mli, sli) # impact abundances resp_ratio <- function(a, b) mean(b) / mean(a) theta_hat <- resp_ratio(xc, xi) ``` The two samples give control and impact means of `r sprintf("%.1f", mean(xc))` and `r sprintf("%.1f", mean(xi))`, so the observed response ratio is `r sprintf("%.3f", theta_hat)` against a true value of `r sprintf("%.3f", rr_true)`. A single sample can land some way from the truth; that is exactly why we want an interval rather than a point. ## Resampling by hand Each bootstrap replicate resamples *within each group* with replacement, keeping the two sample sizes fixed, and recomputes the ratio. The standard error is just the standard deviation of those replicates. ```{r} #| label: boot set.seed(101) B <- 4000 boot <- replicate(B, mean(sample(xi, ni, replace = TRUE)) / mean(sample(xc, nc, replace = TRUE))) se_boot <- sd(boot) bias_boot <- mean(boot) - theta_hat skew_boot <- mean((boot - mean(boot))^3) / sd(boot)^3 ``` The bootstrap standard error is `r sprintf("%.3f", se_boot)`, and the bootstrap distribution is distinctly right-skewed (skewness `r sprintf("%.2f", skew_boot)`). That skew is the whole reason the interval flavours will disagree: a ratio of positive quantities has a longer upper tail, and any method that forces symmetry will get the shape wrong. ## Three intervals, three shapes The **normal** interval takes `theta_hat` plus or minus `1.96` bootstrap standard errors, so it is symmetric by construction. The **percentile** interval reads the 2.5th and 97.5th percentiles straight off the bootstrap distribution, inheriting its skew. The **basic** (reverse-percentile) interval reflects those percentiles through the estimate. The **BCa** interval adjusts the percentiles for two things a bare percentile ignores: median bias (through `z0`, the fraction of replicates below the estimate, mapped to a normal score) and how fast the standard error changes with the parameter (the acceleration `a`, from the jackknife skewness). ```{r} #| label: intervals ci_norm <- theta_hat + c(-1, 1) * qnorm(0.975) * se_boot ci_perc <- as.numeric(quantile(boot, c(0.025, 0.975))) ci_basic <- as.numeric(2 * theta_hat - quantile(boot, c(0.975, 0.025))) # BCa: bias-correction z0 and acceleration a (jackknife) z0 <- qnorm(mean(boot < theta_hat)) jack <- c(sapply(seq_len(nc), function(i) resp_ratio(xc[-i], xi)), sapply(seq_len(ni), function(i) resp_ratio(xc, xi[-i]))) jbar <- mean(jack) acc <- sum((jbar - jack)^3) / (6 * (sum((jbar - jack)^2))^1.5) zc <- qnorm(c(0.025, 0.975)) adj <- pnorm(z0 + (z0 + zc) / (1 - acc * (z0 + zc))) ci_bca <- as.numeric(quantile(boot, adj)) ``` ```{r} #| label: fig-dist #| echo: false #| fig-cap: "The bootstrap distribution of the response ratio (4000 resamples) with the point estimate (solid) and true value (dashed). The four intervals beneath it show how the normal interval stays symmetric while the percentile and BCa intervals lean into the upper tail." #| fig-alt: "A histogram of resampled response ratios with a long right tail. Below it, four horizontal interval bars are stacked: the normal interval is centred and symmetric, the percentile and basic intervals are mildly asymmetric, and the BCa interval is shifted furthest to the right." #| fig-width: 7 #| fig-height: 4.6 suppressMessages(library(ggplot2)) te_ink <- "#275139"; te_green <- "#4e7a5a"; te_gold <- "#c98a1b" te_rose <- "#a5524a"; te_slate <- "#3f5c7a"; te_grey <- "#8a8a80" theme_te <- function(base_size = 12) { theme_minimal(base_size = base_size) + theme(text = element_text(colour = "#2f3e35"), plot.title = element_text(face = "bold", colour = "#233a2c", size = rel(1.02)), plot.subtitle = element_text(colour = "#54655b", size = rel(0.86)), panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e3d8", linewidth = 0.3), plot.background = element_rect(fill = "#faf7f0", colour = NA), panel.background = element_rect(fill = "#faf7f0", colour = NA), axis.text = element_text(colour = "#54655b"), legend.position = "none", plot.margin = margin(8, 12, 8, 12)) } ints <- data.frame( method = factor(c("normal", "percentile", "basic", "BCa"), levels = c("BCa", "basic", "percentile", "normal")), lo = c(ci_norm[1], ci_perc[1], ci_basic[1], ci_bca[1]), hi = c(ci_norm[2], ci_perc[2], ci_basic[2], ci_bca[2]), col = c(te_slate, te_green, te_grey, te_ink)) ymax <- max(hist(boot, breaks = 40, plot = FALSE)$density) yb <- -ymax * c(0.10, 0.20, 0.30, 0.40) ints$y <- yb[as.integer(factor(ints$method, levels = c("normal","percentile","basic","BCa")))] ggplot() + geom_histogram(aes(x = boot, y = after_stat(density)), bins = 40, fill = "#dfe8e0", colour = "#c4d3c6", linewidth = 0.2) + geom_vline(xintercept = theta_hat, colour = te_ink, linewidth = 0.7) + geom_vline(xintercept = rr_true, colour = te_rose, linewidth = 0.6, linetype = "dashed") + geom_segment(data = ints, aes(x = lo, xend = hi, y = y, yend = y, colour = col), linewidth = 1.6) + geom_point(data = ints, aes(x = theta_hat, y = y, colour = col), size = 1.6) + geom_text(data = ints, aes(x = lo, y = y, label = method, colour = col), hjust = 1.08, size = 3.1) + scale_colour_identity() + scale_x_continuous(limits = c(0.4, 5.2)) + labs(title = "One estimate, four confidence intervals", subtitle = "Solid line: estimate. Dashed line: true ratio. Bars: 95% intervals by method.", x = "Response ratio (impact mean / control mean)", y = "Density") + theme_te() ``` The normal interval runs `r sprintf("%.2f", ci_norm[1])` to `r sprintf("%.2f", ci_norm[2])`, with equal arms of `r sprintf("%.2f", theta_hat - ci_norm[1])` on each side. The BCa interval runs `r sprintf("%.2f", ci_bca[1])` to `r sprintf("%.2f", ci_bca[2])`, with a short lower arm of `r sprintf("%.2f", theta_hat - ci_bca[1])` and a long upper arm of `r sprintf("%.2f", ci_bca[2] - theta_hat)`. The bias correction here is mild (`z0` is `r sprintf("%.3f", z0)`, acceleration `r sprintf("%.3f", acc)`), yet it still pushes the interval into the tail where the real uncertainty lives. The basic interval reflects the skew the *other* way, which is a warning sign: reflecting a skewed distribution can send the interval away from where the mass actually is. ## Which interval covers the truth? Shapes are suggestive; coverage is the test. We know the true ratio, so we can redraw the whole study many times, build each interval on each dataset, and count how often it traps the truth. A 95% interval should succeed 95% of the time. ```{r} #| label: coverage covsim <- function(R = 1500, Bb = 800) { hit <- matrix(0, R, 4); mlo <- matrix(0, R, 4); mhi <- matrix(0, R, 4) for (r in seq_len(R)) { a <- rlnorm(nc, mlc, slc); b <- rlnorm(ni, mli, sli) th <- resp_ratio(a, b) bt <- replicate(Bb, mean(sample(b, ni, TRUE)) / mean(sample(a, nc, TRUE))) se <- sd(bt) c_n <- th + c(-1, 1) * qnorm(0.975) * se c_p <- quantile(bt, c(0.025, 0.975)) c_b <- 2 * th - quantile(bt, c(0.975, 0.025)) z0b <- qnorm(mean(bt < th)) jk <- c(sapply(seq_len(nc), function(i) resp_ratio(a[-i], b)), sapply(seq_len(ni), function(i) resp_ratio(a, b[-i]))) jbb <- mean(jk); ac <- sum((jbb - jk)^3) / (6 * (sum((jbb - jk)^2))^1.5) ad <- pnorm(z0b + (z0b + zc) / (1 - ac * (z0b + zc))) c_c <- quantile(bt, ad) L <- list(c_n, as.numeric(c_p), as.numeric(c_b), as.numeric(c_c)) for (k in 1:4) { hit[r, k] <- (rr_true >= L[[k]][1] & rr_true <= L[[k]][2]) mlo[r, k] <- (rr_true < L[[k]][1]) # missed below mhi[r, k] <- (rr_true > L[[k]][2]) # missed above } } list(cov = colMeans(hit), miss_lo = colMeans(mlo), miss_hi = colMeans(mhi)) } set.seed(202) cvr <- covsim(R = 1500, Bb = 800) cov_norm <- cvr$cov[1]; cov_perc <- cvr$cov[2]; cov_basic <- cvr$cov[3]; cov_bca <- cvr$cov[4] hi_norm <- cvr$miss_hi[1]; lo_norm <- cvr$miss_lo[1] hi_basic <- cvr$miss_hi[3]; lo_basic <- cvr$miss_lo[3] lo_bca <- cvr$miss_lo[4]; hi_bca <- cvr$miss_hi[4] ``` ```{r} #| label: fig-coverage #| echo: false #| fig-cap: "Left: empirical coverage of each 95% interval over 1500 simulated studies, all short of the nominal level (dashed). Right: where each interval misses, split into failures below and above the true ratio. A well-calibrated interval should miss about 2.5% on each side. BCa is nearly balanced; the normal and basic intervals miss almost entirely on the high side because they ignore the right skew." #| fig-alt: "Two panels. Left: a bar chart of overall coverage by method, all bars below a dashed line at 0.95, basic lowest. Right: grouped bars of miss-below and miss-above for each method; BCa has near-equal small bars, while normal and especially basic have a tiny miss-below bar and a large miss-above bar." #| fig-width: 7 #| fig-height: 3.9 methods <- c("normal", "percentile", "basic", "BCa") cols <- c("#3f5c7a", "#4e7a5a", "#8a8a80", "#275139") dcov <- data.frame(method = factor(methods, levels = methods), cov = as.numeric(cvr$cov), col = cols) pA <- ggplot(dcov, aes(method, cov, fill = col)) + geom_col(width = 0.68) + geom_hline(yintercept = 0.95, linetype = "dashed", colour = "#a5524a", linewidth = 0.6) + geom_text(aes(label = sprintf("%.3f", cov)), vjust = -0.5, size = 3, colour = "#233a2c") + scale_fill_identity() + coord_cartesian(ylim = c(0.80, 0.965)) + labs(title = "Overall coverage", subtitle = "Nominal 95% (dashed)", x = NULL, y = "coverage") + theme_te() + theme(axis.text.x = element_text(angle = 20, hjust = 1)) dmiss <- data.frame( method = factor(rep(methods, 2), levels = methods), side = rep(c("misses below", "misses above"), each = 4), rate = c(as.numeric(cvr$miss_lo), as.numeric(cvr$miss_hi))) pB <- ggplot(dmiss, aes(method, rate, fill = side)) + geom_col(position = position_dodge(width = 0.72), width = 0.66) + geom_hline(yintercept = 0.025, linetype = "dotted", colour = "#8a8a80", linewidth = 0.5) + scale_fill_manual(values = c("misses below" = "#c98a1b", "misses above" = "#a5524a")) + labs(title = "Where the misses land", subtitle = "Target: 2.5% each side (dotted)", x = NULL, y = "miss rate", fill = NULL) + theme_te() + theme(axis.text.x = element_text(angle = 20, hjust = 1), legend.position = "bottom") grid::grid.newpage() grid::pushViewport(grid::viewport(layout = grid::grid.layout(1, 2))) print(pA, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 1)) print(pB, vp = grid::viewport(layout.pos.row = 1, layout.pos.col = 2)) ``` Every method falls short of 95%: percentile reaches `r sprintf("%.1f", 100*cov_perc)`%, normal `r sprintf("%.1f", 100*cov_norm)`%, BCa `r sprintf("%.1f", 100*cov_bca)`%, and the reflected basic interval is clearly worst at `r sprintf("%.1f", 100*cov_basic)`%. But the overall number hides the more telling failure, which is *where* each interval misses. A calibrated 95% interval should miss about 2.5% above the truth and 2.5% below. BCa comes closest to that split, missing `r sprintf("%.1f", 100*lo_bca)`% low and `r sprintf("%.1f", 100*hi_bca)`% high. The symmetric normal interval misses high `r sprintf("%.1f", 100*hi_norm)`% of the time against only `r sprintf("%.1f", 100*lo_norm)`% low, and the basic interval is worse still, with almost all of its `r sprintf("%.1f", 100*(1-cov_basic))`% failures (`r sprintf("%.1f", 100*hi_basic)`% high versus `r sprintf("%.1f", 100*lo_basic)`% low) on the high side. Both sit too far left because they refuse to lean into the right skew. That imbalance, not the headline coverage, is the real reason to prefer BCa: it splits the risk correctly even when no method reaches the promised 95%. ## What to take away The nonparametric bootstrap is close to free: a few lines of resampling give you a standard error and an interval for an estimator that has no clean formula. The price is that the interval is an approximation whose accuracy depends on the flavour. For a skewed statistic, a symmetric normal interval mislabels the shape, the basic interval can reflect the skew the wrong way, and BCa is worth its small extra width because it corrects for bias and for how the spread changes with the parameter. Whichever you use, remember that reported coverage is optimistic, and always say which interval you computed and how many resamples you used. The next post keeps the resampling idea but changes what we resample from: instead of the data, we resample from a fitted model. ## Related tutorials - [The parametric bootstrap](../the-parametric-bootstrap/) - [Bootstrapping dependent data](../bootstrapping-dependent-data/) - [Checking a bootstrap](../checking-a-bootstrap/) - [Checking a diversity estimate](../checking-a-diversity-estimate/) ## References Efron 1979 Annals of Statistics 7(1):1-26 (10.1214/aos/1176344552) Efron 1987 Journal of the American Statistical Association 82(397):171-185 (10.1080/01621459.1987.10478410) DiCiccio and Efron 1996 Statistical Science 11(3):189-228 (10.1214/ss/1032280214) Puth, Neuhauser and Ruxton 2015 Journal of Animal Ecology 84(4):892-897 (10.1111/1365-2656.12382) Efron and Tibshirani 1993 An Introduction to the Bootstrap. ISBN 978-0-412-04231-7