library(tidyverse)
library(infer)

gss_survey <- read_csv("http://www.cknuckles.com/statsprogramming/data/gss_survey_data.csv")

# Example - Bootstrap CI for a Mean - Using the Infer Package
########################################################################################

workhours_2018 <- gss_survey |>
  filter(year == "2018") |>
  select(workhours) |>
  drop_na(workhours)

bootstrap_distribution <- workhours_2018 |>
  specify(response = workhours) |>
  generate(reps = 10000, type = "bootstrap") |>
  calculate(stat = "mean")
bootstrap_distribution

# The same bootstrap distribution computed without Infer, as in the previous scripts:
# bootstrap_dist_dplyr <- workhours_2018 |>
#   rep_slice_sample(n=828, reps=10000, replace=TRUE) |>
#   group_by(replicate) |>
#   summarize( resample_mean = mean(workhours) )

bootstrap_distribution |> visualize()

bootstrap_distribution |>
  get_ci(level=.95, type = "percentile")

bootstrap_distribution |>
  get_ci(level=.95, type = "se", point_estimate=41.53)

# Calculate original sample mean (by leaving out call to generate() bootstrap samples)
workhours_2018 |>
  specify(response = workhours) |>
  calculate(stat = "mean")

bootstrap_distribution |>
  visualize() +                            # Histogram only
  shade_ci(endpoints = c(40.549,42.511))   # Add shaded plot layer for CI

# Store the get_ci() results into a variable
percentile_ci <- bootstrap_distribution |>
  get_ci(level=.95, type = "percentile")
# Send the stored get_ci() results to shade_ci()
bootstrap_distribution |>
  visualize() +
  shade_ci(endpoints = percentile_ci)      # percentile_ci contains the CI endpoints


# Example - Bootstrap CI for a Proportion - Using the Infer Package
########################################################################################

tvhours_2018 <- gss_survey |>
  filter(year == "2018") |>
  select(tvhours) |>
  drop_na(tvhours)

# Can't use logical expression to define success:
# specify(response = tvhours, success = "tvhours >= 3")    # OOPS - Error
# So first add a TRUE/FALSE column
tvhours_2018_with_potatoes <- tvhours_2018 |>
  mutate(is_potato = (tvhours>=3))  # new column is_potato = TRUE/FALSE
tvhours_2018_with_potatoes

set.seed(5)   # Same random results as in the book

# Compute the single-proportion bootstrap distribution
bootstrap_distribution <- tvhours_2018_with_potatoes |>
  specify(response = is_potato, success = "TRUE") |>
  generate(reps = 10000, type = "bootstrap") |>
  calculate(stat = "prop")

# Compute the SE method bootstrap CI
se_ci <- bootstrap_distribution |>
  get_ci(level=.99, type = "se", point_estimate=.340)
se_ci

# Visualize the results
bootstrap_distribution |>
  visualize() +
  shade_ci(endpoints = se_ci)

# Calculate original sample proportion (by leaving out call to generate() bootstrap samples)
tvhours_2018_with_potatoes |>
  specify(response = is_potato, success = "TRUE") |>
  calculate(stat = "prop")


# Example - Bootstrap CI for a Difference of Proportions - Using the Infer Package
########################################################################################

gss_survey_2018_wrangled <- gss_survey |>
  filter(year == "2018") |>
  drop_na(tvhours,status) |>
  mutate(is_potato = (tvhours>=3)) |>
  mutate(married_status = case_when(
                              (status == "Married") ~ "married",
                              (status != "Married") ~ "single"
                          )
  ) |>
  select(tvhours,is_potato,status,married_status)
gss_survey_2018_wrangled

bootstrap_distribution <- gss_survey_2018_wrangled |>
  specify(response = is_potato, success = "TRUE", explanatory = married_status) |>
  generate(reps = 10000, type = "bootstrap") |>
  calculate(stat = "diff in props", order = c("single","married"))

bootstrap_distribution |>
  get_ci(level=.95, type = "percentile")

bootstrap_distribution |>
  get_ci(level=.95, type = "se", point_estimate=(.383-.291))

# Equivalent specify() calls
# specify(formula = is_potato ~ married_status, success = "TRUE")
# specify(response = is_potato, success = "TRUE", explanatory = married_status)

# Template for a difference of means: response var1 is numerical, explanatory var2 has two categories
# bootstrap_distribution <- data_table |>
#   specify(response = var1, explanatory = var2) |>
#   generate(reps = 10000, type = "bootstrap") |>
#   calculate(stat = "diff in means", order = c("categoryA","categoryB"))