12  Calibrating Interval Estimates using the Bootstrap

Review

$$ \newcommand{X}{} \newcommand{Y}{} % A hex digit after # inside a macro body is eaten by build-time expansion, so digit-leading hexes spell # via .

$$

To calibrate an interval estimate, we add ‘arms’ of equal length to our point estimate. We choose the arm length so that, in 95% of surveys done like ours, the estimation target is within them. We can’t choose that length without the sampling distribution of our point estimate, and we don’t have it. So we estimate it. That’s what ?sec-calibration-with-the-binomial does: plug a sample frequency \hat\theta into the Binomial formula from ?sec-sampling-with-replacement and read the width off the result.

\hat P\p*{\sum_{i=1}^n Y_i = s} = \binom{n}{s} \hat\theta^{s} (1-\hat\theta)^{n-s}

Both of the sampling distributions we’ve worked out came to us as a parametric form. By parametric form, I mean a formula in terms of some parameters. That tells us what parameters we need to estimate to estimate the sampling distribution. Sampling with replacement, the Binomial’s one parameter is the population frequency \theta, which is what we’re trying to estimate anyway. Sampling without replacement, ?sec-sampling-without-replacement, it’s the Hypergeometric, and there are two parameters we don’t know—the population size m and the number of ones m_1—so we plug in m\hat\theta for the second.

That’s what we’re about to give up.

An Exercise

Bootstrapping by Hand

Here’s an exercise you can do with classmates, friends, or—if you don’t have any—imaginary characters you act out yourself.

Suppose you’ve run a small survey. You’ve polled 6 people, asking whether they’ll vote. Here’s what you got.

\begin{array}{r|rrrrrr|r} i & 1 & 2 & 3 & 4 & 5 & 6 & \bar{Y}_{6} \\ Y_i & 0 & 0 & 0 & 1 & 1 & 1 & 0.5 \\ \end{array}

We’re going to estimate the sampling distribution in a new way. It’ll be a lot like calculating the actual sampling distribution: we’ll draw a sample of size n with replacement, use it to calculate our estimator, and repeat. What’s different is the population we’re drawing our sample from. We don’t have responses from the whole population, so we’ll draw from the closest thing we’ve got: the sample. We call this bootstrapping. The estimate of the sampling distribution we get is called the bootstrap sampling distribution.

Exercise 14.1 (Drawing from the Bootstrap Sampling Distribution)  

Exercise

Here’s how you draw a sample from the bootstrap sampling distribution of our point estimator.

1. Roll a 6-sided die 6 times to draw a sample of size n=6 from the sample.
2. Calculate the sample frequency. That’s one draw from the bootstrap sampling distribution. Write it down.

If you do this 5 times each, you’ll have a bunch of draws. Tally them up and visualize the result as a histogram—a histogram of draws from the bootstrap sampling distribution.

Simulating a Few More Draws

If you have access to a computer, you can simulate many more draws.

resample = function(Y) {
  Jstar = sample(seq_along(Y), length(Y), replace=TRUE)
  Y[Jstar]
}
estimator = function(Y) mean(Y)

draws = rep(list(Yexercise), 10000) |>
  map(resample) |>
  map_vec(estimator) |>
  with_seed()
computer.tally = as.data.frame(table(draws)) |>
  rename(theta.hat=draws, count=Freq)

Programming the bootstrap

The bootstrap pipeline has the same three parts as the sampling-distribution pipeline: repeat a source, draw from it, and calculate an estimator. The difference is the source and the drawing function. A sampling distribution uses the population and draw.sample(); a bootstrap distribution uses the observed sample and resample().

Naming resample and estimator makes that substitution visible in the code. The programming reference places the two pipelines side by side.

Discussion

Put the tally from your rolls into count and run it. Your histogram goes up against our binomial estimate of the sampling distribution.

Figure 14.1

If all has gone to plan, your histogram looks a lot like our binomial estimate of the sampling distribution. Substituting in a computer tally of 10,000 draws from the bootstrap sampling distribution, Figure 14.1, we nail it. It looks like the bootstrap sampling distribution is the binomial estimate.

Exercise 14.2 (Why the Bootstrap Matches the Binomial)  

Exercise

Why is that? Think about it before reading on.

TipHint

We know the distribution of the frequency of 1s in a sample of size n drawn with replacement from a population y_1\ldots y_m of binary responses with frequency \theta. It’s \text{Binomial}(n, \theta). To estimate it, we plug in \hat\theta, the frequency of 1s in our sample.

So our Binomial estimate is the distribution of the frequency of 1s in a sample of size n drawn with replacement from ‘a population’ Y_1 \ldots Y_n with frequency \hat\theta.

What’s a draw from the bootstrap sampling distribution?

SolutionWhy the Bootstrap Matches the Binomial

Each bootstrap sample is the frequency of 1s in a sample of size n drawn with replacement from ‘a population’ Y_1 \ldots Y_n in which the frequency of 1s is \hat\theta—the frequency of 1s in our sample. Because that ‘population’ is our sample.

The Picture

What we just did draws as a rectangle.

\begin{array}{ccc} \underset{\textstyle\text{drawing from the population}}{\text{sampling distribution}} & = & \text{Binomial}(n, \theta)\big/n \\[18pt] {\overset{\normalsize ?}{\approx}} & & {\Big\downarrow}\ {\footnotesize\text{plug in } \hat\theta} \\[18pt] \underset{\textstyle\text{drawing from the sample}}{\text{bootstrap sampling distribution}} & = & \text{Binomial}(n, \hat\theta)\big/n \end{array}

The top edge is the derivation from ?sec-sampling-with-replacement: sampling with replacement, the sample frequency is Binomial. The right edge is our estimate—one approximation, \hat\theta in place of \theta. The bottom edge is the exercise you just did: the bootstrap sampling distribution is the plug-in estimate. And the left edge is the one we care about: does the bootstrap estimate the real sampling distribution? Follow the other three edges around the rectangle and you have the answer. The left edge inherits the right edge’s error, and only that.

The Bootstrap

The Bootstrap Interpretation in Our Turnout Poll

Figure 15.1

Population. Every registered voter in Georgia, one faint point each; the fraction who will actually vote is the number we want.

Sample. 625 people drawn from the population — so large that nobody gets called twice, so every dot keeps the same size. The sample’s frequency is our estimate.

Bootstrap Sample. 625 draws taken from the sample itself, with replacement. People get called more than once (dot area grows with calls) and some never get called; this resampled group stands in for the population we cannot see.

The sampling distribution estimate we’ve used was \text{Binomial}(n,\hat\theta) for n=625 and \hat\theta \approx 0.68. It’s the distribution of the proportion of heads in 625 flips of a coin with probability \hat\theta \approx 0.68 of heads, where \hat\theta \approx 0.68 is the proportion of voters we’ve polled who will vote.

That is, it’s the sampling distribution of a ‘poll’ of the people in our sample: roll a 625-sided die 625 times, call up the corresponding person in our sample, and count up the yeses we hear. This is random because we’re drawing with replacement. Each time we run this poll, we call each person in our sample 0, 1, 2, … times, and the number of times we call them is random.

If we plot our voters on a map, you can see the idea in visual terms. On the left, we have the population. In the middle, we have our sample, drawn with replacement from the population. On the right, we have a bootstrap sample—a new sample drawn with replacement from the sample.

Each ‘call’ that a person receives increases the size of their dot: circle area is proportional to number of calls. In the sample, even though it’s drawn with replacement, all dots are the same size—because we draw from such a large population, nobody gets called twice. In the bootstrap sample, dots vary in size—because we draw n people from a sample of size n, it’s almost impossible not to call somebody twice.

Bootstrapping

Before the election, we don’t observe the population. But we do observe the sample. And we can sample from the sample, acting as if it were the population. We’ll take repeated random samples of size 625, with replacement, from our sample of size 625.

We call these bootstrap samples and estimates based on them bootstrap estimates. The distribution of these estimates is called the bootstrap sampling distribution. If the sample is like the population, this should be like our estimator’s actual sampling distribution.

The Sample \begin{array}{r|rrrr|r} i & 1 & 2 & \dots & 625 & \bar{Y}_{625} \\ Y_i & 1 & 1 & \dots & 1 & 0.68 \\ \end{array}

The Bootstrap Sample

\begin{array}{r|rrrr|r} i & 1 & 2 & \dots & 625 & \bar{Y}_{625}^* \\ Y_i^* & 1 & 1 & \dots & 0 & 0.72 \\ \end{array}

The Population

\begin{array}{r|rrrr|r} j & 1 & 2 & \dots & 7.23M & \bar{y}_{7.23M} \\ y_{j} & 1 & 1 & \dots & 1 & 0.70 \\ \end{array}

The ‘Bootstrap Population’ — The Sample \begin{array}{r|rrrr|r} j & 1 & 2 & \dots & 625 & \bar{y}^*_{625} \\ y_j^* & 1 & 1 & \dots & 1 & 0.68 \\ \end{array}

We use stars to distinguish the bootstrap sample from our original sample—we write Y_i^* and \bar Y^*.

The Bootstrap is Nonparametric

resample = function(Y) {
  J.star = sample(seq_along(Y), length(Y), replace=TRUE)
  Y[J.star]
}
estimator = function(Y) mean(Y)

bootstrap.samples = rep(list(Y), 10000) |>
  map(resample) |>
  map_vec(estimator) |>
  with_seed()

We do not need to know the parametric form of our sampling distribution to use the bootstrap. All we do is re-run our poll acting as if our sample were the population. We can do this no matter what we’re estimating. So let’s try it out on something where we don’t know the parametric form.

Beyond Binary Responses

Vote History

So far we’ve asked a binary question: will you vote? But voter files contain richer information. For each registered voter, we can see how many of the last 5 elections they voted in. Let’s call this vote history, and denote it H_i for the people in our sample and h_j for the people in the population.

\begin{array}{r|rrrr|r} i & 1 & 2 & \dots & 625 & \bar{H}_{625} \\ H_i & 1 & 4 & \dots & 2 & 2.8 \\ \end{array}

We might want to estimate the average vote history in the population. Our point estimate would be the sample mean, \bar H_n = 2.8. But what about the sampling distribution?

When the response was binary, we knew the parametric form: Binomial. The mean of n binary responses, each equal to 1 with probability \theta, has a \text{Binomial}(n, \theta)/n distribution. But now our responses take values in \{0, 1, 2, 3, 4, 5\}. There’s no simple parametric form for the distribution of their mean.

We could work one out. A response H_i has some probability p_0 of being 0, p_1 of being 1, and so on. The sum H_1 + \ldots + H_n is a sum of n such random variables, and its distribution could be computed by convolution. But that’s a pain. And it would be an even bigger pain if our responses were continuous.

The bootstrap gives us a way out. We don’t need to know the parametric form. We just resample from the sample.

The Sampling Distribution of Average Vote History

draw.sample = function(h) {
  I = sample(seq_along(h), n, replace=TRUE)
  h[I]
}
history.estimator = function(H) mean(H)

mean.history.samples = rep(list(h), 10000) |>
  map(draw.sample) |>
  map_vec(history.estimator) |>
  with_seed()

Here’s what the sampling distribution of average vote history looks like. As before, if we can estimate it we can use that estimate to get a 95% confidence interval. But, unlike before, we don’t have a nice parametric form like the Binomial. So we’ll use the bootstrap to estimate it.

The Bootstrap Sampling Distribution of Average Vote History

resample = function(H) {
  I = sample(seq_along(H), length(H), replace=TRUE)
  H[I]
}

history.bootstrap.samples = rep(list(H), 10000) |>
  map(resample) |>
  map_vec(history.estimator) |>
  with_seed()

The actual sampling distribution

The bootstrap sampling distribution

It looks like it works in this case. The bootstrap sampling distribution is a good approximation to the actual sampling distribution, and the interval widths are similar.

But we’re no longer able to argue that it should work the way we did before. For the binary case, we took advantage of our knowledge of the parametric form—we knew the bootstrap sampling distribution was the plug-in binomial estimate. And we don’t have that now.

We’ll get there. But we’ll need a few new tools we’ll develop in the coming weeks: normal approximation, which gives us a parametric form for an approximation to our estimator’s sampling distribution, and techniques for variance calculation, which help us understand the parameters that go into it.

The Picture, Broken

\begin{array}{ccc} \underset{\textstyle\text{drawing from the population}}{\text{sampling distribution}} & \approx & \text{normal approximation} \\[18pt] {\overset{\normalsize ?}{\approx}} & & {\Big\downarrow}\ {\footnotesize\text{plug in estimates}} \\[18pt] \underset{\textstyle\text{drawing from the sample}}{\text{bootstrap sampling distribution}} & \approx & \text{its plug-in estimate} \end{array}

Here’s the same failure drawn on the rectangle. For binary responses, the right side was exact: the parametric form. With vote history there is no parametric form, so the right side is broken. The plan is to replace it with the normal approximation—a parametric form for an approximation—which turns the top and bottom edges from = into \approx too. A question to sit with until we get there: we already had one approximation step, plugging in. Once one edge of the rectangle is approximate, how much does it cost that the others are?

Scope of the Bootstrap

The bootstrap resamples with replacement, no matter how the sample itself was drawn. So what it estimates is the with-replacement sampling distribution—the plug-in Binomial. But suppose our poll had sampled without replacement. Then the real sampling distribution is the Hypergeometric from ?sec-sampling-without-replacement:

P_\theta(\hat\theta = t) = \binom{n}{nt}\, \frac{\bigl[m_1(m_1-1)\ldots(m_1-nt+1)\bigr]\,\bigl[m_0(m_0-1)\ldots(m_0-n(1-t)+1)\bigr]}{m(m-1)\ldots(m-n+1)} \quad\text{ for }\quad m_1 = m\theta, \ \ m_0 = m(1-\theta).

Same two stages as before, but calls use people up, so each factor falls as its response repeats. If we knew we’d sampled without replacement, and knew m, we could plug \hat\theta into this instead. The question for the bootstrap is how far apart these two distributions really are.

Comparing Sampling with and without Replacement

Each panel shows the distribution of the estimated population proportion. Rows change the number of calls n; columns change the population size m. The diagonal keeps n/m=1/5. The four cases below change the population proportion \theta_1.

1. Across a row: keep the number of calls fixed. What changes as the population grows? What stays the same?
2. Down the diagonal: increase the number of calls and the population size together. What changes when the fraction sampled stays the same?
3. Between cases: what changes as the population proportion moves away from one half? How do the last two compare?

Even: \theta_1=0.5

Figure 17.1

Mid: \theta_1=0.7

Figure 17.2

Extreme-ish: \theta_1=0.9

Figure 17.3

Flipped: \theta_1=0.1

Figure 17.4

The Hypergeometric looks a lot like the Binomial when n is small relative to m. That’s because sampling without replacement is almost like sampling with replacement when you’re unlikely to pick the same person twice anyway. When n is a meaningful fraction of m, the Hypergeometric is narrower—you get better precision by not wasting calls on people you’ve already talked to.

If you called everybody you’d know \theta_1 exactly, and there’d be nothing left to be uncertain about.