Lab 6: Jackknife and BCa Bootstrap

STAT 536 — Statistical Computing

Author

Soumojit Das

Published

Fall 2026

Learning Objectives

By the end of today’s lab, you will:

Implement the jackknife from scratch and read bias, standard error, influence values, and pseudo-values off the leave-one-out estimates
Explain why the jackknife succeeds for the mean and the variance
Diagnose the jackknife’s failure for the median and state it precisely as inconsistency
Implement the BCa interval from scratch, with its bias correction and acceleration, and cross-check it against boot
Apply the jackknife and BCa to the correlation in Efron’s law-school data, and identify the observation that moves the estimate most
Estimate the coverage of the jackknife, percentile, and BCa intervals by simulation, and compare the three

Setup

Figure labels are typeset by LaTeX. tikz_png() draws a ggplot through tikzDevice, compiles it with pdflatex, and returns a PNG. It needs a LaTeX installation (TinyTeX is enough) and installs tikzDevice, pdftools, and png on first use. The boot package is used only to cross-check the BCa interval you write.

library(ggplot2)
library(patchwork)
library(boot)
set.seed(536)
theme_set(theme_bw(base_size = 14))

pal <- c(red = "#D62828", blue = "#2E86AB", green = "#06A77D", dark = "#023047", grey = "#6a6a6a")

find_pdflatex <- function() {
  cand <- c(Sys.which("pdflatex"),
            path.expand("~/Library/TinyTeX/bin/universal-darwin/pdflatex"),
            file.path(Sys.getenv("APPDATA"), "TinyTeX/bin/windows/pdflatex.exe"),
            path.expand("~/.TinyTeX/bin/x86_64-linux/pdflatex"))
  cand <- cand[nzchar(cand) & file.exists(cand)]
  if (!length(cand)) stop("pdflatex not found; install TinyTeX with tinytex::install_tinytex()")
  cand[[1]]
}

tikz_png <- function(p, png_file, width = 9, height = 5, dpi = 150) {
  for (pkg in c("tikzDevice", "pdftools", "png"))
    if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
  options(tikzLatex = find_pdflatex(),
          tikzDocumentDeclaration = r"(\documentclass[10pt]{article})",
          tikzLatexPackages = c(r"(\usepackage{tikz})",
                                r"(\usepackage[active,tightpage,psfixbb]{preview})",
                                r"(\PreviewEnvironment{pgfpicture})",
                                r"(\setlength\PreviewBorder{0pt})",
                                r"(\usepackage[T1]{fontenc})", r"(\usepackage{lmodern})",
                                r"(\usepackage{amsmath})"))
  work <- file.path(tempdir(), "tikz_png"); dir.create(work, showWarnings = FALSE)
  tex <- file.path(work, "fig.tex"); pdf <- file.path(work, "fig.pdf")
  tikzDevice::tikz(tex, width = width, height = height, standAlone = TRUE)
  print(p); dev.off()
  ok <- system2(getOption("tikzLatex"),
                c("-interaction=batchmode", "-halt-on-error",
                  paste0("-output-directory=", shQuote(work)), shQuote(tex)),
                stdout = FALSE, stderr = FALSE)
  if (ok != 0 || !file.exists(pdf)) stop("pdflatex failed; see ", file.path(work, "fig.log"))
  dir.create(dirname(png_file), recursive = TRUE, showWarnings = FALSE)
  png::writePNG(pdftools::pdf_render_page(pdf, page = 1, dpi = dpi), png_file)
  png_file
}

# tick labels in math mode, so minus signs and decimals are set by LaTeX
tex_num <- function(v) ifelse(is.na(v), NA,
  paste0("$", format(v, scientific = FALSE, drop0trailing = TRUE, trim = TRUE), "$"))
tex_theme <- theme(text = element_text(size = 17), plot.title = element_text(hjust = 0.5),
                   legend.position = "top", legend.key.width = unit(1.2, "cm"))

The data are \(n = 50\) draws from \(\Exp(\text{rate} = 2)\), generated as \(\Gam(1, \text{rate} = 2)\). The distribution is right-skewed, like the \(\Exp(1)\) sample of Lab 5. Its mean is \(1/2\). Its median is \(m = \ln 2 / 2 \approx 0.347\). Its density at the median is \(f(m) = 2e^{-2m} = 1\).

The law-school data come from the bootstrap package: law is the 15 schools of the lectures, and law82 is all 82 schools. The code reads them as bootstrap::law and does not attach the package with library(). The package has its own jackknife() and bootstrap() functions, and this lab writes both.

set.seed(536)
n <- 50
true_rate <- 2
data <- rgamma(n, shape = 1, rate = true_rate)

true_mean   <- 1 / true_rate        # 0.5
true_median <- log(2) / true_rate   # 0.3466
true_var    <- 1 / true_rate^2      # 0.25

# Efron's law-school data, from the bootstrap package (Efron and Tibshirani, 1993).
# law: the 15 schools of both resampling lectures, average LSAT and GPA of the entering class.
# law82: all 82 schools, the population the 15 were sampled from; its correlation is 0.760.
if (!requireNamespace("bootstrap", quietly = TRUE)) install.packages("bootstrap", repos = "https://cloud.r-project.org")
law   <- bootstrap::law
law82 <- bootstrap::law82
LSAT <- law$LSAT
GPA  <- law$GPA

print(round(c(mean = mean(data), median = median(data), sd = sd(data),
              cor_law = cor(LSAT, GPA)), 4))

   mean  median      sd cor_law 
 0.4995  0.4241  0.5227  0.7764

The Jackknife

Recall the jackknife from lecture? It leaves out one observation at a time. The notation is the lecture’s.

\(\hat\theta\): the estimate, computed from \(X_1, \dots, X_n\)
\(\hat\theta_{(-i)}\): the leave-one-out estimate, the same estimate computed without \(X_i\)
\(\hat\theta_{(\cdot)} = \frac1n\sum_{i=1}^n \hat\theta_{(-i)}\): the mean of the leave-one-out estimates

The empirical influence value and the pseudo-value are

\[\widehat{\text{IF}}_i = (n-1)\big(\hat\theta - \hat\theta_{(-i)}\big), \qquad \tilde\theta_i = \hat\theta + \widehat{\text{IF}}_i = n\hat\theta - (n-1)\hat\theta_{(-i)}.\]

Essentially, \(\widehat{\text{IF}}_i\) says how far deleting observation \(i\) moves the estimate. The factor \(n-1\) keeps it an order-one quantity. The pseudo-value adds that amount back onto \(\hat\theta\). The bias and variance estimates are

\[\widehat{\text{Bias}}_{jack} = (n-1)\big(\hat\theta_{(\cdot)} - \hat\theta\big),\]

\[\widehat{\Var}_{jack} = \frac{1}{n(n-1)}\sum_{i=1}^n \big(\widehat{\text{IF}}_i - \overline{\text{IF}}\big)^2 = \frac{n-1}{n}\sum_{i=1}^n \big(\hat\theta_{(-i)} - \hat\theta_{(\cdot)}\big)^2,\]

where \(\overline{\text{IF}} = \frac1n\sum_i \widehat{\text{IF}}_i\). The two forms of the variance are algebraically identical, since \(\widehat{\text{IF}}_i - \overline{\text{IF}} = -(n-1)(\hat\theta_{(-i)} - \hat\theta_{(\cdot)})\). Use the centered form. The uncentered \(\frac{1}{n(n-1)}\sum_i \widehat{\text{IF}}_i^2\) agrees with it only when \(\overline{\text{IF}} = 0\).

Fill in the four blanks in the code chunk below.

The empirical influence values.
The pseudo-values.
The jackknife bias.
The standard error, from the centered form.

Every later exercise calls this function.

jackknife <- function(x, stat) {
  n <- length(x)
  theta_hat <- stat(x)
  theta_i   <- vapply(seq_len(n), function(i) stat(x[-i]), numeric(1))
  theta_dot <- mean(theta_i)
  # TODO: empirical influence values and pseudo-values
  infl   <- ___________
  pseudo <- ___________
  list(estimate = theta_hat, theta_i = theta_i, theta_dot = theta_dot,
       infl = infl, pseudo = pseudo,
       # TODO: the jackknife bias and the centered standard error
       bias = ___________,
       se   = ___________)
}

Part 1: The Jackknife for Smooth Statistics

The jackknife works when the statistic is a smooth function of the data, meaning a small change in one observation moves it by a proportional amount. The mean and the variance are the clean cases.

Exercise 1.1: The Mean

Apply the jackknife to the sample mean. A \(1-\alpha\) jackknife interval uses the bias-corrected estimate and a \(t\) quantile,

\[\hat\theta - \widehat{\text{Bias}}_{jack} \pm t_{n-1,\,1-\alpha/2}\,\widehat{\text{SE}}_{jack}.\]

jk_mean <- jackknife(data, mean)

alpha <- 0.05
t_crit <- qt(1 - alpha / 2, df = n - 1)
ci_mean <- jk_mean$estimate - jk_mean$bias + c(-1, 1) * t_crit * jk_mean$se

print(round(c(estimate = jk_mean$estimate, bias = jk_mean$bias, se = jk_mean$se), 6))
print(round(c(ci_lower = ci_mean[1], ci_upper = ci_mean[2]), 4))

# theoretical SE of the mean of Exp(rate): 1 / (rate * sqrt(n))
print(round(c(se_theory = 1 / (true_rate * sqrt(n)),
              ratio = jk_mean$se / (1 / (true_rate * sqrt(n)))), 4))

Report.

The estimate, the jackknife bias, and the jackknife standard error.
The interval.
One sentence on how close the standard error is to \(1/(\text{rate}\sqrt n)\).

Derive. The jackknife bias for the mean should come out numerically zero. Use \(\hat\theta_{(-i)} = (n\bar x - x_i)/(n-1)\).

Show that \(\hat\theta_{(\cdot)} = \bar x\), so the bias is exactly zero.
Show that the influence value is \(\widehat{\text{IF}}_i = x_i - \bar x\).
Show that the pseudo-value is \(\tilde\theta_i = x_i\).
Confirm 2 and 3 with the code chunk below.

print(all.equal(jk_mean$infl, data - mean(data)))
print(all.equal(jk_mean$pseudo, data))

The influence plot puts \(\widehat{\text{IF}}_i\) against \(x_i\). If your identity is right, the points fall on the line \(\widehat{\text{IF}} = x - \bar x\), which has slope one. The code chunk below is complete. Run it as is.

infl_df <- data.frame(x = data, infl = jk_mean$infl)

p_infl <- ggplot(infl_df, aes(x, infl)) +
  geom_abline(intercept = -mean(data), slope = 1, colour = pal[["dark"]], linewidth = 1) +
  geom_point(colour = pal[["blue"]], size = 2.6, alpha = 0.85) +
  geom_hline(yintercept = 0, linetype = "dashed", colour = pal[["grey"]]) +
  scale_x_continuous(labels = tex_num) + scale_y_continuous(labels = tex_num) +
  labs(x = "$x_i$", y = r"($\widehat{\mathrm{IF}}_i$)",
       title = r"(Influence values for the mean fall on $\widehat{\mathrm{IF}}_i = x_i - \bar x$)") +
  tex_theme

knitr::include_graphics(tikz_png(p_infl, "figures/lab06/p1-mean-influence.png", width = 8, height = 5))

Coverage. Recall coverage from Lab 5? The coverage of an interval method is \(P_F(L \leq \theta \leq U)\), the probability over samples from \(F\) that its interval contains \(\theta\). One sample gives one interval, so it cannot say whether the method keeps its promise of \(0.95\). Algorithm 7 of Lab 5 can. Draw many datasets from a known \(F\), build the interval on each, and count how often it contains \(\theta\).

The code chunk below does this for the jackknife interval of the mean, with \(R = 1000\) datasets of size \(n = 50\) from \(\Exp(\text{rate} = 2)\) and \(\theta = 1/2\). It is complete. Run it as is. Exercises 3.1 and 4.1 reuse jack_ci(), coverage_study(), and draw_exp().

set.seed(536)
jack_ci <- function(x, stat, alpha = 0.05) {      # the jackknife interval, for any statistic
  jk <- jackknife(x, stat)
  jk$estimate - jk$bias + c(-1, 1) * qt(1 - alpha / 2, length(x) - 1) * jk$se
}

# Algorithm 7 of Lab 5: draw() returns one dataset, make_ci() returns one interval per row
coverage_study <- function(R, draw, make_ci, theta) {
  hit <- below <- above <- NULL
  for (r in seq_len(R)) {
    ci_r <- make_ci(draw())
    if (is.null(hit)) hit <- below <- above <- setNames(numeric(nrow(ci_r)), rownames(ci_r))
    hit   <- hit   + (ci_r[, 1] <= theta & theta <= ci_r[, 2])
    below <- below + (theta < ci_r[, 1])
    above <- above + (theta > ci_r[, 2])
  }
  data.frame(method = names(hit), coverage = hit / R, mc_se = sqrt(hit / R * (1 - hit / R) / R),
             theta_below_L = below / R, theta_above_U = above / R, row.names = NULL)
}

draw_exp <- function() rgamma(50, shape = 1, rate = true_rate)   # one dataset, n = 50
cov_mean <- coverage_study(1000, draw_exp, function(x) rbind(jackknife = jack_ci(x, mean)), true_mean)
knitr::kable(cov_mean, digits = 3)

In the table, theta_below_L is the proportion of datasets with \(\theta < L\), and theta_above_U is the proportion with \(\theta > U\). Each is nominally \(0.025\). The column mc_se is the Monte Carlo standard error of the coverage.

Report.

The coverage and its Monte Carlo standard error.
The side on which the interval misses more often.
One sentence on why this interval is the \(t\) interval of Lab 5. Then its coverage at \(n = 17\) in Lab 5 next to this one.

Exercise 1.2: The Variance

Apply the jackknife to two variance estimators. R’s var uses the divisor \(n-1\) and is unbiased. The plug-in variance uses the divisor \(n\), and its bias is \(-\sigma^2/n\). Fill in the plug-in variance in the code chunk below.

# TODO: the plug-in variance, mean squared deviation (divisor n)
plugin_var <- function(x) ___________

jk_var    <- jackknife(data, var)
jk_plugin <- jackknife(data, plugin_var)

tab_var <- rbind(
  unbiased = c(estimate = jk_var$estimate,    bias = jk_var$bias,
               corrected = jk_var$estimate - jk_var$bias,    se = jk_var$se),
  plugin   = c(estimate = jk_plugin$estimate, bias = jk_plugin$bias,
               corrected = jk_plugin$estimate - jk_plugin$bias, se = jk_plugin$se))
print(round(tab_var, 5))
print(round(c(minus_sigma2_over_n = -true_var / n,
              minus_s2_over_n = -var(data) / n), 5))

Report.

The table for both estimators.
Which estimator has zero jackknife bias, and why.
The plug-in estimator’s jackknife bias next to \(-s^2/n\) and next to its true bias \(-\sigma^2/n\).
Whether the bias-corrected plug-in variance equals the unbiased \(s^2\).

Part 2: Where the Jackknife Fails: the Median

The median is not a smooth function of the data. Delete one observation and the median jumps to one of a few fixed values. So the jackknife has almost nothing to work with.

Exercise 2.1: The Leave-One-Out Median Collapses

Jackknife the median and count the distinct leave-one-out values.

jk_med <- jackknife(data, median)
xs <- sort(data)

print(c(n_unique = length(unique(round(jk_med$theta_i, 12)))))
print(round(sort(unique(jk_med$theta_i)), 4))
print(round(c(x_25 = xs[n / 2], x_26 = xs[n / 2 + 1], full_median = median(data)), 4))
print(round(c(estimate = jk_med$estimate, bias = jk_med$bias, se = jk_med$se), 4))

Report.

The number of distinct leave-one-out medians.
The values they take.

Explain. With \(n\) even, the full median is \((x_{(n/2)} + x_{(n/2+1)})/2\).

Argue that deleting any observation below the median leaves one fixed order statistic as the new median.
Argue that deleting any observation above the median leaves the other one.
Say what changes when \(n\) is odd.

So the leave-one-out median takes only these values. That pins the jackknife standard error to a single data spacing. Fill in the closed form in the code chunk below and confirm that it matches.

# TODO: the closed-form jackknife SE for the median, n even:
# sqrt(n - 1) times half the gap between the two central order statistics
print(round(c(se_jack = jk_med$se,
              closed_form = ___________), 6))

The code chunk below draws the two bands. It is complete. Run it as is.

med_df <- data.frame(x = data, loo = jk_med$theta_i)

p_med <- ggplot(med_df, aes(x, loo)) +
  geom_hline(yintercept = median(data), linetype = "dashed", colour = pal[["grey"]]) +
  geom_point(colour = pal[["red"]], size = 2.6, alpha = 0.85) +
  scale_x_continuous(labels = tex_num) + scale_y_continuous(labels = tex_num) +
  labs(x = r"(removed value $x_i$)", y = r"($\hat\theta_{(-i)}$)",
       title = "The leave-one-out median takes two values") +
  tex_theme

knitr::include_graphics(tikz_png(p_med, "figures/lab06/p2-median-collapse.png", width = 8, height = 5))

Exercise 2.2: The Jackknife Variance for the Median Is Inconsistent

Recall from lecture that the jackknife variance for the median is inconsistent. As \(n\) grows, the estimate does not converge to the true sampling variance. We compare it with the mean, where the jackknife variance is consistent.

For large \(n\) the standard deviation of the sample median is \(1/(2 f(m)\sqrt n)\), and \(f(m) = 1\) here. The standard deviation of the sample mean is \(\sigma/\sqrt n\), with \(\sigma = 1/2\). The code chunk below does three things.

It draws \(2000\) samples at each of \(n = 25, 100, 400\).
On each sample it takes the ratio of the jackknife SE to its target, for the median and for the mean.
It reports the mean and the standard deviation of each ratio.

For a consistent estimator the ratio concentrates at one, so its standard deviation shrinks as \(n\) grows. The code chunk is complete. Run it as is.

set.seed(536)

jack_se_median <- function(x) {   # closed form, n even; loop otherwise
  x <- sort(x); m <- length(x)
  if (m %% 2 == 0) return(sqrt(m - 1) * (x[m / 2 + 1] - x[m / 2]) / 2)
  ti <- vapply(seq_len(m), function(i) median(x[-i]), numeric(1))
  sqrt((m - 1) / m * sum((ti - mean(ti))^2))
}
jack_se_mean <- function(x) sd(x) / sqrt(length(x))

ns <- c(25, 100, 400)
D  <- 2000
incons <- do.call(rbind, lapply(ns, function(nn) {
  se_med_true  <- 1 / (2 * sqrt(nn))          # 1 / (2 f(m) sqrt(n)), f(m) = 1
  se_mean_true <- (1 / true_rate) / sqrt(nn)
  r_med  <- numeric(D); r_mean <- numeric(D)
  for (d in seq_len(D)) {
    x <- rgamma(nn, shape = 1, rate = true_rate)
    r_med[d]  <- jack_se_median(x) / se_med_true
    r_mean[d] <- jack_se_mean(x)   / se_mean_true
  }
  data.frame(n = nn,
             median_ratio_mean = mean(r_med),  median_ratio_sd = sd(r_med),
             mean_ratio_mean   = mean(r_mean), mean_ratio_sd   = sd(r_mean))
}))
print(round(incons, 3))

Report.

The table.
The statistic whose ratio has a shrinking standard deviation as \(n\) grows, and the rate.
The statistic whose ratio does not settle.
One sentence on what that says about consistency.

The BCa Interval

The bootstrap is consistent for the median, where the jackknife is not. The BCa (bias-corrected and accelerated) interval is a bootstrap interval. Draw \(B\) resamples, compute \(\hat\theta^{*(b)}\) on resample \(b\), and read off two adjusted quantiles of the replicates. Recall the percentile interval from Lab 5? It is \([\hat\theta^*_{(\alpha/2)}, \hat\theta^*_{(1-\alpha/2)}]\), where \(\hat\theta^*_{(q)}\) is the \(q\)-quantile of the replicates. The percentile interval respects monotone transformations. The interval for \(\phi(\theta)\) is \(\phi\) applied to the interval for \(\theta\) whenever \(\phi\) is monotone. BCa keeps that property and moves the two levels,

\[\begin{aligned} \alpha_1 &= \Phi\!\left(\hat z_0 + \frac{\hat z_0 + z_{\alpha/2}}{1 - \hat a(\hat z_0 + z_{\alpha/2})}\right), \\[6pt] \alpha_2 &= \Phi\!\left(\hat z_0 + \frac{\hat z_0 + z_{1-\alpha/2}}{1 - \hat a(\hat z_0 + z_{1-\alpha/2})}\right), \end{aligned}\]

where \(\Phi\) is the standard normal CDF and \(z_q = \Phi^{-1}(q)\). The interval is \([\hat\theta^*_{(\alpha_1)}, \hat\theta^*_{(\alpha_2)}]\). The bias correction is

\[\hat z_0 = \Phi^{-1}\!\left(\frac{\#\{\hat\theta^{*(b)} < \hat\theta\}}{B}\right),\]

the normal quantile of the proportion of replicates below \(\hat\theta\). The acceleration comes from the jackknife,

\[\hat a = \frac{\sum_{i=1}^n (\hat\theta_{(\cdot)} - \hat\theta_{(-i)})^3}{6\big[\sum_{i=1}^n (\hat\theta_{(\cdot)} - \hat\theta_{(-i)})^2\big]^{3/2}}.\]

Essentially:

\(\hat z_0\) measures how far the bootstrap distribution sits off-center from \(\hat\theta\).
\(\hat a\) measures how fast the standard error changes with \(\theta\). It is the skewness of the influence values \(\widehat{\text{IF}}_i\) divided by \(6\sqrt n\).
If \(\hat z_0 = 0\) and \(\hat a = 0\), the two levels return to \(\alpha/2\) and \(1-\alpha/2\), and BCa is the percentile interval.
BCa is second-order accurate. Each one-sided coverage error is \(O(1/n)\), against \(O(1/\sqrt n)\) for the percentile interval.

This is Algorithm 5 of Lab 5. There the statistic was the mean, and \(\hat a\) had a closed form. Here \(\hat a\) comes from the jackknife, so the same code works for any statistic.

Part 3: The BCa Bootstrap

Exercise 3.1: BCa for the Median

Implement the bootstrap and BCa for a general statistic, and apply them to the median, where the jackknife failed. Fill in the four blanks in the code chunk below.

One bootstrap replicate: resample \(n\) observations with replacement and apply stat.
The bias correction \(\hat z_0\).
The acceleration \(\hat a\), from the leave-one-out estimates.
The two adjusted levels \(\alpha_1\) and \(\alpha_2\).

set.seed(536)
bootstrap <- function(x, stat, B) {
  n <- length(x)
  # TODO: one bootstrap replicate, resample n with replacement and apply stat
  vapply(seq_len(B), function(b) ___________, numeric(1))
}

bca_ci <- function(x, stat, B = 2000, alpha = 0.05) {
  n <- length(x)
  theta_hat  <- stat(x)
  boot_stats <- bootstrap(x, stat, B)

  # TODO: bias correction z0, from the fraction of replicates below theta_hat
  z0 <- ___________

  jk   <- jackknife(x, stat)
  U    <- jk$theta_dot - jk$theta_i
  # TODO: acceleration a from U
  a    <- ___________

  z   <- qnorm(c(alpha / 2, 1 - alpha / 2))
  # TODO: the two BCa-adjusted percentile levels
  adj <- ___________
  list(estimate = theta_hat, boot = boot_stats, z0 = z0, a = a, se = sd(boot_stats),
       bca  = unname(quantile(boot_stats, adj)),
       perc = unname(quantile(boot_stats, c(alpha / 2, 1 - alpha / 2))))
}

bca_med <- bca_ci(data, median, B = 2000)
print(round(c(estimate = bca_med$estimate, se = bca_med$se,
              z0 = bca_med$z0, a = bca_med$a), 4))
print(round(rbind(percentile = bca_med$perc, BCa = bca_med$bca), 4))

The next code chunk draws the bootstrap distribution with both intervals. It is complete. Run it as is.

boot_df <- data.frame(theta = bca_med$boot)
lines_df <- data.frame(
  x = c(bca_med$estimate, bca_med$perc, bca_med$bca),
  kind = c("estimate", "percentile", "percentile", "BCa", "BCa"))
lines_df$kind <- factor(lines_df$kind, levels = c("estimate", "percentile", "BCa"))

p_boot <- ggplot(boot_df, aes(theta)) +
  geom_histogram(bins = 45, fill = pal[["blue"]], colour = "white", alpha = 0.85) +
  geom_vline(data = lines_df, aes(xintercept = x, colour = kind, linetype = kind),
             linewidth = 1.1) +
  scale_colour_manual(values = c(estimate = pal[["dark"]], percentile = pal[["green"]],
                                 BCa = pal[["red"]])) +
  scale_linetype_manual(values = c(estimate = "solid", percentile = "dashed", BCa = "dotted")) +
  scale_x_continuous(labels = tex_num) + scale_y_continuous(labels = tex_num) +
  labs(x = r"($\hat\theta^*$)", y = "count", colour = NULL, linetype = NULL,
       title = "Bootstrap distribution of the median with two intervals") +
  tex_theme

knitr::include_graphics(tikz_png(p_boot, "figures/lab06/p3-median-boot.png", width = 9, height = 5))

Cross-check your interval against boot::boot.ci. The code chunk below is complete. Run it as is.

set.seed(536)
bt <- boot(data, function(d, i) median(d[i]), R = 2000)
print(round(boot.ci(bt, type = c("perc", "bca"))$bca[4:5], 4))

Report.

The estimate, the bootstrap SE, \(\hat z_0\), and \(\hat a\).
The percentile and BCa intervals.
One sentence on why \(\hat a = 0\) here. Look at the two leave-one-out values from Part 2.
The bootstrap SE next to the jackknife SE from Part 2.
The boot.ci interval next to yours.

Coverage. The code chunk below estimates the coverage of three intervals for the median.

The jackknife interval of Exercise 1.1, with the median as the statistic.
The percentile interval.
The BCa interval.

It uses \(R = 1000\) datasets of size \(n = 50\) from \(\Exp(\text{rate} = 2)\), with \(B = 999\) resamples each and \(\theta = \ln 2/2\). It is complete and calls your bca_ci(). Run it as is. It takes about half a minute.

set.seed(536)
three_intervals <- function(x, stat, B = 999) {   # one interval per row
  bc <- bca_ci(x, stat, B)
  rbind(jackknife = jack_ci(x, stat), percentile = bc$perc, BCa = bc$bca)
}
cov_med <- coverage_study(1000, draw_exp, function(x) three_intervals(x, median), true_median)
knitr::kable(cov_med, digits = 3)

Report.

The table.
The intervals whose coverage is within two Monte Carlo standard errors of \(0.95\).
The interval that fails, and the exercise that explains why.

Part 4: Efron’s Law-School Data

The law-school sample is the running example of both resampling lectures. It has 15 schools, each with an average LSAT score and GPA. The question is the correlation between them. The sample correlation is about 0.78. The population value across all 82 schools is about 0.76. The correlation is a smooth statistic, so the jackknife works. Its sampling distribution is skewed, so BCa has real corrections to make.

Exercise 4.1: The Jackknife and BCa for the Correlation

Jackknife the correlation. The code chunk below computes the 15 leave-one-out values and the influence values. Fill in the index of the school that moves the estimate most.

rho <- function(idx) cor(LSAT[idx], GPA[idx])
jk_law <- jackknife(seq_len(length(LSAT)), rho)

law_tab <- data.frame(school = seq_along(LSAT), LSAT = LSAT, GPA = GPA,
                      loo = round(jk_law$theta_i, 4), influence = round(jk_law$infl, 3))
print(law_tab, row.names = FALSE)
print(round(c(estimate = jk_law$estimate, bias = jk_law$bias, se = jk_law$se), 4))

# TODO: index of the school with the largest absolute influence value
top <- ___________
print(c(most_influential_school = top,
        influence = round(jk_law$infl[top], 3),
        loo_without_it = round(jk_law$theta_i[top], 4)))

Report.

The 15 leave-one-out correlations and influence values.
The jackknife SE and bias.
The most influential school.
One sentence on why removing it moves the correlation the way it does. Look at its LSAT and GPA against the others.
What the sign of its influence value says.

The next code chunk plots the influence values. It is complete. Run it as is.

infl_law <- data.frame(school = factor(seq_along(LSAT)), influence = jk_law$infl)
infl_law$hi <- seq_along(LSAT) == top

p_law_infl <- ggplot(infl_law, aes(school, influence, fill = hi)) +
  geom_col(width = 0.7) +
  geom_hline(yintercept = 0, colour = pal[["grey"]]) +
  scale_fill_manual(values = c("FALSE" = pal[["blue"]], "TRUE" = pal[["red"]]), guide = "none") +
  scale_y_continuous(labels = tex_num) +
  labs(x = "school", y = r"($\widehat{\mathrm{IF}}_i$)",
       title = "Empirical influence values for the law-school correlation") +
  tex_theme

knitr::include_graphics(tikz_png(p_law_infl, "figures/lab06/p4-influence.png", width = 9, height = 5))

Now the BCa interval. The sampling distribution of the correlation is skewed, so \(\hat z_0\) and \(\hat a\) are nonzero and BCa should differ from the percentile interval.

set.seed(536)
bca_law <- bca_ci(seq_len(length(LSAT)), rho, B = 2000)
print(round(c(z0 = bca_law$z0, a = bca_law$a, se = bca_law$se), 4))
print(round(rbind(percentile = bca_law$perc, BCa = bca_law$bca), 4))

# the acceleration is the skewness of the influence values over 6 sqrt(n); confirm it
skew <- function(v) { v <- v - mean(v); mean(v^3) / mean(v^2)^(3 / 2) }
print(round(c(a_from_bca = bca_law$a,
              a_from_skew = skew(jk_law$infl) / (6 * sqrt(length(LSAT)))), 5))

The next code chunk draws the bootstrap distribution with both intervals and the population value. It is complete. Run it as is.

law_boot_df <- data.frame(theta = bca_law$boot)
law_lines <- data.frame(
  x = c(bca_law$estimate, bca_law$perc, bca_law$bca, 0.760),
  kind = c("estimate", "percentile", "percentile", "BCa", "BCa", "population"))
law_lines$kind <- factor(law_lines$kind, levels = c("estimate", "percentile", "BCa", "population"))

p_law_boot <- ggplot(law_boot_df, aes(theta)) +
  geom_histogram(bins = 45, fill = pal[["blue"]], colour = "white", alpha = 0.85) +
  geom_vline(data = law_lines, aes(xintercept = x, colour = kind, linetype = kind),
             linewidth = 1.1) +
  scale_colour_manual(values = c(estimate = pal[["dark"]], percentile = pal[["green"]],
                                 BCa = pal[["red"]], population = pal[["grey"]])) +
  scale_linetype_manual(values = c(estimate = "solid", percentile = "dashed",
                                   BCa = "dotted", population = "longdash")) +
  scale_x_continuous(labels = tex_num) + scale_y_continuous(labels = tex_num) +
  labs(x = r"($\hat\rho^*$)", y = "count", colour = NULL, linetype = NULL,
       title = "Bootstrap distribution of the law-school correlation") +
  tex_theme

knitr::include_graphics(tikz_png(p_law_boot, "figures/lab06/p4-correlation-boot.png", width = 9, height = 5))

Cross-check against boot::boot.ci. The code chunk below is complete. Run it as is.

set.seed(536)
bt_law <- boot(cbind(LSAT, GPA), function(d, i) cor(d[i, 1], d[i, 2]), R = 2000)
print(round(boot.ci(bt_law, type = c("perc", "bca"))$bca[4:5], 4))

Report.

\(\hat z_0\) and \(\hat a\).
The check that \(\hat a\) is the skewness of the influence values over \(6\sqrt n\).
The percentile and BCa intervals.
Which side BCa widens, and why.
The boot.ci interval next to yours.
Whether each interval contains the population value \(0.76\).

Coverage. The 15 schools are a sample from the 82 law schools in law82, and the correlation over all 82 is \(\theta = 0.760\). The code chunk below treats the 82 as the population. A dataset is 15 schools drawn from the 82 with replacement, stored as their row indices. The statistic rho82() takes the indices, like rho() above. The study uses \(R = 1000\) datasets and \(B = 999\) resamples each, for the same three intervals as in Exercise 3.1. The code chunk is complete. Run it as is. It takes about half a minute.

set.seed(536)
theta82 <- cor(law82$LSAT, law82$GPA)                         # the population correlation, 0.760
rho82 <- function(idx) cor(law82$LSAT[idx], law82$GPA[idx])   # the statistic on 15 row indices
draw_law <- function() sample.int(82, 15, replace = TRUE)     # one dataset: 15 of the 82 schools

cov_cor <- coverage_study(1000, draw_law, function(idx) three_intervals(idx, rho82), theta82)
knitr::kable(cov_cor, digits = 3)

Report.

The table.
The side on which the percentile interval misses more often, and what BCa does to that miss rate.
Whether any interval reaches \(0.95\) at \(n = 15\).

Summary and Reflection

Takeaways

Leave one out. The jackknife reads bias and variance off the \(n\) leave-one-out estimates \(\hat\theta_{(-i)}\), through the influence values and pseudo-values, in one reusable function.
Smooth statistics. For the mean and the variance the jackknife gives a reliable standard error. Its bias estimate turns a biased variance estimator into an unbiased one.
The median breaks it. The leave-one-out median takes only a few distinct values. The resulting variance estimate is inconsistent, which is worse than imprecise.
BCa. The bias correction \(\hat z_0\) and the acceleration \(\hat a\) move the two percentile levels. BCa is second-order accurate and reduces to the percentile interval when both are zero.
Influence in practice. Jackknifing a correlation names the observations that move the estimate most. BCa corrects a percentile interval that is too short on one side.
Coverage decides. One sample cannot say whether an interval method keeps its promise. A coverage study can, for one statistic and one sample size at a time.

Self-Assessment Questions

Why does the jackknife work well for the mean but fail for the median?
What diagnostic immediately reveals jackknife failure?
What are the two correction factors in BCa, and what does each one do?
When would you choose the jackknife over the BCa bootstrap?
What happens to the jackknife’s number of distinct leave-one-out values as \(n\) increases for the median?
The jackknife SE for the median is the right size on average. Why does the jackknife interval still cover far less than \(0.95\)?