library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata
library(infer)
data(yrbss, package = "openintro")

Exercise 1

What are the counts within each category for the amount of days these students have texted while driving within the past 30 days?

texting_counts <- yrbss %>%
  count(text_while_driving_30d)
print(texting_counts)
## # A tibble: 9 × 2
##   text_while_driving_30d     n
##   <chr>                  <int>
## 1 0                       4792
## 2 1-2                      925
## 3 10-19                    373
## 4 20-29                    298
## 5 3-5                      493
## 6 30                       827
## 7 6-9                      311
## 8 did not drive           4646
## 9 <NA>                     918

Exercise 2

What is the proportion of people who have texted while driving every day in the past 30 days and never wear helmets?

no_helmet <- yrbss %>%
  filter(helmet_12m == "never")

no_helmet <- no_helmet %>%
  mutate(text_ind = ifelse(text_while_driving_30d == "30", "yes", "no"))

head(no_helmet)
## # A tibble: 6 × 14
##     age gender grade hispanic race                      height weight helmet_12m
##   <int> <chr>  <chr> <chr>    <chr>                      <dbl>  <dbl> <chr>     
## 1    14 female 9     not      Black or African American  NA      NA   never     
## 2    14 female 9     not      Black or African American  NA      NA   never     
## 3    15 female 9     hispanic Native Hawaiian or Other…   1.73   84.4 never     
## 4    15 female 9     not      Black or African American   1.6    55.8 never     
## 5    14 male   9     not      Black or African American   1.88   71.2 never     
## 6    15 male   9     not      Black or African American   1.75   63.5 never     
## # ℹ 6 more variables: text_while_driving_30d <chr>, physically_active_7d <int>,
## #   hours_tv_per_school_day <chr>, strength_training_7d <int>,
## #   school_night_hours_sleep <chr>, text_ind <chr>
total_no_helmet <- nrow(no_helmet)


text_every_day_no_helmet <- no_helmet %>%
  filter(text_ind == "30") %>%
  nrow()

proportion_text_every_day_no_helmet <- text_every_day_no_helmet / total_no_helmet
print(proportion_text_every_day_no_helmet)
## [1] 0

Exercise 3

What is the margin of error for the estimate of the proportion of non-helmet wearers that have texted while driving each day for the past 30 days based on this survey?

no_helmet <- yrbss %>%
  filter(helmet_12m == "never")

no_helmet <- no_helmet %>%
  mutate(text_ind = ifelse(text_while_driving_30d == "30", "yes", "no"))

total_no_helmet <- nrow(no_helmet)

text_every_day_no_helmet <- no_helmet %>%
  filter(text_ind == "yes") %>%
  nrow()

p_hat <- text_every_day_no_helmet / total_no_helmet

z <- 1.96

moe <- z * sqrt((p_hat * (1 - p_hat)) / total_no_helmet)

print(moe)
## [1] 0.005840733

Exercise 4

Using the infer package, calculate confidence intervals for two other categorical variables (you’ll need to decide which level to call “success”, and report the associated margins of error. Interpet the interval in context of the data. It may be helpful to create new data sets for each of the two countries first, and then use these data sets to construct the confidence intervals.

helmet_success <- "never"
not_hispanic <- "not"

helmet_filtered <- yrbss %>%
  filter(helmet_12m %in% c("never", "rarely"))

helmet_ci <- helmet_filtered %>%
  specify(response = helmet_12m, success = helmet_success) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_ci(level = 0.95)

print(helmet_ci)
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1    0.900    0.914

How does the proportion affect the margin of error?

n <- 1000

p <- seq(from = 0, to = 1, by = 0.01)
me <- 2 * sqrt(p * (1 - p)/n)

dd <- data.frame(p = p, me = me)
ggplot(data = dd, aes(x = p, y = me)) + 
  geom_line() +
  labs(x = "Population Proportion", y = "Margin of Error")

Exercise 5

Describe the relationship between p and me. Include the margin of error vs. population proportion plot you constructed in your answer. For a given sample size, for which value of p is margin of error maximized?

# The X axis of the plot represents the population proportion (p) from 0 to 1, the Y axis represents the margin error (me). The curve is symetric and reaches its peak when p is 0.5, showing that the margin error is maximized at that point, and when p moves away from 0.5 towards 0 or 1 the margin error decreases. This shows that the margin error is dependant on the sample size and on the porportion.

Exercise 6

Describe the sampling distribution of sample proportions at n=300 and p=0.1. Be sure to note the center, spread, and shape.

# Center: The mean of the sampling distribution is 0.1

# Spread: The standar error of the sampling is approximately 0.01732

# Shape: The sampling distribution can be summed to be nrealy normal because both values are greater than 10
# np = 300 * 0.1 = 30
# n(1-p) = 300 * 0.9 = 270

Exercise 7

Keep n constant and change p. How does the shape, center, and spread of the sampling distribution vary as p changes. You might want to adjust min and max for the x-axis for a better view of the distribution.

library(ggplot2)

n <- 300

p_values <- c(0.1, 0.3, 0.5, 0.7, 0.9)

df_list <- list()

for (p in p_values) {
  mean_p <- p
  se_p <- sqrt(p * (1 - p) / n)
  x <- seq(mean_p - 4 * se_p, mean_p + 4 * se_p, length.out = 100)
  y <- dnorm(x, mean = mean_p, sd = se_p)
  df <- data.frame(x = x, y = y, p = as.factor(p))
  df_list[[as.character(p)]] <- df
}

df_all <- do.call(rbind, df_list)

ggplot(df_all, aes(x = x, y = y, color = p)) +
  geom_line() +
  labs(title = "Different Values of p",
       x = "Sample Proportion",
       y = "Density",
       color = "Population Proportion (p)") +
  theme_minimal()

# Center: The mean of the sampling distribution is the value of p


# Spread: Looking at the results of the SE, as p changes the value of p(1-p) changes; when p es close to 0 or 1, the value of p(1-p) is small resulting is a smaller SE and a more narrow distribution. When p is close to 0.5, the value of p(1-p) is closer to 0.25, resulting in a larger SE and a wider distribution

# Shape: The shape of the sampling distribution remains aproximately normal as long as the success-failure contidition is met, this condition ensures that the binomial distribution can be sproximated by the normal distribution

Exercise 8

library(ggplot2)

create_df <- function(n, p) {
  mean_p <- p
  se_p <- sqrt(p * (1 - p) / n)
  x <- seq(mean_p - 4 * se_p, mean_p + 4 * se_p, length.out = 100)
  y <- dnorm(x, mean = mean_p, sd = se_p)
  df <- data.frame(x = x, y = y, n = as.factor(n))
  return(df)
}

n_values <- c(100, 300, 1000)
p <- 0.5

df_list <- list()

for (n in n_values) {
  df <- create_df(n, p)
  df_list[[as.character(n)]] <- df
}

df_all <- do.call(rbind, df_list)

ggplot(df_all, aes(x = x, y = y, color = n)) +
  geom_line() +
  labs(title = "Sampling Distributions for Different Sample Sizes (n) with p = 0.5",
       x = "Sample Proportion",
       y = "Density",
       color = "Sample Size (n)") +
  theme_minimal()

# Center: The mean of the sampling distribution is 0.5 

# Spread: The spread of the distribution decreases as n increases leading to a narrower distribution, and as n decreases the SE increases resulting in a wider distribution.

# Shape: The shape of the sampling distribution is assumed to be normal if the succes-failure condition is met.