library(tidyverse)
library(infer)   # Contains rep_slice_sample() function

# Example - Bootstrap CIs for the Proportion of Couch Potatoes (tvhours >= 3)
########################################################################################

gss_survey <- read_csv("http://www.cknuckles.com/statsprogramming/data/gss_survey_data.csv")
tvhours_2018 <- gss_survey |>
  filter(year == "2018") |>
  select(tvhours) |>
  drop_na(tvhours)

# Bootstrap Distribution of 10,000 Proportions
bootstrap_distribution_of_proportions <- tvhours_2018 |>
  rep_slice_sample(n=833 , reps=10000 ,  replace=TRUE) |>
  group_by(replicate) |>
  summarize(
    num_couch_potatoes = sum(tvhours >= 3),
    n = n(),
    resample_proportion = num_couch_potatoes/n   # proportion of couch potatoes
  )

ggplot(bootstrap_distribution_of_proportions, mapping = aes(x = resample_proportion)) +
  geom_histogram(bins=52, color="white") +
  labs(x = "Bootstrap Sample Proportions")

# Bootstrap Statistics from the Bootstrap Distribution
bootstrap_distribution_of_proportions |>
  summarize(
      boot_dist_mean = mean(resample_proportion),
      boot_dist_SE = sd(resample_proportion)
  )

# Bootstrap SE method 99% CI
SE <- .0164
z_alpha2 <- 2.58   # 2.58 = qnorm(.995, mean=0, sd=1)
c(.340 - z_alpha2*SE, .340 + z_alpha2*SE)

# Bootstrap Percentile method 99% CI
quantile(bootstrap_distribution_of_proportions$resample_proportion, c(.005,.995))

summary(bootstrap_distribution_of_proportions$resample_proportion)