illustrating-cis-for-percentages.R

# load packages
library(tidyverse)

# change these values
n_1s <- 10  # number of 1s in the box
n_0s <- 10  # number of 0s in the box
sample_size <- 25  # number of times to draw from the box (w/ replacement)

# number of times to repeat the study 
# note: in practice, we do this only once
n_sims <- 2500

# the code below performs the hypothetical study of drawing 'sample_size' times
# from the 'box' (with replacement) and calculating the percentage of ones and 
# using the bootstrap method to calculate the SE and 95% confidence interval.
df <- NULL
box <- c(rep(1, n_1s), rep(0, n_0s))
for (i in 1:n_sims) {
  sample <- sample(box, size = sample_size, replace = TRUE)
  p <- 100*mean(sample)
  frac_1s <- mean(sample)
  se_number <- sqrt(sample_size)*sqrt(frac_1s*(1-frac_1s))
  se_percent <- 100*se_number/sample_size
  lwr <- p - 2*se_percent
  upr <- p + 2*se_percent
  df0 <- data.frame(iter = i, p, lwr, upr)
  df <- rbind(df, df0)
}
true_percent <- 100*mean(box)
df <- mutate(df, covers = ifelse(lwr <= true_percent & upr >= true_percent, "Yes", "No"))
message <- paste0("In the ", n_sims, " studies, ", round(100*mean(df$covers ==  "Yes")), "% cover the true value of ", round(100*mean(box)), "%.")
ggplot(filter(df, iter <= 100), aes(x = p, y = iter, xmin = lwr, xmax = upr, color = covers)) + 
  geom_vline(xintercept = true_percent) + 
  geom_errorbarh(height = 0) + 
  geom_point() + 
  labs(x = "Estimated Percent of Ones in the Box",
       y = "Study Number",
       title = "The Confidence Intervals from 100 Studies",
       caption = message)