R-Lab Categorical Data

Load Packages

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)

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?

unique(yrbss$text_while_driving_30d)

[1] "0"             NA              "30"            "did not drive"
[5] "1-2"           "3-5"           "20-29"         "10-19"        
[9] "6-9"          

yrbss |>
  group_by(text_while_driving_30d) |>
  count()

# A tibble: 9 × 2
# Groups:   text_while_driving_30d [9]
  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

Answer 1: 4792 students claim that they did not text and drive within the past 30 days. 925 students texted 1-2 days, 493 students texted 3-5 days, 311 students texted 6-9 days, 373 students texted 10-19 days, and 298 students texted 20-29 days. 827 students texted and drove all of the past 30 days while 4646 students did not drive at all.

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") |>
  mutate(text_ind = ifelse (text_while_driving_30d == "30", "yes", "no")) |>
  filter(!is.na(text_ind))
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    15 female 9     hispanic Native Hawaiian or Other…   1.73   84.4 never     
3    15 female 9     not      Black or African American   1.6    55.8 never     
4    16 male   9     not      Black or African American   1.68   74.8 never     
5    14 male   9     not      Black or African American   1.73   73.5 never     
6    15 male   9     not      Black or African American   1.83   67.6 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>

no_helmet |>
  summarise(prop_text_no_helmet = mean(text_ind == "yes")) |>
  pull()

[1] 0.07119791

Answer 2: The propotion of people who have texted while driving every day in the past 30 days and never wear helmets is 0.071.

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?

boot <- no_helmet |>
  specify(response = text_ind, success = "yes") |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "prop")
boot

Response: text_ind (factor)
# A tibble: 1,000 × 2
   replicate   stat
       <int>  <dbl>
 1         1 0.0737
 2         2 0.0727
 3         3 0.0677
 4         4 0.0717
 5         5 0.0723
 6         6 0.0681
 7         7 0.0686
 8         8 0.0718
 9         9 0.0697
10        10 0.0684
# ℹ 990 more rows

ci <- get_ci(boot, level = .95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1   0.0652   0.0777

se <- boot |>
  summarize(se = sd(stat)) |>
  pull()
se

[1] 0.003130602

me <- 1.96*se
me

[1] 0.00613598

Answer 3: The margin of error is 0.006.

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.

yrbss

# A tibble: 13,583 × 13
     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 Americ…  NA      NA   never     
 2    14 female 9     not      Black or African Americ…  NA      NA   never     
 3    15 female 9     hispanic Native Hawaiian or Othe…   1.73   84.4 never     
 4    15 female 9     not      Black or African Americ…   1.6    55.8 never     
 5    15 female 9     not      Black or African Americ…   1.5    46.7 did not r…
 6    15 female 9     not      Black or African Americ…   1.57   67.1 did not r…
 7    15 female 9     not      Black or African Americ…   1.65  132.  did not r…
 8    14 male   9     not      Black or African Americ…   1.88   71.2 never     
 9    15 male   9     not      Black or African Americ…   1.75   63.5 never     
10    15 male   10    not      Black or African Americ…   1.37   97.1 did not r…
# ℹ 13,573 more rows
# ℹ 5 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>

yrbss |>
  group_by(gender) |>
count()

# A tibble: 3 × 2
# Groups:   gender [3]
  gender     n
  <chr>  <int>
1 female  6621
2 male    6950
3 <NA>      12

only_hispanic <- yrbss |>
  filter(hispanic == "hispanic") |>
  mutate(female_hispanic = ifelse(gender == "female", "yes", "no")) |>
  filter(!is.na(female_hispanic))
head(only_hispanic)

# A tibble: 6 × 14
    age gender grade hispanic race                      height weight helmet_12m
  <int> <chr>  <chr> <chr>    <chr>                      <dbl>  <dbl> <chr>     
1    15 female 9     hispanic Native Hawaiian or Other…   1.73   84.4 never     
2    12 female 9     hispanic Black or African American   1.8    90.7 did not r…
3    17 female 11    hispanic <NA>                        1.68   86.2 did not r…
4    12 male   12    hispanic White                      NA      NA   never     
5    18 male   12    hispanic <NA>                        1.65   72.6 never     
6    15 male   10    hispanic Black or African American   1.78   68.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>, female_hispanic <chr>

only_hispanic |>
  summarise(prop_female_hispanic = mean(female_hispanic == "yes")) |>
  pull()

[1] 0.4906542

boot <- only_hispanic |>
  specify(response = female_hispanic, success = "yes") |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "prop")
boot

Response: female_hispanic (factor)
# A tibble: 1,000 × 2
   replicate  stat
       <int> <dbl>
 1         1 0.471
 2         2 0.492
 3         3 0.492
 4         4 0.500
 5         5 0.498
 6         6 0.491
 7         7 0.498
 8         8 0.495
 9         9 0.502
10        10 0.508
# ℹ 990 more rows

ci <- get_ci(boot, level = .95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1    0.474    0.507

Answer 4 part i): The confidence intervals for being both female and hispanic are 0.474 and 0.507.

se <- boot |>
  summarize(se = sd(stat)) |>
  pull()
se

[1] 0.008243454

me <- 1.96*se
me

[1] 0.01615717

Answer 4 part ii): The margin of error is 0.017.

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?

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")

Answer 5: The margin of error seems to increase as the population proportion increases. However, once the population proportion reaches 50 % (maximum margin of error) the margin of error starts to decrease.