Skip to main content

Statistical Inference via Data Science A ModernDive into R and the Tidyverse (Second Edition)

Chester Ismay, Albert Y. Kim, Arturo Valdivia

Section 8.2 Estimation with the bootstrap

In 1979, Brad Efron published an article [8] introducing a method called the bootstrap that is next summarized. A random sample of size \(n\) is taken from the population. This sample is used to find another sample, with replacement, also of size \(n\text{.}\) This is called resampling with replacement and the resulting sample is called a bootstrap sample. For example, if the original sample is \(\{4, 2, 5, 4, 1, 3, 7, 4, 6, 1\}\text{,}\) one particular bootstrap sample could be \(\{6, 4, 7, 4, 2, 7, 2, 5, 4, 1\}\text{.}\) Observe that the number 7 appears once in the original sample, but twice in the bootstrap sample; similarly, the number 3 in the original sample does not appear in the bootstrap sample. This is not uncommon for a bootstrap sample, some of the numbers in the original sample are repeated and others are not included.
The basic idea of the bootstrap is to gain a large number of bootstrap samples, all drawn from the same original sample. Then, we use all these bootstrap samples to find estimates of population parameters, standard errors, or even the density curve of the population. Using them we can construct confidence intervals, perform hypothesis testing, and other inferential methods.
This method takes advantage of the large number of bootstrap samples that can be determined. In several respects, this exercise is not different from the sampling distribution explained in Chapter 7. The only difference, albeit an important one, is that we are not sampling from the population, we are sampling from the original sample. How many different bootstrap samples could we get from a single sample? A very large number, actually. If the original sample has 10 numbers, as the one shown above, each possible bootstrap sample of size 10 is determined by sampling 10 times with replacement, so the total number of bootstrap samples is \(10^{10}\) or 10 billion different bootstrap samples. If the original sample has 20 numbers, the number of bootstrap samples is \(20^{20}\text{,}\) a number greater than the total number of stars in the universe. Even with modern powerful computers, it would be an onerous task to calculate every possible bootstrap sample. Instead, a thousand or so bootstrap samples are retrieved, similar to the simulations performed in Chapter 7, and this number is often large enough to provide useful results.
Since Efron [8] proposed the bootstrap, the statistical community embraced this method. During the 1980s and 1990s, many theoretical and empirical results were presented showing the strength of bootstrap methods. As an illustration, Efron [8], Hall [9], Efron and Tibshirani [10], and Hall [11] showed that bootstrapping was at least as good if not better than existent methods, when the goal was to estimate the standard error of an estimator or find the confidence intervals of a parameter. Modifications were proposed to improve the algorithm in situations where the basic method was not producing accurate results. With the continuous improvement of computing power and speed, and the advantages of having ready-to-use statistical software for its implementation, the use of the bootstrap has become more and more popular in many fields.
As an illustration, if we are interested in the mean of the population, \(\mu\text{,}\) and we have collected one random sample, we can gain a large number of bootstrap samples from this original sample, use them to calculate sample means, order the sample means from smallest to largest, and choose the interval that contains the middle 95% of these sample means. This will be the simplest way to find a confidence interval based on the bootstrap. In the next few subsections, we explore how to incorporate this and similar methods to construct confidence intervals.

Subsection 8.2.1 Bootstrap samples: revisiting the almond activity

To study and understand the behavior of bootstrap samples, we return to our example of the chocolate-covered almonds in a bowl. Recall that the bowl is considered the population of almonds, and we are interested in estimating the population mean weight of almonds.
As we did before, we only have access to a single random sample. In this section, we use the data frame almonds_sample_100, a random sample of 100 almonds taken earlier. We call this the original sample, and it is used in this section to create the bootstrap samples. The first 10 rows are shown: 
almonds_sample_100

# A tibble: 100 Ă— 3
# Groups:   replicate [1]
   replicate    ID weight
       <int> <int>  <dbl>
 1         1   166    4.2
 2         1  1215    4.2
 3         1  1899    3.9
 4         1  1912    3.8
 5         1  4637    3.3
 6         1   511    3.5
 7         1   127    4  
 8         1  4419    3.5
 9         1  4729    4.2
10         1  2574    4.1
# â„ą 90 more rows

Constructing a bootstrap sample: resampling once.

We start by constructing one bootstrap sample of 100 almonds from the original sample of 100 almonds. These are the steps needed to perform this task manually:
Step 1: Place the original sample of 100 almonds into a bag or hat.
Step 2: Mix the bag contents, draw one almond, weigh it, and record the weight as seen in Figure 8.2.1.
Photo of a single chocolate-covered almond being weighed on a scale.
Figure 8.2.1. Step 2: Weighing one almond at random.
Step 3: Put the almond back into the bag! In other words, replace it as seen in Figure 8.2.2.
Photo illustrating the act of replacing an item back into a container after sampling.
Figure 8.2.2. Step 3: Replacing almond.
Step 4: Repeat Steps 2 and 3 a total of 99 more times, resulting in 100 weights.
These steps describe resampling with replacement, and the resulting sample is called a bootstrap sample. This procedure results in some almonds being chosen more than once and other almonds not being chosen at all. Resampling with replacement induces sampling variation, so every bootstrap sample can be different than any other.
This activity can be performed manually following the steps described above. We can also take advantage of the R code we have introduced in Chapter 7 and do this virtually. The data frame almonds_sample_100 contains the random sample of almonds taken from the population. We show selected rows from this sample. 
almonds_sample_100 <- almonds_sample_100 |>
  ungroup() |>
  select(-replicate)
almonds_sample_100

# A tibble: 100 Ă— 2
      ID weight
   <int>  <dbl>
 1   166    4.2
 2  1215    4.2
 3  1899    3.9
 4  1912    3.8
 5  4637    3.3
 6   511    3.5
 7   127    4  
 8  4419    3.5
 9  4729    4.2
10  2574    4.1
# â„ą 90 more rows
We use ungroup() and select to eliminate the variable replicate from the almonds_sample_100 as this variable may create clutter when resampling. We can now create a bootstrap sample also of size 100 by resampling with replacement once. 
set.seed(202)
boot_sample <- almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1)
We have used this type of R syntax many times in Chapter 7. We first select the data frame almonds_sample_100 that contains the almonds’ weights in the original sample. We then perform resampling with replacement once: we resample by using rep_sample_n(), a sample of size 100 by setting size = 100, with replacement by adding the argument replace = TRUE, and one time by setting reps = 1. The object boot_sample is a bootstrap sample of 100 almonds’ weights gained from the original sample of 100 almonds’ weights. We show the first ten rows of boot_sample: 
boot_sample

# A tibble: 100 Ă— 3
# Groups:   replicate [1]
   replicate    ID weight
       <int> <int>  <dbl>
 1         1   511    3.5
 2         1   166    4.2
 3         1  2574    4.1
 4         1  1899    3.9
 5         1   127    4  
 6         1   166    4.2
 7         1  4637    3.3
 8         1  4729    4.2
 9         1   511    3.5
10         1  1215    4.2
# â„ą 90 more rows
We can also study some of the characteristics of this bootstrap sample, such as its sample mean: 
boot_sample |>
  summarize(mean_weight = mean(weight))

# A tibble: 1 Ă— 2
  replicate mean_weight
      <int>       <dbl>
1         1        3.69
By using summarize() and mean() on the bootstrap sample boot_sample, we determine that the mean weight is 3.69 grams. Recall that the sample mean of the original sample was found in the previous subsection as 3.68 grams. So, the sample mean of the bootstrap sample is likely different than the sample mean of the original sample. This variation is induced by resampling with replacement, the method for finding the bootstrap sample. We can also compare the histogram of weights for the bootstrap sample with the histogram of weights for the original sample. 
ggplot(boot_sample, aes(x = weight)) +
  geom_histogram(binwidth = 0.1, color = "white") +
  labs(title = "Resample of 100 weights")
ggplot(almonds_sample_100, aes(x = weight)) +
  geom_histogram(binwidth = 0.1, color = "white") +
  labs(title = "Original sample of 100 weights")
Two stacked histograms: top shows the bootstrap resample of 100 almonds’ weights, bottom shows the original sample. The shapes are roughly similar but not identical.
Figure 8.2.3. Comparing weight in the resampled boot_sample with the original sample almonds_sample_100.
Observe in Figure 8.2.3 that while the general shapes of both distributions of weights are roughly similar, they are not identical. This is the typical behavior of bootstrap samples. They are samples that have been determined from the original sample, but because replacement is used before each new observation is attained, some values often appear more than once while others often do not appear at all.

Subsection 8.2.2 Many bootstrap samples: resampling multiple times

In this subsection, we take full advantage of resampling with replacement by taking many bootstrap samples and study relevant information, such as the variability of their sample means. We can start by using the R syntax we used before, this time for 35 replications. 
set.seed(20)
bootstrap_samples_35 <- almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 35)
bootstrap_samples_35

# A tibble: 3,500 Ă— 3
# Groups:   replicate [35]
   replicate    ID weight
       <int> <int>  <dbl>
 1         1  1215    4.2
 2         1   166    4.2
 3         1  2574    4.1
 4         1  1899    3.9
 5         1   127    4  
 6         1  4637    3.3
 7         1  4419    3.5
 8         1  4729    4.2
 9         1   511    3.5
10         1  4729    4.2
# â„ą 3,490 more rows
The syntax is the same as before, but this time we set reps = 35 to get 35 bootstrap samples. The resulting data frame, bootstrap_samples_35, has \(35 \times 100 = 3500\) rows corresponding to 35 resamples of 100 almonds’ weights. Let’s now compute the resulting 35 sample means: 
boot_means <- bootstrap_samples_35 |>
  summarize(mean_weight = mean(weight))
boot_means

# A tibble: 35 Ă— 2
   replicate mean_weight
       <int>       <dbl>
 1         1        3.68
 2         2        3.69
 3         3        3.63
 4         4        3.68
 5         5        3.68
 6         6        3.67
 7         7        3.68
 8         8        3.71
 9         9        3.64
10        10        3.68
# â„ą 25 more rows
Observe that boot_means has 35 rows, corresponding to the 35 bootstrap sample means. Furthermore, observe that the values of mean_weight vary as shown in Figure 8.2.4. 
ggplot(boot_means, aes(x = mean_weight)) +
  geom_histogram(binwidth = 0.01, color = "white") +
  labs(x = "sample mean weight in grams")
A coarse histogram of 35 bootstrap sample mean weights, centered around 3.68 grams.
Figure 8.2.4. Distribution of 35 sample means from 35 bootstrap samples.
This histogram highlights the variation of the sample mean weights. Since we have only used 35 bootstrap samples, the histogram looks a little coarse. To improve our perception of this variation, we find 1000 bootstrap samples and their sample means: 
# Retrieve 1000 bootstrap samples
bootstrap_samples <- almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1000)

# Compute sample means from the bootstrap samples
boot_means <- bootstrap_samples |>
  summarize(mean_weight = mean(weight))
We can combine these two operations into a single chain of pipe (|>) operators: 
set.seed(20)
boot_means <- almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1000) |>
  summarize(mean_weight = mean(weight))
boot_means

# A tibble: 1,000 Ă— 2
   replicate mean_weight
       <int>       <dbl>
 1         1        3.68
 2         2        3.69
 3         3        3.63
 4         4        3.68
 5         5        3.68
 6         6        3.67
 7         7        3.68
 8         8        3.71
 9         9        3.64
10        10        3.68
# â„ą 990 more rows
The data frame boot_means contains 1000 sample mean weights. Each is calculated from a different bootstrap sample and visualized in Figure 8.2.5. 
ggplot(boot_means, aes(x = mean_weight)) +
  geom_histogram(binwidth = 0.01, color = "white") +
  labs(x = "sample mean weight in grams")
A bell-shaped histogram of 1000 bootstrap sample mean weights, centered around 3.68 grams.
Figure 8.2.5. Histogram of 1000 bootstrap sample mean weights of almonds.
The histogram is a graphical approximation of the bootstrap distribution of the sample mean. This distribution is constructed by getting all the sample means from every bootstrap sample constructed based on the original sample. Since the total number of possible bootstraps is really large, we have not used all of them here, but 1000 of them already provides a good visual approximation.
Observe also that the bootstrap distribution itself can approximate the sampling distribution of the sample mean, a concept we studied in Chapter 7 where we took multiple samples from the population. The key difference here is that we resample from a single sample, the original sample, not from the entire population.
By inspecting the histogram in Figure 8.2.5, the bell-shape is apparent. We can also approximate the center and the spread of this distribution by computing the mean and the standard deviation of these 1000 bootstrap sample means: 
boot_means |>
  summarize(mean_of_means = mean(mean_weight),
            sd_of_means = sd(mean_weight))

# A tibble: 1 Ă— 2
  mean_of_means sd_of_means
          <dbl>       <dbl>
1          3.68      0.0357
Everything we learned in Chapter 7 when studying the sampling distribution of the sample mean applies here. For example, observe that the mean of these bootstrap sample means is near 3.68 grams, very close to the mean of the original sample. Our intention is not to study the distribution of the bootstrap samples, but rather to use them to estimate population values, such as the population mean. In the next section, we discuss how can we use these bootstrap samples to construct confidence intervals.

Exercises Exercises

1.
What is the chief difference between a bootstrap distribution and a sampling distribution?
3.
Which of the following is true about the confidence level when constructing a confidence interval?
  • A. The confidence level determines the width of the interval and affects how likely it is to contain the population parameter.
  • B. The confidence level is always fixed at 95% for all statistical analyses involving confidence intervals.
  • C. A higher confidence level always results in a narrower confidence interval, making it more useful for practical purposes.
  • D. The confidence level is only relevant when the population standard deviation is known.
4.
How does increasing the sample size affect the width of a confidence interval for a given confidence level?
  • A. It increases the width of the confidence interval, making it less precise.
  • B. It has no effect on the width of the confidence interval since the confidence level is fixed.
  • C. It decreases the width of the confidence interval, making it more precise by reducing the margin of error.
  • D. It changes the confidence level directly, regardless of other factors.

Subsection 8.2.3 Confidence intervals and the bootstrap: original workflow

The process of determining bootstrap samples and using them for estimation is called bootstrapping. We can estimate population parameters such as the mean, median, or standard deviation. We can also construct confidence intervals.
In this subsection, we focus on the latter and construct confidence intervals based on bootstrap samples. For this, we review the R syntax and workflow we have already used in previous sections and also introduce a new package: the infer package for tidy and transparent statistical inference.

Original workflow.

Recall that we took bootstrap samples, then calculated the sample means from these samples. Let’s revisit the original workflow using dplyr verbs and the |> operator.
First, we use rep_sample_n() to resample from the original sample almonds_sample_100 of 100 almonds. We set size = 100 to generate bootstrap samples of the same size as the original sample, and we resample with replacement by setting replace = TRUE. We create 1000 bootstrap samples by setting reps = 1000: 
almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1000)
Second, we add another pipe followed by summarize() to compute the sample mean() weight for each replicate: 
almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1000) |>
  summarize(mean_weight = mean(weight))
For this simple case, all we needed was to use the rep_sample_n() function and a dplyr verb. However, using only dplyr verbs provides us with a limited set of tools that is not ideal when working with more complicated situations. This is the reason we introduce the infer package.

Subsection 8.2.4 The infer package workflow

The infer package is an R package for statistical inference. It makes efficient use of the |> pipe operator we introduced in Chapter 3 to spell out the sequence of steps necessary to perform statistical inference in a “tidy” and transparent fashion. Just as the dplyr package provides functions with verb-like names to perform data wrangling, the infer package provides functions with intuitive verb-like names to perform statistical inference, such as constructing confidence intervals or performing hypothesis testing. We have discussed the theory-based implementation of the former in Subsection 8.1.4 and we introduce the latter in Chapter 9.
Using the example of almonds’ weights, we introduce infer first by comparing its implementation with dplyr. Recall that to calculate a sample statistic or point estimate from a sample, such as the sample mean, when using dplyr we use summarize() and mean(): 
almonds_sample_100 |>
  summarize(stat = mean(weight))

# A tibble: 1 Ă— 1
   stat
  <dbl>
1  3.68
If we want to use infer instead, we use the functions specify() and calculate() as shown: 
almonds_sample_100 |>
  specify(response = weight) |>
  calculate(stat = "mean")

# A tibble: 1 Ă— 1
   stat
  <dbl>
1  3.68
The new structure using infer seems slightly more complicated than the one using dplyr for this simple calculation. These functions will provide three chief benefits moving forward.
  • First, the infer verb names better align with the overall simulation-based framework you need to understand to construct confidence intervals and to conduct hypothesis tests (in Chapter 9). We will see flowchart diagrams of this framework in the upcoming Figure 8.2.11 and in Chapter 9.
  • Second, you can transition seamlessly between confidence intervals and hypothesis testing with minimal changes to your code. This becomes apparent in Subsection 9.4.2 when we compare the infer code for both of these inferential methods.
  • Third, the infer workflow is much simpler for conducting inference when you have more than one variable. We introduce two-sample inference where the sample data is collected from two groups, such as in Section 8.4 where we study the contagiousness of yawning and in Section 9.2 where we compare the popularity of music genres.
We now illustrate the sequence of verbs necessary to construct a confidence interval for \(\mu\text{,}\) the population mean weight of almonds.

Subsubsection 8.2.4.1 1. specify variables

A flowchart diagram showing the specify step in the infer workflow, highlighting the response variable selection.
Figure 8.2.6. Diagram of the specify() verb.
As shown in Figure 8.2.6, the specify() function is used to choose which variables in a data frame are the focus of our statistical inference. We do this by specifying the response argument. For example, in our almonds_sample_100 data frame of the 100 almonds sampled from the bowl, the variable of interest is weight: 
almonds_sample_100 |>
  specify(response = weight)

Response: weight (numeric)
# A tibble: 100 Ă— 2
      ID weight
   <int>  <dbl>
 1   166    4.2
 2  1215    4.2
 3  1899    3.9
 4  1912    3.8
 5  4637    3.3
 6   511    3.5
 7   127    4  
 8  4419    3.5
 9  4729    4.2
10  2574    4.1
# â„ą 90 more rows
Notice how the data itself does not change, but the Response: weight (numeric) meta-data does. This is similar to how the group_by() verb from dplyr doesn’t change the data, but only adds “grouping” meta-data, as we saw in Chapter 3.
We can also specify which variables are the focus of the study by introducing a formula = y ~ x in specify(). This is the same formula notation you saw in Chapter 5 and Chapter 6 on regression models: the response variable y is separated from the explanatory variable x by a ~ (“tilde”). The following use of specify() with the formula argument yields the same result seen previously: 
almonds_sample_100 |>
  specify(formula = weight ~ NULL)
In the case of almonds we only have a response variable and no explanatory variable of interest. Thus, we set the x on the right-hand side of the ~ to be NULL. In cases where inference is focused on a single sample, as it is the almonds’ weights example, either specification works. In the examples we present in later sections, the formula specification is simpler and more flexible. In particular, this comes up in the upcoming Section 8.4 on comparing two proportions and Section 10.3 on inference for regression.

Subsubsection 8.2.4.2 2. generate replicates

A flowchart diagram showing the generate step in the infer workflow, illustrating bootstrap resampling.
Figure 8.2.7. Diagram of generate() replicates.
After we specify() the variables of interest, we pipe the results into the generate() function to generate replicates. This is the function that produces the bootstrap samples or performs the similar resampling process a large number of times, based on the variable(s) specified previously, as shown in Figure 8.2.7. Recall we did this 1000 times.
The generate() function’s first argument is reps, which sets the number of replicates we would like to generate. Since we want to resample the 100 almonds in almonds_sample_100 with replacement 1000 times, we set reps = 1000. The second argument type determines the type of computer simulation used. Setting this to type = "bootstrap" produces bootstrap samples using resampling with replacement. We present different options for type in Chapter 9. 
almonds_sample_100 |>
  specify(response = weight) |>
  generate(reps = 1000, type = "bootstrap")

# A tibble: 100,000 Ă— 2
# Groups:   replicate [1000]
   replicate weight
       <int>  <dbl>
 1         1    3.6
 2         1    3.4
 3         1    4.2
 4         1    3.4
 5         1    3.5
 6         1    3.4
 7         1    3.8
 8         1    3.7
 9         1    4.1
10         1    3.4
# â„ą 99,990 more rows
Observe that the resulting data frame has 100,000 rows. This is because we have found 1000 bootstrap samples, each with 100 rows.
The variable replicate indicates the bootstrap sample each row belongs to, from 1 to 1000, each replicate repeated 100 times. The default value of the type argument is "bootstrap" in this scenario, so the inclusion was only made for completeness. If the last line was written simply as generate(reps = 1000), the result would be the same.
Comparing with original workflow: Note that the steps of the infer workflow so far produce the same results as the original workflow using the rep_sample_n() function we saw earlier. In other words, the following two code chunks produce similar results: 
# infer workflow:
almonds_sample_100 |>
  specify(response = weight) |>
  generate(reps = 1000)

# Original workflow:
almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1000)

Subsubsection 8.2.4.3 3. calculate summary statistics

A flowchart diagram showing the calculate step in the infer workflow, where a summary statistic is computed for each bootstrap sample.
Figure 8.2.8. Diagram of calculate() summary statistics.
After we generate() 1000 bootstrap samples, we want to summarize each of them, for example, by calculating the sample mean of each one of them. As Figure 8.2.8 shows, the calculate() function does this.
In our example, we calculate the mean weight for each bootstrap sample by setting the stat argument equal to "mean" inside the calculate() function. The stat argument can be used for other common summary statistics such as "median", "sum", "sd" (standard deviation), and "prop" (proportion). To see a list of other possible summary statistics you can use, type ?calculate and read the help file.
Let’s save the result in a data frame called bootstrap_means and explore its contents: 
bootstrap_means <- almonds_sample_100 |>
  specify(response = weight) |>
  generate(reps = 1000) |>
  calculate(stat = "mean")
bootstrap_means

# A tibble: 1,000 Ă— 2
   replicate  stat
       <int> <dbl>
 1         1  3.68
 2         2  3.69
 3         3  3.63
 4         4  3.68
 5         5  3.68
 6         6  3.68
 7         7  3.68
 8         8  3.71
 9         9  3.64
10        10  3.68
# â„ą 990 more rows
Observe that the resulting data frame has 1000 rows and 2 columns corresponding to the 1000 replicate values. It also has the mean weight for each bootstrap sample saved in the variable stat.
Comparing with original workflow: You may have recognized at this point that the calculate() step in the infer workflow produces the same output as the summarize() step in the original workflow. 
# infer workflow:
almonds_sample_100 |>
  specify(response = weight) |>
  generate(reps = 1000) |>
  calculate(stat = "mean")

# Original workflow:
almonds_sample_100 |>
  rep_sample_n(size = 100, replace = TRUE, reps = 1000) |>
  summarize(stat = mean(weight))

Subsubsection 8.2.4.4 4. visualize the results

A flowchart diagram showing the visualize step in the infer workflow, displaying the bootstrap distribution as a histogram.
Figure 8.2.9. Diagram of visualize() results.
The visualize() verb provides a quick way to visualize the bootstrap distribution as a histogram of the numerical stat variable’s values as shown in Figure 8.2.10. The pipeline of the main infer verbs used for exploring bootstrap distribution results is shown in Figure 8.2.11. 
visualize(bootstrap_means)
Histogram of 1000 bootstrap sample means of almond weights, produced by the visualize() function from infer.
Figure 8.2.10. Bootstrap distribution.
Comparing with original workflow: In fact, visualize() is a wrapper function for the ggplot() function that uses a geom_histogram() layer. Recall that we illustrated the concept of a wrapper function in Figure 5.1.7 in Subsection 5.1.2. 
# infer workflow:
visualize(bootstrap_means)

# Original workflow:
ggplot(bootstrap_means, aes(x = stat)) +
  geom_histogram()
The visualize() function can take many other arguments to customize the plot further. In later sections, we take advantage of this flexibility. In addition, it works with helper functions to add shading of the histogram values corresponding to the confidence interval values. We have introduced the different elements of the infer workflow for constructing a bootstrap distribution and visualizing it. A summary of these steps is presented in Figure 8.2.11.
Flowchart showing the four steps of the infer workflow for confidence intervals: specify, generate, calculate, and visualize.
Figure 8.2.11. infer package workflow for confidence intervals.

Subsection 8.2.5 Confidence intervals using bootstrap samples with infer

We are ready to introduce confidence intervals using the bootstrap via infer. We present two different methods for constructing 95% confidence intervals as interval estimates of an unknown population parameter: the percentile method and the standard error method. Let’s now check out the infer package code that explicitly constructs these.

Subsubsection 8.2.5.1 Percentile method

Recall that in Subsection 8.2.4 we have generated 1000 bootstrap samples and stored them in data frame bootstrap_means: 
bootstrap_means

# A tibble: 1,000 Ă— 2
   replicate  stat
       <int> <dbl>
 1         1  3.68
 2         2  3.69
 3         3  3.63
 4         4  3.68
 5         5  3.68
 6         6  3.68
 7         7  3.68
 8         8  3.71
 9         9  3.64
10        10  3.68
# â„ą 990 more rows
The sample means stored in bootstrap_means represent a good approximation to the bootstrap distribution of all possible bootstrap samples. The percentile method for constructing 95% confidence intervals sets the lower endpoint of the confidence interval at the 2.5th percentile of bootstrap_means and similarly sets the upper endpoint at the 97.5th percentile. The resulting interval captures the middle 95% of the values of the sample mean weights of almonds in bootstrap_means. This is the interval estimate of the population mean weight of almonds in the entire bowl.
We can compute the 95% confidence interval by piping bootstrap_means into the get_confidence_interval() function from the infer package, with the confidence level set to 0.95 and the confidence interval type to be "percentile". We save the results in percentile_ci. 
percentile_ci <- bootstrap_means |>
  get_confidence_interval(level = 0.95, type = "percentile")
percentile_ci

# A tibble: 1 Ă— 2
  lower_ci upper_ci
     <dbl>    <dbl>
1     3.61     3.76
Alternatively, we can visualize the interval by piping the bootstrap_means data frame into the visualize() function and adding a shade_confidence_interval() layer. We set the endpoints argument to be percentile_ci. 
visualize(bootstrap_means) +
  shade_confidence_interval(endpoints = percentile_ci)
Histogram of 1000 bootstrap sample mean weights with the middle 95% shaded, representing the percentile method confidence interval. Vertical lines mark the 2.5th and 97.5th percentiles.
Figure 8.2.12. Percentile method 95% confidence interval shaded corresponding to potential values.
Observe in Figure 8.2.12 that 95% of the sample means stored in the stat variable in bootstrap_means fall between the two endpoints marked with the darker lines, with 2.5% of the sample means to the left of the shaded area and 2.5% of the sample means to the right. You also have the option to change the colors of the shading using the color and fill arguments.
The infer package has incorporated a shorter named function shade_ci() that produces the same results. Try out the following code: 
visualize(bootstrap_means) +
  shade_ci(endpoints = percentile_ci, color = "hotpink", fill = "khaki")

Subsubsection 8.2.5.2 Standard error method

In Subsection 8.1.4 we introduced theory-based confidence intervals. We show that a 95% confidence interval can be constructed as
\begin{equation*} \left(\overline{x} - 1.96 \cdot SE(\bar{x}),\quad \overline{x} + 1.96 \cdot SE(\bar{x})\right) \end{equation*}
where \(\overline{x}\) is the sample mean of the original sample, 1.96 is the number of standard errors around the mean needed to account for 95% of the area under the density curve (when the distribution is normal), and \(SE(\bar{x})\) is the standard error of the sample mean that can be computed as \(\sigma/\sqrt{n}\) if the population standard deviation is known, or estimated as \(s/\sqrt{n}\) if we have to use the sample standard deviation, \(s\text{,}\) and the sample size, \(n\text{.}\)
We use the same structure to construct confidence intervals but using the bootstrap sample means to estimate the standard error of \(\overline{x}\text{.}\) Thus, the 95% confidence interval for the population mean, \(\mu\text{,}\) using the standard error estimated via bootstrapping, \(SE_{\text{boot}}\text{,}\) is:
\begin{equation*} \left(\overline{x} - 1.96 \cdot SE_{\text{boot}},\quad \overline{x} + 1.96 \cdot SE_{\text{boot}}\right) \end{equation*}
We can compute this confidence interval using dplyr. First, we calculate the estimated standard error: 
SE_boot <- bootstrap_means |>
  summarize(SE = sd(stat)) |>
  pull(SE)
SE_boot

[1] 0.03566146
and then use the original sample mean to calculate the 95% confidence interval: 
almonds_sample_100 |>
  summarize(lower_bound = mean(weight) - 1.96 * SE_boot,
            upper_bound = mean(weight) + 1.96 * SE_boot)

# A tibble: 1 Ă— 2
  lower_bound upper_bound
        <dbl>       <dbl>
1        3.61        3.75
Alternatively, computation of the 95% confidence interval can once again be done via infer. We find the sample mean of the original sample and store it in variable x_bar: 
x_bar <- almonds_sample_100 |>
  specify(response = weight) |>
  calculate(stat = "mean")
x_bar

# A tibble: 1 Ă— 1
   stat
  <dbl>
1  3.68
Now, we pipe the bootstrap_means data frame we created into the get_confidence_interval() function. We set the type argument to be "se" and specify the point_estimate argument to be x_bar in order to set the center of the confidence interval to the sample mean of the original sample. 
standard_error_ci <- bootstrap_means |>
  get_confidence_interval(type = "se", point_estimate = x_bar, level = 0.95)
standard_error_ci

# A tibble: 1 Ă— 2
  lower_ci upper_ci
     <dbl>    <dbl>
1     3.61     3.75
The results are the same whether dplyr or infer is used, but as explained earlier, the latter provides more flexibility for other tests.
If we would like to visualize the interval, we can once again pipe the bootstrap_means data frame into the visualize() function and add a shade_confidence_interval() layer to our plot. We set the endpoints argument to be standard_error_ci. The resulting standard-error method based on a 95% confidence interval for \(\mu\) can be seen in Figure 8.2.13. 
visualize(bootstrap_means) +
  shade_confidence_interval(endpoints = standard_error_ci)
Histogram of 1000 bootstrap sample mean weights with the standard-error method 95% confidence interval shaded. The interval is centered at the sample mean with width determined by 1.96 bootstrap standard errors.
Figure 8.2.13. Standard-error method 95% confidence interval.
Because we are using bootstrap samples to construct these intervals, we call the percentile and standard error methods simulation-based methods. We can compare the 95% confidence intervals from using both simulation-based methods as well as the one attained using the theory-based method described in Subsection 8.1.4:
  • Percentile method: approximately (3.61, 3.76)
  • Standard error method: approximately (3.61, 3.75)
  • Theory-based method: (3.653, 3.707)

Exercises 8.2.5.3 Exercises

1.
Construct a 95% confidence interval for the median weight of all almonds with the percentile method. Is it appropriate to also use the standard-error method?
3.
What is the main purpose of bootstrapping in statistical inference?
  • A. To visualize data distributions and identify outliers.
  • B. To generate multiple samples from the original data for estimating parameters.
  • C. To replace missing data points with the mean of the dataset.
  • D. To validate the assumptions of a regression model.
5.
What is a key difference between the percentile method and the standard error method for constructing confidence intervals using bootstrap samples?
  • A. The percentile method requires the population standard deviation.
  • B. The percentile method uses the middle 95% of bootstrap statistics, while the standard error method relies on the estimated standard error.
  • C. The standard error method always results in a narrower confidence interval.
  • D. The percentile method requires more bootstrap samples.
PrevTopNext