## Learning objectives {.unnumbered}
- Distinguish sampling *with* replacement from sampling *without* replacement, and identify which is appropriate in each setting.
- State and justify the two fundamental rules of bootstrapping: bootstrap samples must be drawn *with* replacement and must be the *same size* as the original sample.
- Explain why the bootstrap distribution and the sampling distribution share approximately the same spread even though they are centred in different places.
- Write R code with `rep_sample_n()` to generate a bootstrap distribution for the mean, the median, the standard deviation, and a proportion.
- Compute the bootstrap standard error and interpret it as an estimate of how much estimates typically vary from sample to sample.
- Identify the main failure modes of bootstrapping: small samples, unrepresentative samples, wrong bootstrap sample size, and statistics defined by extreme order values.
## Introduction {.unnumbered}
In @sec-module01, we discussed the problem of using a random sample to estimate a population parameter: the *sampling variability*. Since statistics depend on the random sample, statistics are also random, meaning that the value of a statistic varies from sample to sample. To quantify that variability, we introduced the *sampling distribution*, which describes how a statistic varies across all possible samples of a given size from a population.
Unfortunately, generally speaking, we cannot obtain the exact sampling distribution.
So, we need to approximate it. The best approximation we can obtain is by drawing
numerous samples from the population and computing the statistic of interest for each sample.
However, since this approach requires us to have access to the population,
it is not feasible in practice. So, we called it the *surreal approach*.
In practice, we have only one sample, and we want to know how much the statistic of interest would vary if we could take many samples. Let us return to the Apple quality-control team. After testing $n = 300$ screens (all now destroyed), they have one estimate of the average crack pressure: $\bar{x} \approx 1{,}008$ psi. They cannot run the test again on those screens; they are gone. Neither they can use more screens from the same shipment.
Yet the question cannot be avoided: *how reliable is their estimate?* A different random draw of $300$ screens from the same shipment would have produced a different number. The true population mean $\mu$ could be $995$ psi, or $1{,}020$ psi, or something else entirely. To answer this precisely, they would need the *sampling distribution* of $\bar{X}$, and building that requires knowing the population. They do not.
> *Can we somehow approximate the sampling distribution using only the one sample we have?*
This is the problem bootstrapping solves. The idea is deceptively simple: since the sample is our best available information about the population, we treat it as a proxy for the population and redraw from it, many times over, studying the distribution of the resulting estimates. The collection of those estimates, called the *bootstrap distribution*, closely mimics the spread of the sampling distribution. From it, we can extract a principled estimate of uncertainty without ever returning to the population.
## Sampling With and Without Replacement
Before we can understand bootstrapping, we need to nail down a fundamental choice that arises in any sampling procedure: after we select an element, do we put it back before making the next draw?
This is the distinction between sampling *with* and *without* replacement, and it turns out to be the single most important technical detail in bootstrapping.
### Sampling without replacement
To draw a simple random sample of size $n$ *without* replacement from a population of size $N$:
1. List all $N$ elements of the population.
2. Select one element uniformly at random.
3. *Remove* it from the population.
4. Repeat steps 2–3 until you have drawn exactly $n$ elements.
Because each selected element is removed before the next draw, no element can appear more than once. The draws are not strictly independent: removing one element changes the pool composition for the next draw. When $n$ is small relative to $N$ (say, $n/N < 0.05$), this dependence is negligible in practice.
*Sampling without replacement is what we do when taking a real sample from a population.* In the quality-control test, each of the $300$ screens was tested once and destroyed in the process.
### Sampling with replacement
To draw a simple random sample of size $n$ *with* replacement:
1. List all $N$ elements of the population.
2. Select one element uniformly at random.
3. *Return* it to the population.
4. Repeat steps 2–3 until you have drawn exactly $n$ elements.
Because each selected element is returned before the next draw, the same element can appear multiple times in the sample. Each draw is independent: the pool composition never changes.
*Sampling with replacement is what bootstrapping uses.* We will see exactly why shortly.
:::{.callout-important}
## Common mistake: mixing up the two settings
When sampling from a *population*, we sample *without* replacement: each individual is unique and should not be counted twice in the data. When bootstrapping from a *sample*, we sample *with* replacement: this is the mechanism that generates variability between bootstrap samples. Mixing up these two contexts is the most common bootstrapping error.
:::
#### Exercises
:::{#exr-bs-2-1}
A coin collector has a jar containing exactly $9$ Canadian dimes, each minted in a different year: $1982$, $1985$, $1989$, $1994$, $1997$, $2004$, $2010$, $2015$, and $2021$. She shakes the jar and draws $9$ dimes one at a time.
When she lays them out, she notices the $2015$ dime appears twice. Which sampling scheme did she use?
```{ojs}
//| echo: false
viewof answer_bs_2_1 = Inputs.radio(
["Sampling without replacement: she removed each dime before the next draw.",
"Sampling with replacement: she returned each dime to the jar before the next draw.",
"It is impossible to tell from this information alone."],
{label: "Sampling scheme: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_2_1
if (!is.na(answer_bs_2_1)) {
if (grepl("with replacement", answer_bs_2_1)) {
cat("Correct!\U1F389 If she had sampled without replacement, each of the 9 dimes would appear exactly once; no repeats are possible. The only way to draw the same 2015 dime twice is to return it to the jar between draws. That is sampling with replacement.")
} else if (grepl("without replacement", answer_bs_2_1)) {
cat("Not quite.\U1F4AA Sampling without replacement removes each dime from the jar after drawing it. Since the jar had exactly 9 dimes and she drew 9, sampling without replacement would give each dime exactly once. The 2015 dime appearing twice is the giveaway that dimes were returned between draws.")
} else {
cat("Oopsy, try again!\U1F4AA There is enough information here to decide. The jar had exactly 9 distinct dimes. Ask yourself: if no dime was ever returned to the jar, is it possible for the same dime to be drawn twice?")
}
}
```
:::
:::{#exr-bs-2-2}
## Explore: Sampling With vs. Without Replacement in R
In R, `slice_sample(n, replace = FALSE)` samples without replacement (the default), and `slice_sample(n, replace = TRUE)` samples with replacement.
Below is a small dataset with the crack pressures (in psi) of $10$ screens. Draw a sample of size $10$ under both schemes and compare the outputs.
```{webr}
##############################################
# A small dataset of 10 screens' crack psi #
##############################################
small_screens <- tibble(
screen_id = paste("Screen", 1:10),
crack_pressure = c(820, 975, 1050, 890, 1130, 760, 980, 840, 1210, 910)
)
# Sample WITHOUT replacement (default)
cat("--- Without replacement ---\n")
small_screens |> slice_sample(n = 10, replace = FALSE)
# Sample WITH replacement
cat("\n--- With replacement ---\n")
small_screens |> slice_sample(n = 10, replace = TRUE)
```
Re-run the cell a few times. What do you notice about the two outputs?
\noindent (a) Which sampling scheme can produce a repeated screen in the output?
```{ojs}
//| echo: false
viewof answer_bs_2_2_a = Inputs.radio(
["Without replacement: screens can be selected more than once.",
"With replacement: screens can be selected more than once.",
"Both can produce repeated screens."],
{label: "Your answer: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_2_2_a
if (!is.na(answer_bs_2_2_a)) {
if (grepl("With replacement", answer_bs_2_2_a)) {
cat("Correct!\U1F389 Without replacement, each of the 10 screens appears exactly once; the output is simply the dataset in a random order. With replacement, the same screen can be drawn multiple times, while others may not appear at all. This is the variability that makes bootstrapping work.")
} else if (grepl("without replacement", answer_bs_2_2_a)) {
cat("Not quite.\U1F4AA With 10 screens and a sample of size 10, sampling without replacement always gives you every screen exactly once, just in a different order. It is the 'with replacement' scheme that allows a single screen to appear multiple times.")
} else {
cat("Oopsy, try again!\U1F4AA Only one of the two schemes can repeat screens. Run the cell a few more times and watch for rows with identical screen IDs.")
}
}
```
\noindent (b) Re-run the with-replacement sample several times. Is the mean crack pressure always the same across runs?
```{ojs}
//| echo: false
viewof answer_bs_2_2_b = Inputs.radio(
["Yes: the mean is always the same, since the same $10$ screens are involved.",
"No: the mean changes because some screens are counted twice while others are missing.",
"It depends on the random seed."],
{label: "Your answer: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_2_2_b
if (!is.na(answer_bs_2_2_b)) {
if (grepl("mean changes", answer_bs_2_2_b)) {
cat("Exactly right.\U1F389 With replacement, each run produces a different mix of screens; some appear twice or more, others are absent. The mean crack pressure therefore varies from run to run, and that variability is the bootstrap distribution being born. This is precisely the mechanism that lets us estimate uncertainty.")
} else if (grepl("always the same", answer_bs_2_2_b)) {
cat("Not quite.\U1F4AA That would be the case if we sampled without replacement; we'd always get the same 10 screens, just reordered. With replacement, the composition of the sample changes: a screen with 1,210 psi might appear twice while the 760 psi screen is absent, or vice versa. This shifts the mean each time.")
} else {
cat("Oopsy, try again!\U1F4AA The key here is not the random seed but the replacement mechanism. Each with-replacement draw produces a genuinely different mix of screens, so the mean crack pressure does change from run to run.")
}
}
```
:::
## You Only Have One Shot
Let us return to the quality-control team. Their single test of $n = 300$ screens consumed those screens permanently. Even if testing were free and instantaneous, the number of distinct samples of size $n = 300$ from a shipment of $N = 50{,}000$ screens is
$$\binom{50{,}000}{300} \approx 10^{926}$$
a number so astronomically large that no computer could enumerate them all. Even the *surreal approach* from Chapter 1 (Sampling Distributions), simulating thousands of samples from the population, was only feasible because we were pretending to know the population. With only our single sample, neither enumeration nor simulation from the population is possible.
[The challenge, then, is this: how do we approximate the spread of the sampling distribution from a single sample, without access to the population?]{.highlight}
The answer is bootstrapping.
## Bootstrapping: The Core Idea
The key insight is this: since the sample is our best available approximation of the population, we can *treat the sample as if it were the population* and draw new samples from it.
### The algorithm
The bootstrap procedure has five steps:
1. Start with your original sample of size $n$.
2. Draw a new sample of size $n$ *with replacement* from the original sample. This is a *bootstrap sample*.
3. Compute the statistic of interest (e.g., the sample mean) from the bootstrap sample. This is the *bootstrap statistic*.
4. Repeat steps 2 and 3 many times (at least $1{,}000$ to $2{,}000$ repetitions).
5. Collect all the bootstrap statistics. Their distribution is the *bootstrap distribution*.
Each step in the procedure has a name. The sample drawn in step 2 is a **bootstrap sample**; the statistic computed from it in step 3 is a **bootstrap statistic**; and the collection of all bootstrap statistics from step 5 is the **bootstrap distribution**.
:::{#def-bootstrap-sample}
## Bootstrap Sample and Bootstrap Statistic
A *bootstrap sample* is a sample of size $n$ drawn *with replacement* from the original sample. The statistic computed from a bootstrap sample is its *bootstrap statistic*.
:::
:::{#def-bootstrap-distribution}
## Bootstrap Distribution
The distribution of a statistic computed across many bootstrap samples. It approximates the spread of the sampling distribution of that statistic.
:::
Two questions arise naturally at this point. First: does it matter *how many* bootstrap samples we draw? Second: does it matter *how large* the original sample was? These turn out to have very different answers, and the next exercise is designed to show you why.
:::{#exr-bs-explore-B-vs-n}
## Explore: More Replicates or More Data?
Two dials control a bootstrap distribution: the number of bootstrap replicates $B$ and the original sample size $n$. Use the sliders below to change each one independently and watch what happens to the shape and width of the bootstrap distribution.
```{ojs}
//| echo: false
viewof explore_n = Inputs.range([10, 100], {step: 10, value: 20,
label: "Sample size (n): "})
viewof explore_B = Inputs.range([50, 3000], {step: 50, value: 200,
label: "Bootstrap replicates (B): "})
```
```{webr}
#| autorun: true
#| echo: false
#| output: false
# 1. Generate the toy population
set.seed(99)
toy_pop <- tibble(crack_pressure = round(rlnorm(5000, meanlog = 6.9, sdlog = 0.15), 1))
# 2. Pre-generate all 10 sample size scenarios up to B = 3000 replicates
possible_n <- seq(10, 100, by = 10)
max_B <- 3000
precomputed_list <- list()
for (i in seq_along(possible_n)) {
current_n <- possible_n[i]
set.seed(current_n * 7)
toy_sample <- toy_pop$crack_pressure[sample(1:nrow(toy_pop), current_n)]
# Fast vectorized resampling
boot_means <- colMeans(matrix(
sample(toy_sample, size = current_n * max_B, replace = TRUE),
nrow = current_n,
ncol = max_B
))
precomputed_list[[i]] <- tibble(
n = current_n,
replicate = 1:max_B,
boot_mean = boot_means
)
}
precomputed_boot <- bind_rows(precomputed_list)
```
```{webr}
#| echo: false
#| input:
#| - explore_n
#| - explore_B
# 1. Filter to current slider values
filtered_boot <- precomputed_boot |>
filter(n == explore_n, replicate <= explore_B)
boot_se_ex <- sd(filtered_boot$boot_mean)
# 2. Get original sample mean and true population mean for reference lines
set.seed(explore_n * 7)
toy_sample_ref <- toy_pop$crack_pressure[sample(1:nrow(toy_pop), explore_n)]
sample_mean_val <- mean(toy_sample_ref)
pop_mean_val <- mean(toy_pop$crack_pressure)
# 3. Render the plot (instantaneous)
filtered_boot |>
ggplot(aes(x = boot_mean)) +
geom_histogram(bins = 40, fill = "steelblue", color = "white") +
geom_vline(xintercept = sample_mean_val, color = "darkorange", linetype = "dashed", linewidth = 1) +
geom_vline(xintercept = pop_mean_val, color = "red", linetype = "dotted", linewidth = 1) +
coord_cartesian(xlim = c(850, 1200)) +
labs(
title = paste0("Bootstrap distribution (n = ", explore_n, ", B = ", explore_B, ")"),
subtitle = paste0("Bootstrap SE ≈ ", round(boot_se_ex, 1), " psi (Orange: x̄, Red: μ)"),
x = "Mean Crack Pressure (psi)",
y = "Count"
) +
theme_classic() +
theme(text = element_text(size = 10))
```
*Instructions:*
1. Set $n = 20$ and $B = 200$. Note the shape and the bootstrap SE in the subtitle.
2. Keep $n = 20$. Increase $B$ to $3{,}000$. Observe what changes and what does not.
3. Reset $B$ to $200$. Now increase $n$ to $100$. What happens to the width?
4. Answer the questions below.
\noindent (a) When you increased $B$ from $200$ to $3{,}000$ (keeping $n = 20$ fixed), what happened to the bootstrap SE?
```{ojs}
//| echo: false
viewof answer_explore_a = Inputs.radio(
["It decreased substantially: more replicates produce a more accurate SE.",
"It stayed roughly the same: the bootstrap SE barely changed.",
"It increased: more replicates spread the distribution out wider."],
{label: "Effect of increasing B: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_explore_a
if (!is.na(answer_explore_a)) {
if (grepl("stayed roughly", answer_explore_a)) {
cat("Exactly right.\U1F389 More replicates makes the histogram smoother and more stable, but the bootstrap SE barely moves. The SE is determined by the original sample size n, not by how many times we resample from it. B controls the precision of our approximation to the bootstrap distribution; n controls the width of that distribution itself.")
} else if (grepl("decreased substantially", answer_explore_a)) {
cat("Not quite.\U1F4AA Increasing B does reduce some random noise in the reported SE, but not by much; the change is small. What you should see is that the bootstrap SE at B = 200 and B = 3,000 are close to each other. The distribution becomes smoother, but its width stays approximately the same.")
} else {
cat("Oopsy, try again!\U1F4AA Increasing B does not widen the distribution; it fills in more of the shape, making the histogram smoother. The width (and therefore the SE) stays approximately constant. Changing n is what actually changes the width.")
}
}
```
\noindent (b) When you increased $n$ from $20$ to $100$ (keeping $B = 200$ fixed), what happened to the bootstrap SE?
```{ojs}
//| echo: false
viewof answer_explore_b = Inputs.radio(
["It stayed roughly the same: only increasing B reduces the SE.",
"It increased: larger samples have more variability.",
"It decreased: larger samples produce more precise estimates."],
{label: "Effect of increasing n: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_explore_b
if (!is.na(answer_explore_b)) {
if (grepl("decreased", answer_explore_b)) {
cat("Correct!\U1F389 Larger samples produce narrower bootstrap distributions; the SE drops noticeably as n grows. This is the same law of diminishing returns from the sampling distribution module: to halve the SE, you need to quadruple n. More data is the only lever that genuinely reduces estimation uncertainty.")
} else if (grepl("stayed roughly", answer_explore_b)) {
cat("Not quite.\U1F4AA Look at the bootstrap SE in the subtitle as you move the n slider. Going from n = 20 to n = 100 should roughly halve the SE; the distribution narrows substantially. It is increasing B, not n, that leaves the SE approximately unchanged.")
} else {
cat("Oopsy, try again!\U1F4AA Larger samples are more precise, not less. Each additional observation gives us more information about the population, so the estimates cluster more tightly; the bootstrap distribution narrows.")
}
}
```
\noindent (c) In your own words, what does $B$ control and what does $n$ control?
:::{.reasoning-question}
$B$, the number of bootstrap replicates, controls the *smoothness and stability* of the bootstrap distribution. More replicates reduce random jitter in the reported SE, but they do not reduce the SE itself. You are just taking more measurements of the same spread.
$n$, the original sample size, controls the *actual width* of the bootstrap distribution, and therefore the true precision of the estimator. A larger sample captures more of the population's variability, so bootstrap resamples are more homogeneous and the distribution narrows. This directly mirrors the law of diminishing returns we explored in Chapter 1: to halve the standard error, we must quadruple the sample size $n$.
The practical upshot: if your bootstrap SE is too large, the only fix is more data. More replicates is just computational effort.
:::
□
:::
### The two rules
Two rules define a legitimate bootstrap procedure. Breaking either defeats the purpose.
:::{.callout-important}
## The Two Rules of Bootstrapping
*Rule 1: Draw with replacement.*
If we sampled the original sample *without* replacement, every bootstrap sample of size $n$ would be a permutation of the original; the statistic would be the same every time and the bootstrap distribution would have zero variability. Replacement is what creates the spread we want to study.
*Rule 2: Bootstrap samples must be the same size as the original sample.*
The spread of the sampling distribution depends on $n$. If bootstrap samples have the wrong size, the bootstrap distribution simulates a different experiment:
- Larger than $n$ → artificially narrow bootstrap distribution → underestimates the true SE.
- Smaller than $n$ → artificially wide bootstrap distribution → overestimates the true SE.
:::
#### Exercises
:::{#exr-bs-4-1}
\noindent (a) Bootstrap samples are drawn from which source?
```{ojs}
//| echo: false
viewof answer_bs_4_1_a = Inputs.radio(
["The population (bootstrapping requires access to the population).",
"The original sample (bootstrapping treats the sample as a stand-in for the population).",
"A new, independently collected sample."],
{label: "Bootstrap samples are drawn from: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_4_1_a
if (!is.na(answer_bs_4_1_a)) {
if (grepl("original sample", answer_bs_4_1_a)) {
cat("Exactly right.\U1F389 The whole point of bootstrapping is that we have only one sample and no access to the population. Bootstrap samples are drawn, with replacement, from the original sample, which serves as our best approximation of the population.")
} else if (grepl("population", answer_bs_4_1_a)) {
cat("Not quite.\U1F4AA If we had access to the population, we could draw many independent samples directly and build the sampling distribution without bootstrapping. Bootstrapping exists precisely because the population is unavailable.")
} else {
cat("Oopsy, try again!\U1F4AA A second independent sample would require going back to the population; that is what bootstrapping avoids. Bootstrap samples are always drawn from the one sample we already have.")
}
}
```
\noindent (b) Which of the following would violate one of the two bootstrapping rules?
```{ojs}
//| echo: false
viewof answer_bs_4_1_b = Inputs.radio(
["Drawing bootstrap samples of size n with replacement from the original sample.",
"Drawing bootstrap samples of size 2n with replacement from the original sample.",
"Drawing 5,000 bootstrap replicates instead of the usual 1,000 to 2,000."],
{label: "Which violates a rule? "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_4_1_b
if (!is.na(answer_bs_4_1_b)) {
if (grepl("size 2n", answer_bs_4_1_b)) {
cat("Correct!\U1F389 Rule 2 requires bootstrap samples to be the same size as the original sample. Samples of size 2n simulate a study in which twice as many observations were collected, so the bootstrap distribution will be artificially narrow and will underestimate the true SE.")
} else if (grepl("size n with replacement", answer_bs_4_1_b)) {
cat("Not quite.\U1F4AA That is exactly what correct bootstrapping looks like: size n, with replacement, from the original sample. Neither rule is violated.")
} else {
cat("Oopsy, try again!\U1F4AA The number of bootstrap replicates does not affect the validity of the procedure; more replicates simply give a smoother, more stable bootstrap distribution. The two rules concern the size of each individual bootstrap sample and whether it is drawn with replacement.")
}
}
```
:::
## Building a Bootstrap Distribution in R
Let us now construct the bootstrap distribution for the quality-control problem. For learning purposes, we once again use the *surreal approach*: we will pretend to have access to the full shipment of $50{,}000$ screens, the same population from Chapter 1 (Sampling Distributions). This lets us later compare the bootstrap distribution (built from our single sample alone) against the true sampling distribution (built from many population samples).
```{webr}
#| autorun: true
#| echo: false
#| output: false
set.seed(2025)
N_screens <- 50000
apple_threshold <- 750
#########################################################
# Recreate the Apple screens population from Module 1 #
#########################################################
screens_pop <- tibble(
screen_id = sprintf("Screen %05d", 1:N_screens),
crack_pressure = round(rlnorm(N_screens, meanlog = 6.9, sdlog = 0.15), 1)
)
pop_mean <- mean(screens_pop$crack_pressure)
##################################
# Draw one sample (n = 300) #
##################################
set.seed(7)
n_sample <- 300
screens_sample <- screens_pop |> slice_sample(n = n_sample)
sample_mean_cp <- mean(screens_sample$crack_pressure)
```
```{webr}
# Explore the population and our single sample
cat("Population mean (μ): ", round(pop_mean, 1), "psi\n")
cat("Sample mean (x̄): ", round(sample_mean_cp, 1), "psi\n")
cat("Sample size (n): ", nrow(screens_sample), "\n\n")
screens_sample |> head()
```
Running this code displays the population parameter $\mu$, the sample statistic $\bar{x}$, and the sample size $n$ of our single sample. Notice that the sample mean $\bar{x}$ differs slightly from $\mu$ due to sampling variability; this single sample is the only data we would actually have in practice.
The team's estimate $\bar{x}$ differs from $\mu$. That gap is not a mistake; it is sampling variability. The question bootstrapping answers is: *how large could this gap plausibly be?*
### One bootstrap sample
We start with a single bootstrap resample, using `slice_sample()` with `replace = TRUE` and a sample size equal to the original $n$.
```{webr}
set.seed(101)
###################################
# Draw one bootstrap sample #
###################################
bootstrap_sample_1 <- screens_sample |>
slice_sample(n = nrow(screens_sample), replace = TRUE)
cat("Original sample mean: ", round(sample_mean_cp, 1), "psi\n")
cat("Bootstrap sample 1 mean: ", round(mean(bootstrap_sample_1$crack_pressure), 1), "psi\n\n")
# Inspect the first rows; notice some screens appear more than once
bootstrap_sample_1 |> arrange(screen_id) |> head(15)
```
Re-run the cell a few times. Each run produces a different bootstrap mean, because different screens are drawn with replacement. Some screens appear two or three times; others are absent entirely.
### The full bootstrap distribution
To approximate the bootstrap distribution, we need many bootstrap samples. `rep_sample_n()` from the `infer` package [@infer2021] makes this easy: it draws `reps` samples of a specified size and returns a grouped data frame with one group per replicate.
```{webr}
set.seed(1)
n_boot <- 2000 # Number of bootstrap replicates
#####################################################
# Generate 2,000 bootstrap samples and compute #
# the mean crack pressure of each one #
#####################################################
(bootstrap_dist <- screens_sample |>
rep_sample_n(size = nrow(screens_sample), # same size as original
reps = n_boot,
replace = TRUE) |> # with replacement!
summarize(bootstrap_mean = mean(crack_pressure)))
```
The printed output shows the resulting data frame of $2{,}000$ bootstrap statistics, with one mean for each replicate. If you scroll through them, you will notice that the bootstrap means vary slightly from run to run, reflecting the uncertainty of our estimate.
```{webr}
###################################################
# Visualize the bootstrap distribution #
# Orange line marks the original sample mean x̄ #
###################################################
bootstrap_dist |>
ggplot(aes(x = bootstrap_mean)) +
geom_histogram(bins = 40, fill = "steelblue", color = "white") +
geom_vline(xintercept = sample_mean_cp,
color = "darkorange", linetype = "dashed", linewidth = 1) +
annotate("text",
x = sample_mean_cp + 2, y = Inf,
label = paste0("x̅ = ", round(sample_mean_cp, 1), " psi"),
color = "darkorange", hjust = 0, vjust = 1.5, size = 3.5) +
labs(
title = "Bootstrap Distribution of the Sample Mean (n = 300, B = 2,000)",
subtitle = "Orange dashed line = original sample mean x̅",
x = "Bootstrap Mean Crack Pressure (psi)",
y = "Count"
) +
theme_classic() +
theme(text = element_text(size = 10))
```
Notice where the distribution is centered: near $\bar{x}$, the original sample mean, *not* near $\mu$, the true population mean. We will return to this observation in the next section, where it becomes the central conceptual point.
:::: {.callout-note}
## Interactive Exercise: Bootstrapping the Median
::: {#exr-bs-5-1}
The quality-control team also wants to estimate the *median* crack pressure of the shipment. Modify the bootstrap code above to compute the bootstrap distribution of the sample *median* instead of the mean. Only the `summarize()` step needs to change.
```{webr}
#| exercise: bs-5-1
set.seed(1)
n_boot <- 2000
#####################################################
# Generate 2,000 bootstrap samples and compute #
# the median crack pressure of each one #
#####################################################
bootstrap_dist_median <- screens_sample |>
rep_sample_n(size = nrow(screens_sample),
reps = n_boot,
replace = TRUE) |>
summarize(bootstrap_median = ___(crack_pressure))
bootstrap_dist_median |>
ggplot(aes(x = bootstrap_median)) +
geom_histogram(bins = 40, fill = "steelblue", color = "white") +
labs(
title = "Bootstrap Distribution of the Sample Median",
x = "Bootstrap Median Crack Pressure (psi)",
y = "Count"
) +
theme_classic()
```
```{webr}
#| check: true
#| exercise: bs-5-1
if (!exists("bootstrap_dist_median")) {
list(correct = FALSE,
message = "Oopsy, try again!\U1F4AA Make sure bootstrap_dist_median is created.")
} else if (!"bootstrap_median" %in% names(bootstrap_dist_median)) {
list(correct = FALSE,
message = "Not quite.\U1F4AA The summarize() step should create a column called bootstrap_median. Replace ___ with the appropriate function name.")
} else {
med_val <- median(screens_sample$crack_pressure)
boot_center <- mean(bootstrap_dist_median$bootstrap_median)
if (abs(boot_center - med_val) < 30) {
list(correct = TRUE,
message = paste0("Correct!\U1F389 The bootstrap distribution of the median is centered near the sample median (",
round(med_val, 1), " psi). Notice that the shape may differ slightly from the bootstrap distribution of the mean."))
} else {
list(correct = FALSE,
message = "Not quite.\U1F4AA The bootstrap median values do not look right. Make sure you are using median() in the summarize() step.")
}
}
```
::: {.hint exercise="bs-5-1"}
Replace the `___` with `median`.
:::
::: {.solution exercise="bs-5-1"}
```r
set.seed(1)
n_boot <- 2000
#####################################################
# Generate 2,000 bootstrap samples and compute #
# the median crack pressure of each one #
#####################################################
bootstrap_dist_median <- screens_sample |>
rep_sample_n(size = nrow(screens_sample),
reps = n_boot,
replace = TRUE) |>
summarize(bootstrap_median = median(crack_pressure))
bootstrap_dist_median |>
ggplot(aes(x = bootstrap_median)) +
geom_histogram(bins = 40, fill = "steelblue", color = "white") +
labs(
title = "Bootstrap Distribution of the Sample Median",
x = "Bootstrap Median Crack Pressure (psi)",
y = "Count"
) +
theme_classic()
```
:::
:::
::::
## The Bootstrap Distribution and Standard Error
Now that we can build a bootstrap distribution, let us understand what it tells us and, equally important, what it does not.
### Similar spread, different center
The bootstrap distribution and the true sampling distribution are not the same thing. They differ in where they are centered. To understand how these two concepts connect, let us contrast the **Sampling Distribution** from Chapter 1 with the **Bootstrap Distribution** of Chapter 2.
| Property | Sampling Distribution (Chapter 1) | Bootstrap Distribution (Chapter 2) |
|---|---|---|
| **Underlying Concept** | How much a statistic varies if we draw many independent samples from the *population*. | How much a statistic varies if we resample many times from the *original sample*. |
| **Number of Samples** | Simulated from thousands of independent samples. | Resampled from a *single* observed sample. |
| **Ground Truth Needed?** | Yes: requires knowing the population distribution (the surreal approach). | No: only requires the single sample we already collected in real life. |
| **Center** | Centered at the true, unknown population parameter $\mu$. | Centered at the observed sample statistic $\bar{x}$. |
| **Spread** | Standard deviation is the true Standard Error ($\text{SE}$). | Standard deviation is the Bootstrap Standard Error ($\widehat{\text{SE}}$). |
The bootstrap distribution is centered at $\bar{x}$ (our observed estimate), not at the unknown $\mu$. This shift in center is the price we pay for not having the population. The key fact (the reason bootstrapping is useful at all) is that the *spread* of the two distributions is approximately equal.
Let us verify this with the surreal approach. We draw $5{,}000$ independent samples from the population to build the true sampling distribution, then place it side-by-side with the bootstrap distribution.
```{webr}
set.seed(2)
n_reps_surreal <- 5000
######################################################
# Surreal approach: true sampling distribution #
######################################################
sampling_dist <- screens_pop |>
rep_sample_n(size = n_sample, reps = n_reps_surreal) |>
summarize(xbar = mean(crack_pressure))
######################################################
# Place both distributions in one data frame #
# so we can compare them with facet_wrap #
######################################################
comparison <- bind_rows(
sampling_dist |>
rename(value = xbar) |>
mutate(distribution = paste0("Sampling distribution\n(centered at μ = ",
round(pop_mean, 1), " psi)")),
bootstrap_dist |>
rename(value = bootstrap_mean) |>
mutate(distribution = paste0("Bootstrap distribution\n(centered at x̅ = ",
round(sample_mean_cp, 1), " psi)"))
)
comparison |>
ggplot(aes(x = value, fill = distribution)) +
geom_histogram(bins = 40, color = "white") +
facet_wrap(~distribution, ncol = 2, scales = "free_x") +
scale_fill_manual(values = c("steelblue", "darkorange")) +
labs(
title = "Sampling Distribution vs. Bootstrap Distribution",
x = "Mean Crack Pressure (psi)",
y = "Count"
) +
theme_classic() +
theme(legend.position = "none", text = element_text(size = 10))
```
The two histograms have similar shapes and similar widths; they sit at different positions on the x-axis. The sampling distribution is centered at $\mu$; the bootstrap distribution is centered at $\bar{x}$. [Bootstrapping cannot move the center of the distribution toward the truth, but it accurately captures how spread out that distribution is.]{.highlight}
### The bootstrap standard error
The *standard error* (SE) of a statistic is the standard deviation of its sampling distribution, measuring the typical gap between an estimate and the true parameter. Since we cannot compute the sampling distribution directly, we estimate it from the bootstrap distribution. These numerical summaries of the bootstrap distribution are called **bootstrap standard errors**.
:::{#def-bootstrap-se}
## Bootstrap Standard Error
The standard deviation of the bootstrap distribution, used as an estimate of the true standard error:
$$\widehat{\text{SE}} = \text{sd}(\text{bootstrap distribution})$$
:::
```{webr}
##################################################
# Bootstrap SE: sd of the bootstrap distribution #
##################################################
(bootstrap_se <- bootstrap_dist |>
summarize(se = sd(bootstrap_mean)) |>
pull(se))
##################################################
# True SE: sd of the sampling distribution #
##################################################
(true_se <- sampling_dist |>
summarize(se = sd(xbar)) |>
pull(se))
cat("\nBootstrap SE: ", round(bootstrap_se, 1), "psi\n")
cat("True SE: ", round(true_se, 1), "psi\n")
cat("Difference: ", round(abs(bootstrap_se - true_se), 1), "psi\n")
```
The bootstrap SE is close to the true SE, close enough to be practically useful. The gap reflects the randomness of our particular sample; on average, across many possible original samples, the bootstrap SE is an unbiased estimate of the true SE. [What it cannot do is replace more data: the bootstrap SE is determined by $n$, not by how many times we resample.]{.highlight}
#### Exercises
:::{#exr-bs-6-1}
\noindent (a) The center of the bootstrap distribution is approximately equal to which quantity?
```{ojs}
//| echo: false
viewof answer_bs_6_1_a = Inputs.radio(
["The true population parameter μ.",
"The original sample statistic x̄.",
"Zero (bootstrap distributions are always centered at zero).",
"The bootstrap standard error."],
{label: "The bootstrap distribution is centered at: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_6_1_a
if (!is.na(answer_bs_6_1_a)) {
if (grepl("sample statistic", answer_bs_6_1_a)) {
cat("Correct!\U1F389 Every bootstrap sample is drawn from the original sample, so the bootstrap distribution inherits that sample's center. The bootstrap distribution is centered at x̅ (our observed estimate), not at the unknown μ.")
} else if (grepl("population parameter", answer_bs_6_1_a)) {
cat("Not quite.\U1F4AA That describes the sampling distribution, not the bootstrap distribution. Since bootstrap samples are drawn from the original sample (which is centered around x̅, not necessarily μ), the bootstrap distribution is centered at x̅.")
} else if (grepl("Zero", answer_bs_6_1_a)) {
cat("Oopsy, try again!\U1F4AA For crack pressure data with x̅ around 1,000 psi, the bootstrap distribution of the mean is centered near 1,000 psi, not near zero. Bootstrap distributions are centered at the original sample statistic.")
} else {
cat("Oopsy, try again!\U1F4AA The standard error measures the spread of the distribution, not its center. The center of the bootstrap distribution is the original sample statistic x̅.")
}
}
```
\noindent (b) The bootstrap distribution and the true sampling distribution share approximately the same:
```{ojs}
//| echo: false
viewof answer_bs_6_1_b = Inputs.radio(
["Center: both are centered near the true population parameter μ.",
"Spread: both have approximately the same standard deviation.",
"Center and spread: they are essentially the same distribution.",
"Neither center nor spread: they are unrelated."],
{label: "What the two distributions share: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_6_1_b
if (!is.na(answer_bs_6_1_b)) {
if (grepl("Spread", answer_bs_6_1_b)) {
cat("Exactly right.\U1F389 The bootstrap distribution is centered at x̅ (not μ), but its spread, measured by the standard deviation, closely approximates the standard deviation of the sampling distribution, which is the standard error. This shared spread is precisely what we need to quantify uncertainty.")
} else if (grepl("Center and spread", answer_bs_6_1_b)) {
cat("Close, but one part is off.\U1F4AA The spread is indeed shared. But the centers differ: the sampling distribution is centered at μ, the bootstrap distribution at x̅. These coincide only if x̅ = μ exactly, which is never the case.")
} else if (grepl("Neither", answer_bs_6_1_b)) {
cat("Oopsy, try again!\U1F4AA Looking at the side-by-side histograms above, the two distributions have very similar widths. The bootstrap distribution does not share the same center as the sampling distribution, but it does share the same spread, and that spread is exactly what the bootstrap SE estimates.")
} else {
cat("Not quite.\U1F4AA The sampling distribution is centered at μ; the bootstrap distribution is centered at x̅. The center they share is not μ. But their spreads (standard deviations) are approximately equal, and that is what makes bootstrapping useful.")
}
}
```
:::
## Bootstrapping for Any Statistic
One of the greatest strengths of bootstrapping is its generality. The *same* resampling algorithm applies to any statistic: you change only the `summarize()` step. This is in sharp contrast to formula-based methods, which require a different derivation for each statistic, and which sometimes have no closed-form SE at all.
| Statistic | `summarize()` call |
|---|---|
| Mean | `mean(crack_pressure)` |
| Median | `median(crack_pressure)` |
| Standard deviation | `sd(crack_pressure)` |
| IQR | `IQR(crack_pressure)` |
| Proportion failing | `mean(crack_pressure < 750)` |
```{webr}
set.seed(3)
n_boot <- 2000
########################################################
# Bootstrap distributions of multiple statistics #
########################################################
boot_multiple <- screens_sample |>
rep_sample_n(size = nrow(screens_sample), reps = n_boot, replace = TRUE) |>
summarize(
boot_mean = mean(crack_pressure),
boot_median = median(crack_pressure),
boot_sd = sd(crack_pressure)
)
cat("Bootstrap SE of the mean: ", round(sd(boot_multiple$boot_mean), 1), "psi\n")
cat("Bootstrap SE of the median: ", round(sd(boot_multiple$boot_median), 1), "psi\n")
cat("Bootstrap SE of the SD: ", round(sd(boot_multiple$boot_sd), 1), "psi\n")
```
The three bootstrap SEs differ: each statistic has its own sampling variability, and bootstrapping captures each correctly. For instance, the median is often less variable than the mean for right-skewed crack pressure data.
### Bootstrapping for proportions
So far we have focused on the mean, a summary for numerical variables. But the quality-control team's most pressing question is actually about a *proportion*: what fraction of screens in the shipment fails below Apple's $750$ psi threshold? Bootstrap the proportion directly.
```{webr}
set.seed(4)
#######################################################
# Bootstrap distribution of the failure proportion #
# (screens failing below the 750 psi threshold) #
#######################################################
bootstrap_dist_prop <- screens_sample |>
rep_sample_n(size = nrow(screens_sample), reps = 2000, replace = TRUE) |>
summarize(bootstrap_prop = mean(crack_pressure < apple_threshold))
(sample_prop <- mean(screens_sample$crack_pressure < apple_threshold))
(bootstrap_se_prop <- sd(bootstrap_dist_prop$bootstrap_prop))
bootstrap_dist_prop |>
ggplot(aes(x = bootstrap_prop)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
geom_vline(xintercept = sample_prop,
color = "darkorange", linetype = "dashed", linewidth = 1) +
scale_x_continuous(labels = \(x) paste0(round(x * 100, 1), "%")) +
labs(
title = "Bootstrap Distribution of p̂ (proportion failing below 750 psi)",
subtitle = paste0("Bootstrap SE ≈ ",
round(bootstrap_se_prop * 100, 1), " percentage points"),
x = "Bootstrap Proportion",
y = "Count"
) +
theme_classic() +
theme(text = element_text(size = 10))
```
The console output reports our sample proportion $\hat{p} \approx 8.7\%$ and the bootstrap standard error of approximately $1.6$ percentage points. The plot shows the distribution of all simulated bootstrap proportions, centered at the sample proportion $\hat{p}$.
#### Exercises
:::: {.callout-note}
## Practice: UBC Blundstone Boots Survey
::: {#exr-bs-7-1}
A researcher surveys a random sample of $45$ UBC students and asks whether they own a pair of Blundstone boots. The data are stored in `boots_sample`.
(a) Compute the bootstrap distribution of the sample proportion $\hat{p}$ (proportion who own Blundstones).
(b) Report the bootstrap SE.
```{webr}
#| exercise: bs-7-1
set.seed(55)
boots_sample <- tibble(
student = 1:45,
owns_blundstones = sample(c("yes", "no"), 45,
replace = TRUE, prob = c(0.6, 0.4))
)
##################################
# Your code here #
##################################
bootstrap_dist_boots <- boots_sample |>
rep_sample_n(size = ___, reps = ___, replace = ___) |>
summarize(boot_prop = mean(owns_blundstones == "___"))
bootstrap_se_boots <- ___(bootstrap_dist_boots$boot_prop)
cat("Sample proportion:", round(mean(boots_sample$owns_blundstones == "yes"), 3), "\n")
cat("Bootstrap SE: ", round(bootstrap_se_boots, 4), "\n")
```
```{webr}
#| check: true
#| exercise: bs-7-1
if (!exists("bootstrap_dist_boots") || !exists("bootstrap_se_boots")) {
list(correct = FALSE,
message = "Oopsy, try again!\U1F4AA Make sure both bootstrap_dist_boots and bootstrap_se_boots are created.")
} else if (nrow(bootstrap_dist_boots) < 1000) {
list(correct = FALSE,
message = "Not quite.\U1F4AA Use at least 1,000 bootstrap replicates (reps = 1000 or more).")
} else if (abs(bootstrap_se_boots - 0.073) > 0.04) {
list(correct = FALSE,
message = "Not quite.\U1F4AA Check: size = 45, replace = TRUE, and mean(owns_blundstones == 'yes') for the proportion.")
} else {
list(correct = TRUE,
message = paste0("Correct!\U1F389 The bootstrap SE is approximately ",
round(bootstrap_se_boots * 100, 1),
" percentage points; estimates of the proportion who own Blundstones typically vary by about that amount from sample to sample."))
}
```
::: {.hint exercise="bs-7-1"}
- `size = 45` (same as the original sample of 45 students)
- `reps = 2000` (or any number ≥ 1,000)
- `replace = TRUE`
- `mean(owns_blundstones == "yes")` computes the proportion
- `sd(...)` computes the bootstrap SE
:::
::: {.solution exercise="bs-7-1"}
```r
set.seed(55)
boots_sample <- tibble(
student = 1:45,
owns_blundstones = sample(c("yes", "no"), 45,
replace = TRUE, prob = c(0.6, 0.4))
)
#####################################################
# Generate 2,000 bootstrap samples and compute #
# the proportion of students who own Blundstones #
#####################################################
bootstrap_dist_boots <- boots_sample |>
rep_sample_n(size = 45, reps = 2000, replace = TRUE) |>
summarize(boot_prop = mean(owns_blundstones == "yes"))
bootstrap_se_boots <- sd(bootstrap_dist_boots$boot_prop)
cat("Sample proportion:", round(mean(boots_sample$owns_blundstones == "yes"), 3), "\n")
cat("Bootstrap SE: ", round(bootstrap_se_boots, 4), "\n")
```
:::
:::
::::
## When Bootstrapping Fails: Limitations
Bootstrapping is powerful, but it is not magic. Three situations exist in which it can fail badly, and knowing when to distrust it is as important as knowing how to apply it.
### Small samples
Bootstrapping approximates the sampling distribution by treating the original sample as a stand-in for the population. When $n$ is small, that stand-in is a poor approximation: the sample does not capture the diversity of the population, and the bootstrap distribution can be lumpy, discrete, and unreliable.
```{webr}
set.seed(9)
tiny_sample <- screens_pop |> slice_sample(n = 8)
#####################################################
# Generate a bootstrap distribution for a tiny #
# sample of size n = 8 #
#####################################################
boot_tiny <- tiny_sample |>
rep_sample_n(size = 8, reps = 2000, replace = TRUE) |>
summarize(boot_mean = mean(crack_pressure))
boot_tiny |>
ggplot(aes(x = boot_mean)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
labs(
title = "Bootstrap Distribution with n = 8",
subtitle = "Spiky and discrete: few distinct means possible",
x = "Bootstrap Mean Income ($)",
y = "Count"
) +
theme_classic() +
theme(text = element_text(size = 10))
```
With only $8$ distinct values to resample from, very few distinct bootstrap means are possible. The histogram is spiky, does not resemble a bell curve, and the bootstrap SE from this distribution is an unreliable estimate of the true SE. Note how this differs from Chapter 1 (Sampling Distributions): under the central limit theorem (CLT), the true sampling distribution of the mean might still look smooth and symmetric even for small sample sizes if the population is relatively symmetric, but the bootstrap distribution is completely constrained by the small number of unique observations in our sample.
*Rule of thumb*: bootstrapping becomes unreliable below roughly $n = 20$ to $30$. With very small samples, formula-based methods (when available) are more trustworthy.
### Unrepresentative (biased) samples
Bootstrapping cannot fix a biased sample. If the original $300$ screens over-represent stronger units (perhaps because weaker screens were flagged and removed before testing, or because sampling came from a single high-performing production batch), then every bootstrap sample is drawn from that biased pool. The bootstrap distribution will be shifted, and the bootstrap SE will quantify variability around the *wrong* center.
[Bootstrapping quantifies *sampling variability*. It does not correct *sampling bias*.]{.highlight}
This is not unique to bootstrapping: any method applied to a biased sample faces the same problem. But it is worth stating explicitly, because students sometimes hope that many bootstrap replicates will wash away a biased starting point. They will not.
### Wrong bootstrap sample size
The bootstrap sample size must equal $n$. Using a different size produces a distribution with the wrong spread:
| Bootstrap size vs. original | Effect on the distribution | Consequence |
|---|---|---|
| Larger than $n$ | Artificially narrower | Underestimates SE |
| Smaller than $n$ | Artificially wider | Overestimates SE |
| Same as $n$ ✓ | Correct spread | Accurate SE |
### Extreme order statistics: the maximum and minimum
This is one of the most important and counterintuitive failures of bootstrapping.
Consider estimating the maximum crack pressure in the shipment. Our sample of $n = 300$ screens has an observed maximum: the single highest crack pressure in the dataset. Every bootstrap sample is drawn from these same $300$ values, so the bootstrap maximum can *never exceed* the original sample maximum. The bootstrap distribution of the maximum is therefore structurally truncated: it cannot explore the region above the observed maximum, even though the shipment maximum almost certainly lies there.
```{webr}
set.seed(6)
n_boot <- 2000
#####################################################
# Bootstrap distributions of the maximum and mean #
# to compare their behavior and shapes #
#####################################################
# Bootstrap distribution of the MAXIMUM
boot_max <- screens_sample |>
rep_sample_n(size = nrow(screens_sample), reps = n_boot, replace = TRUE) |>
summarize(boot_stat = max(crack_pressure)) |>
mutate(statistic = "Maximum (psi)")
# Bootstrap distribution of the MEAN (for comparison)
boot_mean_cmp <- screens_sample |>
rep_sample_n(size = nrow(screens_sample), reps = n_boot, replace = TRUE) |>
summarize(boot_stat = mean(crack_pressure)) |>
mutate(statistic = "Mean (psi)")
bind_rows(boot_max, boot_mean_cmp) |>
ggplot(aes(x = boot_stat)) +
geom_histogram(bins = 40, fill = "steelblue", color = "white") +
facet_wrap(~statistic, ncol = 2, scales = "free") +
labs(
title = "Bootstrap Distribution: Maximum vs. Mean",
x = "Bootstrap Statistic (psi)",
y = "Count"
) +
theme_classic() +
theme(text = element_text(size = 10))
```
The contrast is stark. The bootstrap distribution of the maximum is heavily concentrated at the observed sample maximum: most bootstrap samples share the same maximum as the original, because the extreme observation is selected repeatedly. This distribution *severely* underestimates the true variability of the sample maximum.
The bootstrap distribution of the mean, by contrast, is smooth and well-spread, because the mean uses all $300$ values and resampling creates genuine variability around that average.
:::{.callout-caution}
## Easy to confuse: smooth statistics vs. extreme-value statistics
Bootstrapping is reliable for statistics that are *smooth* functions of the data: mean, median, standard deviation, proportions. These statistics incorporate all observations and respond well to resampling.
It is unreliable for statistics defined by *one or a few extreme observations*: the maximum, the minimum, and extreme quantiles (such as the 99th percentile with small $n$). The bootstrap distribution is structurally truncated for these, and the bootstrap SE will be a severe underestimate. This failure persists regardless of how large $n$ is; it is a fundamental structural limitation, not a small-$n$ problem.
:::
#### Exercises
:::{#exr-bs-9-1}
\noindent (a) A researcher draws 3,000 bootstrap replicates from a sample of size $n = 8$. Nearly all of the bootstrap means land in only a handful of distinct values, making the histogram look like a discrete bar chart. What is the most likely cause?
```{ojs}
//| echo: false
viewof answer_bs_9_1_a = Inputs.radio(
["She used replace = FALSE instead of replace = TRUE.",
"She used too many bootstrap replicates: 3,000 is excessive for n = 8.",
"The original sample is too small: with only 8 distinct values, very few distinct bootstrap means are possible.",
"She should have used a larger bootstrap sample size than n = 8."],
{label: "Most likely cause: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_9_1_a
if (!is.na(answer_bs_9_1_a)) {
if (grepl("too small", answer_bs_9_1_a)) {
cat("Correct!\U1F389 With only 8 observations, the number of distinct bootstrap samples is severely limited, and so is the number of distinct bootstrap means. The distribution becomes discretized and lumpy, and the bootstrap SE is unreliable. Below roughly n = 20-30, formula-based methods are more trustworthy when available.")
} else if (grepl("replace = FALSE", answer_bs_9_1_a)) {
cat("Sampling without replacement would make things worse, but it is not the primary explanation here.\U1F4AA Without replacement from n = 8, every bootstrap 'sample' is an identical permutation of the original; the mean would be exactly the same every time, producing a bootstrap distribution with zero spread. The description says a 'handful of distinct values,' which is more consistent with the small-n bootstrapping problem.")
} else if (grepl("too many", answer_bs_9_1_a)) {
cat("Not quite.\U1F4AA More bootstrap replicates always helps; they produce a smoother, more stable bootstrap distribution. The number of replicates does not cause the discretization problem. The culprit is the tiny sample size: only 8 distinct values to resample from.")
} else {
cat("Oopsy, try again!\U1F4AA Using a larger bootstrap sample size would violate Rule 2 and produce an artificially narrow distribution, but it would not make the distribution discrete. The discreteness comes from having only 8 original observations to resample from.")
}
}
```
\noindent (b) A colleague computes bootstrap samples of size $n = 200$ from your original sample of size $n = 120$. What effect will this have on the bootstrap SE?
```{ojs}
//| echo: false
viewof answer_bs_9_1_b = Inputs.radio(
["The bootstrap SE will be larger than the true SE: it overestimates variability.",
"The bootstrap SE will be smaller than the true SE: it underestimates variability.",
"The bootstrap SE will be unaffected: only the number of replicates matters.",
"The bootstrap SE will be exactly zero."],
{label: "Effect on bootstrap SE: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_9_1_b
if (!is.na(answer_bs_9_1_b)) {
if (grepl("smaller", answer_bs_9_1_b)) {
cat("Correct!\U1F389 Bootstrap samples of size 200 simulate a study in which the researcher collected 200 observations, not 120. Larger samples have smaller standard errors, so the bootstrap distribution will be artificially narrow; the bootstrap SE will underestimate the true SE of an estimator based on n = 120 observations.")
} else if (grepl("larger", answer_bs_9_1_b)) {
cat("The direction is reversed.\U1F4AA Larger bootstrap samples produce narrower bootstrap distributions and smaller SEs. Smaller bootstrap samples are what overestimate the SE. Here the bootstrap size (200) exceeds the original (120), so the SE will be underestimated, not overestimated.")
} else if (grepl("unaffected", answer_bs_9_1_b)) {
cat("Oopsy, try again!\U1F4AA The bootstrap sample size directly controls the spread of the bootstrap distribution, because it determines which experiment is being simulated. Bootstrap samples of size 200 simulate collecting 200 observations, producing an artificially narrow distribution.")
} else {
cat("Oopsy, try again!\U1F4AA The bootstrap SE will not be zero; different bootstrap samples of size 200 will still yield different means. But those means will cluster more tightly than they should, because larger samples are more precise.")
}
}
```
\noindent (c) Why does the bootstrap distribution of the *maximum* fail to estimate the true SE accurately, even with a large sample?
```{ojs}
//| echo: false
viewof answer_bs_9_1_c = Inputs.radio(
["Bootstrap resamples can never exceed the original sample maximum, so the bootstrap distribution is structurally truncated.",
"The maximum is too computationally expensive to bootstrap with large samples.",
"The maximum is not a valid summary statistic and cannot be bootstrapped.",
"The maximum always equals the true population maximum, so its bootstrap distribution is degenerate."],
{label: "Why the bootstrap fails for the maximum: "}
)
```
```{webr}
#| edit: false
#| echo: false
#| output: true
#| input:
#| - answer_bs_9_1_c
if (!is.na(answer_bs_9_1_c)) {
if (grepl("structurally truncated", answer_bs_9_1_c)) {
cat("Exactly right.\U1F389 Every bootstrap sample is drawn from the same n values in the original sample, so the bootstrap maximum can never exceed the observed sample maximum. The region above the sample maximum (where the true population maximum almost certainly lies) is completely invisible to the bootstrap. This is a fundamental structural limitation that persists regardless of how large n is or how many replicates you use.")
} else if (grepl("computationally expensive", answer_bs_9_1_c)) {
cat("Not quite.\U1F4AA Computing max() is extremely fast; speed is not the issue. The problem is structural: the bootstrap distribution is capped at the observed sample maximum, so it cannot capture the true variability of the maximum as an estimator.")
} else if (grepl("not a valid summary statistic", answer_bs_9_1_c)) {
cat("Not quite.\U1F4AA The maximum is a perfectly valid summary statistic; you can compute it from any sample. The issue is not what it is, but how bootstrapping behaves for statistics that depend entirely on a single extreme observation.")
} else {
cat("Oopsy, try again!\U1F4AA The sample maximum is almost certainly smaller than the population maximum; we are very unlikely to have sampled the single strongest screen in a shipment of 50,000. Bootstrapping cannot recover information about values above the observed maximum, because it only resamples from what we observed.")
}
}
```
:::
## Take-home points
- **Sampling without replacement** removes each element after drawing it (no repeats possible). **Sampling with replacement** returns each element before the next draw (repeats possible). When taking a real sample from a population, use without replacement. Bootstrapping uses *with replacement*, because that is what creates variability between bootstrap samples.
- **The two rules of bootstrapping** require us to draw each bootstrap sample with replacement and use the same sample size $n$ as the original. Breaking either rule produces a bootstrap distribution with the wrong spread.
- **The bootstrap distribution** (@def-bootstrap-distribution) is centered at the original sample statistic $\bar{x}$, not at the true population parameter $\mu$. The two distributions share approximately the same *spread*, and that shared spread is what makes bootstrapping useful for estimating uncertainty.
- **The bootstrap standard error** (@def-bootstrap-se) is the standard deviation of the bootstrap distribution:
$$\widehat{\text{SE}} = \text{sd}(\text{bootstrap distribution})$$
It estimates the true standard error ($\text{SE}$) of the sampling distribution. It is bounded by $n$: more replicates $B$ stabilize our estimate of the bootstrap SE, but more data $n$ is the only way to actually reduce the true uncertainty. Key R code:
```r
screens_sample |>
rep_sample_n(size = nrow(screens_sample), reps = 2000, replace = TRUE) |>
summarize(boot_stat = mean(crack_pressure)) # or median(), sd(), etc.
```
- **Generality of bootstrapping** means the same resampling procedure applies to any statistic (mean, median, standard deviation, proportion, IQR, and more). Simply change the `summarize()` step.
- **Four failure modes to know** when bootstrapping becomes unreliable:
- *Small samples* ($n \lesssim 20$–$30$): too few distinct values to resample from; the bootstrap distribution is lumpy and the SE unreliable.
- *Biased samples*: bootstrapping quantifies sampling variability, not sampling bias. A biased original sample produces a biased bootstrap distribution, regardless of how many replicates you draw.
- *Wrong bootstrap sample size*: bootstrap samples larger than $n$ underestimate the SE; smaller ones overestimate it.
- *Extreme order statistics* (maximum, minimum, extreme percentiles): the bootstrap distribution is structurally capped at the observed extremes and cannot estimate the true SE. This failure persists at any sample size.