Fit the model a thousand times on a thousand might-have-been datasets.
Why is drawing rows with replacement a legitimate way to manufacture new datasets? The JonStats course justifies it with a short Socratic exchange, reproduced here because there is no improving on it:
“Every observation in our dataset is equally likely.”
Why is this?
“Because each specific observation in our dataset has been observed the same number of times.”
Why do you say that?
“Because each observation in the dataset has been observed exactly one time, and 1 = 1!”
And so a resampled dataset in which every row was drawn with equal probability, with replacement, to the same size $n$, is as equally likely and as equally precise as the dataset you actually observed — a dataset that might have happened. Make a thousand of them, fit your model to each, and the spread of the thousand answers estimates the sampling distribution.
Start where a formula already exists, so the loop can be checked. The curvature route priced the heart-rate slope at $\hat\beta_1 = -1.271$ with SE $0.2915$ and Wald interval $(-1.84, -0.70)$. Now bootstrap it: each click refits the line to resampled rows and adds a faint orange line; the histogram collects the slopes.
The handshake. Run the full thousand: SE* $= 0.2964$ against the curvature route's $0.2915$; interval $(-1.840, -0.668)$ against $(-1.843, -0.700)$. Two utterly different arguments — second derivatives of a log-likelihood versus brute re-fitting — arriving at the same place. This is the well-behaved case; when the two disagree, believe the disagreement and investigate (Three Routes stages exactly that).
Now the headline act. Below are 80 resting blood-pressure readings — right-skewed,
as blood pressure is (Tutorial 5 models this shape with a
Gamma GLM). For skewed data the
natural summaries are quantiles — the median, the 10th percentile, the 90th
— and none of them has a Std. Error column anywhere in classical
output. The loop does not care. Pick a quantile, run the resamples, read off its
interval.
Watch the tails. The 10th percentile's interval is about 6 mm Hg wide; the 90th's is over 19 — three times the uncertainty, because the sparse upper tail of a skewed sample pins its high quantiles far more loosely than its low ones. Intervals can also come out asymmetric around the estimate: percentile intervals inherit whatever shape the bootstrap distribution actually has, rather than imposing $\pm 1.96\,\text{SE}$ symmetry.
Bootstrap the small-$n$ logistic model from the Bayesian page (biomarker slope, $n = 30$) and something instructive happens: some resampled datasets are perfectly separated — every low-biomarker patient a non-responder, every high one a responder — and the refit diverges.
Read the breakage. 115 of 1,000 resamples separate — the bootstrap telling you, bluntly, that these 30 observations sit close to the edge of what logistic regression can estimate. And the usable refits give $(1.51, 5.41)$: right-shifted and asymmetric, agreeing with the flat-prior credible interval $(1.37, 5.53)$ far better than with the symmetric Wald $(0.74, 4.12)$. Two assumption-light methods siding with each other against the quadratic — the full line-up is on Three Routes.
Refinements exist beyond the percentile method used here — BCa (bias-corrected and accelerated) intervals adjust for bias and skew in the bootstrap distribution — but the percentile interval is where every treatment starts, and usually where practice ends.
Real runs (R and numpy use their own random streams, so their digits differ slightly from the seeded interactives above — the method, not the digits, is the point):
set.seed(42)
boot_slopes <- replicate(1000, {
idx <- sample(60, replace = TRUE) # rows, with replacement
coef(lm(y ~ x, data = hr[idx, ]))[2] # refit, keep the slope
})
sd(boot_slopes) # bootstrap SE
quantile(boot_slopes, c(0.025, 0.975)) # percentile CI
set.seed(42) # any statistic at all:
boot_med <- replicate(1000, median(sample(bp, replace = TRUE)))
quantile(boot_med, c(0.025, 0.975))
# sd(boot_slopes) 0.3032
# 2.5% 97.5%
# -1.8486 -0.6641
#
# median BP: 2.5% 97.5%
# 115 126
rng = np.random.default_rng(42)
slopes = []
for _ in range(1000):
idx = rng.integers(0, 60, 60) # rows, with replacement
slopes.append(np.polyfit(x[idx], y[idx], 1)[0])
print(np.std(slopes, ddof=1)) # bootstrap SE
print(np.percentile(slopes, [2.5, 97.5])) # percentile CI
rng = np.random.default_rng(42) # any statistic at all:
meds = [np.median(rng.choice(bp, 80)) for _ in range(1000)]
print(np.percentile(meds, [2.5, 97.5]))
# 0.2909
# [-1.8144 -0.652 ]
#
# [115. 126.]