3 - Confidence Intervals with Bootstrapping

Learning Objectives

Upon completing today’s lab activity, students should be able to do the following using R and RStudio:

1. Computing quantiles/percentiles given some data.

2. Using bootstrapping to compute Confidence Intervals (CIs).

library(tidyverse)
library(openintro)
library(infer)
library(ggplot2)
library(gghighlight)

Computing Quantiles given Some Data

We can compute the quantiles or percentiles of a given data using the base R function called quantile.

# given some data (unordered)
x <- c(4,5,5,6,7,7,2,3,4,5,5,5)

# the 50th percentile (also known as the median)
print(quantile(x,0.50))

## 50% 
##   5

# the 2.5th percentile
print(quantile(x,0.025))

##  2.5% 
## 2.275

# the 97.5th percentile
print(quantile(x,0.975))

## 97.5% 
##     7

Case Study - Medical Consultant (Proportion)

One consultant tried to attract patients by noting the average complication rate for liver donor surgeries in the US is about 10%, but her clients have had only 3 complications in the 62 liver donor surgeries she has facilitated.

She claims this is strong evidence that her work meaningfully contributes to reducing complications.

Let \(p\) represent the true complication rate for liver donors working with this consultant. We estimate \(p\) using the data, and label the estimate \(\hat{p}\).

\[\hat{p} = \frac{3}{62} = 0.048 \longrightarrow \text{Observed Statistic}\]

\[p_0 = 0.10 \longrightarrow \text{The Null Statistic - the average complication rate in the US}\] The claim is that there is a causal connection, but the data are observational, so we must be on the lookout for confounding variables. We can’t conclude a causal connection for observational studies.

Objective: Estimate the unknown population proportion by using the sample to approximate the proportion of complications for a client of the medical consultant. Make a 95% confidence interval to give a plausible range of values for the population proportion.

# Even though we don't have the actual data, we have the sample proportion and we can work on that for performing bootstrapping.
set.seed(334422)
bsprop_med <- tibble(
  bsprop = rbinom(10000, size=62, prob=(3/62))/62 # this function simulates the binomial ("success" (complication) or "failure" (no complication))
)

bsprop_med_summary <- bsprop_med %>%
  summarise(
    bsprop_025 = quantile(bsprop, 0.025),
    bsprop_975 = quantile(bsprop, 0.975),
  )

ggplot(bsprop_med, aes(x = bsprop)) +
  geom_histogram(binwidth = 0.0075, fill = "gray") +
  annotate("segment", 
           x = bsprop_med_summary$bsprop_025, y = 0, 
           xend = bsprop_med_summary$bsprop_025, yend = 1000,
           linetype = "dashed", color="blue") +
  annotate("segment", 
           x = bsprop_med_summary$bsprop_975, y = 0, 
           xend = bsprop_med_summary$bsprop_975, yend = 1000,
           linetype = "dashed", color="blue") +
  annotate("text", x = bsprop_med_summary$bsprop_025, y = 1200, label = "2.5th\npercentile") +
  annotate("text", x = bsprop_med_summary$bsprop_975, y = 1200, label = "97.5th\npercentile") +
  labs(
    x = "Bootstrapped proportion of surgical complications",
    y = "Count",
    title = "10,000 bootstrapped proportions (with 95% CI)"
  )

The original medical consultant data is bootstrapped 10,000 times. Each simulation creates a sample from the original data where the proportion of a complication is 3/62. The bootstrap 2.5 percentile proportion is 0 and the 97.5 percentile is 0.113. The result is: we are confident that, in the population, the true probability of a complication is between 0% and 11.3%.

The 95% Confidence Interval (CI): Using bootstrapping the true proportion of a complication is between 0 and 0.1129032. The interval overlaps the null statistic 0.10. There is a possibility that the consultant’s work is associated with a higher risk (\(p > 0.10\)), higher than the US average.

Case Study - Housing prices in Ames, Iowa (Mean)

In this section, we are using the ames data set which is available in the openintro package.

Data set contains information from the Ames Assessor’s Office used in computing assessed values for individual residential properties sold in Ames, IA from 2006 to 2010. Here, we are focusing on the sale price in USD.

prices <- ames %>% select(price)
mean_price <- mean(ames$price)
paste("Mean price = ",mean_price)

## [1] "Mean price =  180796.060068259"

The mean house sale price is

\[\bar{x} = 180796.06 \longrightarrow \text{Observed Statistic.}\]

Let \(\mu\) be the true mean price. We estimate \(\mu\) using the data, and label the estimate \(\bar{x}\).

Objective: Estimate the unknown population proportion by using the sample to approximate the mean house sale price.

set.seed(334422)

k <- 10000 # number of trials
n <- 100 # sample size

bs_prices <- rep_sample_n(prices,n,replace=TRUE,reps=k)

bs_prices_means <- bs_prices %>%
  group_by(replicate) %>% 
  summarize(stat = mean(price))

• Make a 95% confidence interval to give a plausible range of values for the population proportion.

probs_95 <- c(0.025,0.975)

ci_95 <- bs_prices_means %>%
  pull(stat) %>%
  quantile(probs_95)
  
ci_95

##     2.5%    97.5% 
## 165654.6 196754.8

The 95% confidence interval is (165777.5,197069.0). We are 95% confident that the true mean house sale price is in between $165777.5 and $197069.0.

ggplot(bs_prices_means, aes(x = stat)) +
  geom_histogram(bins = 50, fill = "gray") +
  annotate("segment", 
           x = ci_95[1], y = 0, 
           xend = ci_95[1], yend = 400,
           linetype = "dashed", color = "blue") +
  annotate("segment", 
           x = ci_95[2], y = 0, 
           xend = ci_95[2], yend = 400,
           linetype = "dashed", color = "blue") +
  annotate("text", x = ci_95[1], y = 402, label = "2.5th\npercentile") +
  annotate("text", x = ci_95[2], y = 402, label = "97.5th\npercentile") +
  labs(
    x = "Bootstrapped proportion of sale prices",
    y = "Count",
    title = "10,000 bootstrapped means  (with 95% CI)"
  )

• Make a 90% confidence interval to give a plausible range of values for the population proportion.

probs_90 <- c(0.05,0.95)

ci_90 <- bs_prices_means %>%
  pull(stat) %>%
  quantile(probs_90)
  
ci_90

##       5%      95% 
## 167904.5 194218.3

The 90% confidence interval is (168063.4,194202.3 ). We are 90% confident that the true mean house sale sprice is in between $168063.4 and $194202.3.

ggplot(bs_prices_means, aes(x = stat)) +
  geom_histogram(bins = 50, fill = "gray") +
  annotate("segment", 
           x = ci_90[1], y = 0, 
           xend = ci_90[1], yend = 400,
           linetype = "dashed", color = "blue") +
  annotate("segment", 
           x = ci_90[2], y = 0, 
           xend = ci_90[2], yend = 400,
           linetype = "dashed", color = "blue") +
  annotate("text", x = ci_90[1], y = 402, label = "5th\npercentile") +
  annotate("text", x = ci_90[2], y = 402, label = "95th\npercentile") +
  labs(
    x = "Bootstrapped proportion of sale prices",
    y = "Count",
    title = "10,000 bootstrapped means  (with 90% CI)"
  )

Case Study - Climate Change (Proportion)

This section is from the textbook OpenIntro: Section 15.4

A 2019 Pew Research report states the following:

To keep our computation simple, we will assume a total population size of 100,000 (even though that’s smaller than the population size of all US adults).

Roughly six-in-ten U.S. adults (62%) say climate change is currently affecting their local community either a great deal or some, according to a new Pew Research Center survey.

Source: Most Americans say climate change impacts their community, but effects vary by region

In this exercise, you will assume this 62% is a true population proportion and learn about how sample proportions can vary from sample to sample by taking smaller samples from the population. We will first create our population assuming a population size of 100,000. This means 62,000 (62%) of the adult population think climate change impacts their community, and the remaining 38,000 does not think so.

NOTE: Population parameters are unknown. This exercise pretends we know the popoulation proportion for demonstration puposes.

us_adults <- tibble(
  climate_change_affects = c(rep("Yes", 62000), rep("No", 38000))
)

ggplot(us_adults, aes(x = climate_change_affects)) +
  geom_bar() +
  labs(
    x = "", y = "",
    title = "Do you think climate change is affecting your local community?"
  ) +
  coord_flip() 

us_adults %>%
  count(climate_change_affects) %>%
  mutate(p = n /sum(n))

## # A tibble: 2 × 3
##   climate_change_affects     n     p
##   <chr>                  <int> <dbl>
## 1 No                     38000  0.38
## 2 Yes                    62000  0.62

In this section, you’ll start with a simple random sample of size 60 from the population.

n <- 60
samp <- us_adults %>%
  sample_n(size = n)

In essence, bootstrapping assumes that there are more of observations in the populations like the ones in the observed sample. So we “reconstruct” the population by resampling from our sample, with replacement. The bootstrapping scheme is as follows:

• Step 1. Take a bootstrap sample - a random sample taken with replacement from the original sample.

• Step 2. Calculate the bootstrap statistic (e.g. proportion, mean, etc.) computed on the bootstrap samples.

• Step 3. Repeat steps (1) and (2) many times to create a bootstrap distribution - a distribution of bootstrap statistics.

• Step 4. Calculate the bounds of the XX% confidence interval as the middle XX% of the bootstrap distribution.

The two methods to compute the confidence interval is using the percentile method or using the standard error method.

Instead of coding up each of these steps, we will construct confidence intervals using the infer package.

This code will find the 95 percent confidence interval for proportion of US adults who think climate change affects their local community.

samp %>%
  specify(response = climate_change_affects, success = "Yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_ci(level = 0.95, type = "percentile")

## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1    0.516     0.75

• In specify we specify the response variable and the level of that variable we are calling a success.

• In generate we provide the number of resamples we want from the population in the reps argument (this should be a reasonably large number) as well as the type of resampling we want to do, which is “bootstrap” in the case of constructing a confidence interval.

• Then, we calculate the sample statistic of interest for each of these resamples, which is proportion.

• To compute the confidence interval, we use the get_ci function and use 0.95 as the confidence level. Here, wer are using the percentile method.

Feel free to test out the rest of the arguments for these functions, since these commands will be used together to calculate confidence intervals and solve inference problems for the rest of the semester.

To recap: even though we don’t know what the full population looks like, we’re 95% confident that the true proportion of US adults who think climate change affects their local community is between the two bounds reported as result of this pipeline.

Earlier in this section, we assumed that 62% is a true population proportion. In our 95% interval, the 62% is within this interval, which means that we captured the true population proportion. Again, population parameters are unknown, this example is for demonstration purposes. In reality, There is still a 5% chance that the interval is wrong.

Lab Exercises

Note: You must include your code and results to answer the exercise problems.

I. Climate Change

From the Case Study - Climate Change (Proportion) section, use similar codes to answer the following questions. These problems requre a combination of what we have learned in lab since the first week (e.g. for loops, while loops, conditional statements etc).

For convenience, below is the R code used in this section.

# "the population" or "the data"
us_adults <- tibble(
  climate_change_affects = c(rep("Yes", 62000), rep("No", 38000))
)

1. Obtain a random sample of 60 from the us_adults data. Calculate the sample proportion, and use these to calculate and store the lower and upper bounds of the 95% confidence interval (use 1000 trials or resamples). Repeat these steps 100 times.

2. From problem 1, you should have 100 95% confidence intervals. Compute the number of intervals that contain the true population proportion. Remember that we assumed that the true population proportion is 0.62 (62%).

3. A 95% confidence interval gives us a region where, had we redone the same data, then 95% of the time, the true value \(p\) will be contained in the interval. What proportion of your confidence intervals include the true population proportion? Is this proportion exactly equal or close to the confidence level? Explain your reasoning.

4. Choose a different confidence level than 95%. Would you expect a confidence interval at this level to be wider or narrower than the confidence interval you calculated at the 95% confidence level? Explain your reasoning.

5. Find a confidence interval for the proportion of US Adults who think climate change is affecting their local community with a confidence level of your choosing (other than 95%) and interpret it.