Let us say you are a PhD candidate. You have spent months cleaning your data, fitting your model, and preparing the results. You show the results to your supervisors. One of them looks at your confidence interval and asks:
Did you bootstrap this?
Another one adds:
Carry out a bootstrap and see whether the conclusion remains the same.
You write the comment down. But perhaps you are not completely sure what they have asked you to do. Is bootstrap another statistical test? Is it a way of correcting a bad model? And why does a serious statistical method have the name of something attached to a boot?
The same can happen in a management meeting. A statistician says that the confidence intervals were bootstrapped. If you are not a statistician, you may wonder what this means and whether it changes the decision.
In this article I will use a simple example of children’s weight-for-height z-scores to explain bootstrap. We will first do an ordinary test, then carry out the same test using bootstrap.
Baron Münchhausen was a fictional nobleman famous for telling impossible stories. In one of the stories, he and his horse became trapped in a swamp. He claimed that he escaped by pulling himself and the horse out using his own hair.
The expression “pull yourself up by your own bootstraps” came later. It means doing something that appears impossible without help from outside. In the original story the Baron used his hair, not his bootstraps.
Statistics has a similar problem. We want to know what would happen if we repeated our study many times and collected a new sample each time. Usually, we do not have the money or time to repeat the whole study. We have only one sample.
In bootstrap we use the one sample we have to create many new samples. Bradley Efron formally introduced the statistical method in his 1979 paper, “Bootstrap Methods: Another Look at the Jackknife”. The stories of Baron Münchhausen can also be read through Project Gutenberg.
In simple terms:
A bootstrap repeatedly resamples the data we observed, recalculates the result, and uses the variation in those recalculated results to show how uncertain our original result is.
Let us see what this means using data.
Suppose we collected weight and height measurements from children in a community nutrition survey. We then converted each child’s measurement into a weight-for-height z-score, which we will call WHZ.
The WHO child growth standards give us a reference for interpreting the values:
The data below are simulated. This allows us to learn the method without sharing information about real children.
n_children <- 80
whz <- rnorm(
n = n_children,
mean = -0.65,
sd = 1.05
)
whz <- pmax(pmin(whz, 2.4), -3.8)
children <- data.table(
child_id = seq_len(n_children),
whz = whz
)
children[
,
wasted := whz < wasting_cutoff
]
children[1:10] |>
kable(
digits = 2,
col.names = c("Child", "WHZ", "WHZ below -2")
)| Child | WHZ | WHZ below -2 |
|---|---|---|
| 1 | -0.10 | FALSE |
| 2 | -1.78 | FALSE |
| 3 | -0.50 | FALSE |
| 4 | -0.74 | FALSE |
| 5 | -1.35 | FALSE |
| 6 | -3.29 | TRUE |
| 7 | -1.42 | FALSE |
| 8 | -1.72 | FALSE |
| 9 | -0.53 | FALSE |
| 10 | -1.15 | FALSE |
Before running a test, let us look at the values.
The sample mean is -0.74. It is below zero, so the children in this sample are lighter for their height than the WHO reference median on average.
But a sample mean is not a population mean. A different sample of 80 children would give us a different answer. Our first question is:
Is the population mean WHZ below 0, or could the negative sample mean simply be sampling variation?
This can be written as:
\[ H_0: \mu = 0 \qquad \text{against} \qquad H_1: \mu < 0. \]
The null hypothesis says that the population mean WHZ is at the reference median. The alternative says it is below the reference median.
This is not yet a test of how many children are wasted. A mean below 0 and an individual WHZ below -2 are different things. We will return to that distinction.
The familiar approach is a one-sample t-test.
classical_test <- t.test(
children$whz,
mu = null_mean,
alternative = "less"
)
classical_results <- data.table(
`Sample mean` = mean(children$whz),
`t statistic` = unname(classical_test$statistic),
`p-value` = classical_test$p.value
)
kable(classical_results, digits = 4)| Sample mean | t statistic | p-value |
|---|---|---|
| -0.7397 | -6.6449 | 0 |
The p-value is 1.77^{-9}. If the population mean were really 0, a sample result this far below 0 would be unusual under the assumptions of the t-test.
The t-test compares our observed mean with what we would expect if the population mean were zero. It does this using a t-distribution. We now want to see whether bootstrap gives us the same conclusion.
We start with one bootstrap sample before repeating the process 5,000 times.
For a small demonstration, take the first 10 WHZ values.
small_sample <- children$whz[1:10]
set.seed(27)
one_resample <- sample(
small_sample,
size = length(small_sample),
replace = TRUE
)
resample_example <- data.table(
position = seq_along(small_sample),
observed_whz = small_sample,
resampled_whz = one_resample
)
kable(
resample_example,
digits = 2,
col.names = c("Position", "Observed WHZ", "One bootstrap resample")
)| Position | Observed WHZ | One bootstrap resample |
|---|---|---|
| 1 | -0.10 | -1.35 |
| 2 | -1.78 | -1.78 |
| 3 | -0.50 | -3.29 |
| 4 | -0.74 | -0.53 |
| 5 | -1.35 | -1.72 |
| 6 | -3.29 | -0.50 |
| 7 | -1.42 | -0.10 |
| 8 | -1.72 | -0.53 |
| 9 | -0.53 | -0.10 |
| 10 | -1.15 | -0.53 |
We sampled 10 values from the original 10 values with replacement. After drawing a value, we put it back before drawing again. This is why one value can appear several times while another value may not appear at all.
The original mean of these 10 values is -1.26. The mean of this particular resample is -1.04.
One resample does not tell us much. The next step is to repeat the process many times.
We now resample all 80 children 5,000 times. Each time, we calculate the mean WHZ.
set.seed(310)
bootstrap_means <- numeric(n_boot)
for (b in seq_len(n_boot)) {
whz_star <- sample(
children$whz,
size = n_children,
replace = TRUE
)
bootstrap_means[b] <- mean(whz_star)
}
mean_bootstrap_se <- sd(bootstrap_means)
mean_bootstrap_ci <- percentile_interval(bootstrap_means)
This graph is called the bootstrap distribution of the mean. It shows how the mean changes when the observed sample is resampled again and again.
The standard deviation of these 5,000 means is 0.109. We call this the bootstrap standard error. It is an estimate of the typical sampling variation in the mean.
The middle 95% of the bootstrap means run from -0.95 to -0.53. The interval is entirely below zero. In practical language, the data suggest that the population represented by this sample has a mean WHZ below the WHO reference median.
Notice what we have done so far. We resampled the observations without changing them first. This is useful for estimating the uncertainty around the observed mean. It is not yet a valid bootstrap hypothesis test.
Our observed data have a negative mean. When we resample them as they are, the bootstrap samples will also tend to have a negative mean. That is exactly what we want when building a confidence interval around the observed result.
For hypothesis testing we need to ask a different question:
What kind of test statistics would we observe if the null hypothesis were true?
The null hypothesis says that the population mean is 0. We therefore create a world where this is true by shifting every observed value by the same amount:
\[ x_i^{0} = x_i - \bar{x} + \mu_0, \]
where \(\bar{x}\) is the observed mean and \(\mu_0\) is the value specified by the null hypothesis. Here, \(\mu_0 = 0\).
This transformation keeps the shape and spread of the observed values, but moves their mean to zero.
whz_under_null <- children$whz -
mean(children$whz) +
null_mean
null_world_summary <- data.table(
Dataset = c("Observed data", "Recentered null data"),
Mean = c(mean(children$whz), mean(whz_under_null)),
`Standard deviation` = c(sd(children$whz), sd(whz_under_null))
)
kable(null_world_summary, digits = 3)| Dataset | Mean | Standard deviation |
|---|---|---|
| Observed data | -0.74 | 0.996 |
| Recentered null data | 0.00 | 0.996 |
The centre has moved to zero, but the standard deviation has not changed. We can now resample from these recentered values to see what results the null world produces.
We will use the same statistic as the ordinary t-test:
\[ T_{\text{obs}} = \frac{\bar{x} - \mu_0}{s/\sqrt{n}}. \]
For every null bootstrap sample, we calculate:
\[ T^{*(b)} = \frac{\bar{x}^{*(b)} - \mu_0}{s^{*(b)}/\sqrt{n}}. \]
This is called a studentized statistic because the difference in means is scaled by its estimated standard error. A large negative value provides evidence for our alternative hypothesis that \(\mu < 0\).
observed_t <- one_sample_t(
children$whz,
mu = null_mean
)
set.seed(711)
null_bootstrap_t <- numeric(n_boot)
for (b in seq_len(n_boot)) {
whz_star <- sample(
whz_under_null,
size = n_children,
replace = TRUE
)
null_bootstrap_t[b] <- one_sample_t(
whz_star,
mu = null_mean
)
}
bootstrap_p_value <- monte_carlo_p(
simulated_statistics = null_bootstrap_t,
observed_statistic = observed_t,
tail = "left"
)
The observed t-statistic is -6.64. The bootstrap p-value is 2^{-4}.
The p-value is the proportion of bootstrap samples generated under the null hypothesis that produced a t-statistic at least as negative as the one we observed. We use the small Monte Carlo correction:
\[ \widehat{p} = \frac{ 1 + \#\left\{T^{*(b)} \leq T_{\text{obs}}\right\} }{ B + 1 }. \]
The 1 in the numerator and denominator prevents us from reporting a p-value of exactly zero merely because a finite simulation did not produce a more extreme result.
test_comparison <- data.table(
Method = c(
"Classical one-sample t-test",
"Nonparametric bootstrap test"
),
`p-value` = c(
classical_test$p.value,
bootstrap_p_value
)
)
kable(test_comparison, digits = 4)| Method | p-value |
|---|---|
| Classical one-sample t-test | 0e+00 |
| Nonparametric bootstrap test | 2e-04 |
The two p-values do not have to be identical. The classical test uses the t-distribution as its null reference. The bootstrap test estimates the null distribution by resampling the recentered data.
What matters here is the substantive conclusion. Both methods indicate that the negative sample mean would be difficult to explain if the population mean WHZ were really zero.
No. This is an important distinction.
Our mean analysis asked whether children are below the WHO reference median on average. Wasting is an individual classification. According to the WHO definition, a child is wasted when WHZ is below -2.
If a manager asks, “What percentage of children are wasted?”, the statistic we need is:
\[ \widehat{\pi} = \frac{ \text{number of children with WHZ below } -2 }{ \text{number of children measured} }. \]
We can bootstrap this proportion in the same way we bootstrapped the mean: resample children, calculate the proportion wasted, and repeat.
observed_wasting_prevalence <- mean(
children$whz < wasting_cutoff
)
set.seed(912)
bootstrap_wasting_prevalence <- numeric(n_boot)
for (b in seq_len(n_boot)) {
whz_star <- sample(
children$whz,
size = n_children,
replace = TRUE
)
bootstrap_wasting_prevalence[b] <- mean(
whz_star < wasting_cutoff
)
}
wasting_prevalence_ci <- percentile_interval(
bootstrap_wasting_prevalence
)The estimated prevalence of wasting is 6.2%. The bootstrap interval runs from 1.2% to 12.5%.
For a manager, the interval is often more useful than the point estimate alone. Another sample from the same population would not necessarily produce exactly 6.2%. The interval shows a plausible range of sampling uncertainty.
There is an additional warning. If the children came from a complex survey with villages as clusters, strata, or unequal sampling weights, resampling individual children as though they were independently selected may be wrong. The bootstrap must respect the way the sample was collected.
Let us now write the same idea using statistical notation. Suppose the population follows an unknown distribution \(F\). The quantity we want to estimate is:
\[ \theta = t(F). \]
The function \(t(F)\) can represent a population mean, median, prevalence, or regression coefficient. We observe a sample:
\[ X_1, X_2, \ldots, X_n \sim F, \]
but we do not know \(F\). Bootstrap replaces \(F\) with the empirical distribution of the values in our sample:
\[ \widehat{F}_n(x) = \frac{1}{n} \sum_{i=1}^{n} \mathbf{1}(X_i \leq x). \]
This gives each observed value a probability of \(1/n\). Sampling from \(\widehat{F}_n\) is therefore the same as sampling the observed values with replacement.
For bootstrap sample \(b\), we draw:
\[ X_1^{*(b)}, \ldots, X_n^{*(b)} \overset{\text{i.i.d.}}{\sim} \widehat{F}_n, \]
and calculate the statistic again:
\[ \widehat{\theta}^{*(b)} = t\left(\widehat{F}_n^{*(b)}\right). \]
The collection of values \(\widehat{\theta}^{*(1)}, \ldots, \widehat{\theta}^{*(B)}\) is the bootstrap distribution. Its spread tells us how much the statistic changes from one bootstrap sample to another.
For a percentile confidence interval we take the lower and upper percentiles:
\[ \left[ \widehat{\theta}^{*}_{\alpha/2}, \widehat{\theta}^{*}_{1-\alpha/2} \right]. \]
For a 95% interval, these are the 2.5th and 97.5th percentiles. This is the method we used for the mean WHZ and wasting prevalence. Other methods, such as the bootstrap-t and bias-corrected and accelerated intervals, may be more appropriate when the bootstrap distribution is highly skewed or biased.
So far, we have used the nonparametric bootstrap. It resamples the observed values. There is another common approach.
bootstrap_methods <- data.table(
Method = c("Nonparametric bootstrap", "Parametric bootstrap"),
`What creates a new dataset?` = c(
"Resampling the observed data",
"Simulation from a fitted probability model"
),
`Main assumption` = c(
"The sample represents the population and the resampling scheme matches the design",
"The fitted model reasonably describes how the data arise"
),
`Question to ask` = c(
"What exactly was resampled?",
"What model generated the simulated data?"
)
)
kable(bootstrap_methods)| Method | What creates a new dataset? | Main assumption | Question to ask |
|---|---|---|---|
| Nonparametric bootstrap | Resampling the observed data | The sample represents the population and the resampling scheme matches the design | What exactly was resampled? |
| Parametric bootstrap | Simulation from a fitted probability model | The fitted model reasonably describes how the data arise | What model generated the simulated data? |
For our WHZ data, a parametric bootstrap could assume a Normal distribution, estimate its mean and standard deviation, and simulate new values from that fitted distribution. For hypothesis testing, we would simulate under the null:
\[ X_i^* \sim \mathcal{N}\left(\mu_0, \widehat{\sigma}^2\right). \]
This is convenient, but it adds a stronger assumption: that the Normal model is a good description of the WHZ values.
So far we have tested a sample mean. We can use the same bootstrap idea when the hypothesis concerns a regression coefficient.
Let us use the crabs dataset in the MASS package. It contains measurements from 200 crabs of two species. We will use sp, the crab species; CW, the carapace width; and sex, whether the crab is female or male.
Suppose our question is:
After accounting for carapace width, does sex help us distinguish the two crab species?
We fit a logistic regression because the outcome has two categories. The full model is:
\[ \log\left( \frac{P(\text{orange species})} {1 - P(\text{orange species})} \right) = \beta_0 + \beta_1\operatorname{CW} + \beta_2\operatorname{sex}. \]
The coefficient we are interested in is \(\beta_2\):
\[ H_0: \beta_2 = 0 \qquad \text{against} \qquad H_1: \beta_2 \neq 0. \]
In ordinary language, the null hypothesis says that once width is in the model, adding sex does not improve the model.
We first fit two models. The reduced model represents the null hypothesis and contains width only. The full model contains width and sex.
data(crabs, package = "MASS")
crab_data <- as.data.table(crabs)
crab_data[, orange_species := as.integer(sp == "O")]
crab_reduced <- glm(
orange_species ~ CW,
family = binomial(link = "logit"),
data = crab_data
)
crab_full <- glm(
orange_species ~ CW + sex,
family = binomial(link = "logit"),
data = crab_data
)
crab_lrt <- anova(
crab_reduced,
crab_full,
test = "Chisq"
)
observed_crab_lrt <- crab_lrt$Deviance[2]
classical_crab_p <- crab_lrt$`Pr(>Chi)`[2]
crab_model_result <- data.table(
Result = c(
"Observed likelihood-ratio statistic",
"Classical p-value"
),
Value = c(observed_crab_lrt, classical_crab_p)
)
kable(crab_model_result, digits = 4)| Result | Value |
|---|---|
| Observed likelihood-ratio statistic | 0.0529 |
| Classical p-value | 0.8181 |
The likelihood-ratio statistic measures how much fit we gain by adding sex:
\[ T_{\text{obs}} = 2\left\{ \ell(\widehat{\beta}_{\text{full}}) - \ell(\widehat{\beta}_{\text{reduced}}) \right\}. \]
A large value means that the full model improves substantially on the reduced model. The ordinary test compares this statistic with a chi-square reference distribution.
For the WHZ test, we recentered the values so their mean was zero. Here the reduced model already represents the null hypothesis. It says that species probability depends on width but not on sex.
We therefore take the probabilities fitted by the reduced model and use them to simulate a new species outcome for every crab. We then fit both models to the simulated data and recalculate the likelihood-ratio statistic.
set.seed(1503)
crab_null_probability <- fitted(crab_reduced)
crab_bootstrap_lrt <- numeric(n_boot)
for (b in seq_len(n_boot)) {
crab_star <- copy(crab_data)
crab_star[
,
orange_species := rbinom(
n = .N,
size = 1,
prob = crab_null_probability
)
]
reduced_star <- glm(
orange_species ~ CW,
family = binomial(link = "logit"),
data = crab_star
)
full_star <- glm(
orange_species ~ CW + sex,
family = binomial(link = "logit"),
data = crab_star
)
crab_bootstrap_lrt[b] <- as.numeric(
2 * (logLik(full_star) - logLik(reduced_star))
)
}
bootstrap_crab_p <- monte_carlo_p(
simulated_statistics = crab_bootstrap_lrt,
observed_statistic = observed_crab_lrt,
tail = "right"
)
crab_test_comparison <- data.table(
Method = c(
"Classical chi-square test",
"Parametric bootstrap test"
),
`p-value` = c(classical_crab_p, bootstrap_crab_p)
)
kable(crab_test_comparison, digits = 4)| Method | p-value |
|---|---|
| Classical chi-square test | 0.8181 |
| Parametric bootstrap test | 0.8168 |
The observed likelihood-ratio statistic is 0.0529. The classical p-value is 0.8181 and the bootstrap p-value is 0.8168. The two p-values are almost the same and both are large. We therefore do not have evidence that adding sex improves the model after accounting for carapace width.
This follows the same procedure as the WHZ test. We state the null hypothesis, generate data under that hypothesis, recalculate the statistic, and find out how often the simulated statistic is at least as large as the observed one.
This is a parametric bootstrap because the new species outcomes are generated from the reduced logistic model. The result therefore depends on that model being a reasonable description of the data under the null hypothesis.
Bootstrap is useful, but it is not magic.
Resampling a problem 5,000 times gives you 5,000 versions of the same problem.
For independent observations, resampling rows may be reasonable. For a clustered survey, we may need to resample clusters. For longitudinal data, we may need to resample people or blocks of time. For a regression, we must decide whether to resample rows, residuals, clusters, or simulate from a fitted model.
The difficult part of bootstrap is often not writing the loop. It is deciding what represents a valid repeat of the study.
If you are the PhD candidate, ask:
If you are the manager listening to a presentation, ask:
The next time a supervisor or statistician asks whether you bootstrapped the result, you do not need to hear it as an incantation. They are asking how stable your conclusion is when the observed data are made to play the role of many possible samples.
Your next question should be:
What did we resample, and did we resample it in a way that respects how the data were collected?