Inference for Categorical Data Lab

Author

Alexandra Veremeychik

Based on a lab provided by the OpenIntro Introduction to Modern Statistics textbook.

Getting Started

Load packages

library(tidyverse)
library(openintro)
library(infer)

The data

data("yrbss")
head(yrbss)

# A tibble: 6 × 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 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    15 female 9     not      Black or African American   1.5    46.7 did not r…
6    15 female 9     not      Black or African American   1.57   67.1 did not r…
# ℹ 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

starting httpd help server ... done

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

text_drive <- yrbss |>
  group_by(text_while_driving_30d) |>
  count()
text_drive

# 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

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

# Code provided by directions
no_helmet <- yrbss |>
  filter(helmet_12m == "never") |>
  mutate(text_ind = ifelse(text_while_driving_30d == "30", "yes", "no")) |>
  filter(!is.na(text_ind))

# My way
text_no_helmet <- yrbss |>
  filter(helmet_12m == "never" & text_while_driving_30d == "30")
  
text_no_helmet_prop <- nrow(text_no_helmet) / nrow(yrbss)
text_no_helmet_prop

[1] 0.03408673

Inference on proportions

Bootstrap distribution

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.0755
 2         2 0.0735
 3         3 0.0706
 4         4 0.0729
 5         5 0.0667
 6         6 0.0714
 7         7 0.0729
 8         8 0.0697
 9         9 0.0703
10        10 0.0721
# ℹ 990 more rows

Confidence interval

ci <- get_ci(boot, level = 0.95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1   0.0649   0.0770

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?

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

me <- 1.96 * se
me

[1] 0.006099682

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. Interpret 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.

What proportion of female-identifying people did not wear a helmet when biking in the past 12 months?

female_biked <- yrbss |>
  filter(gender == "female" & helmet_12m != "did not ride") # Filter out non-riders

female_no_helmet <- female_biked |>
  mutate(helmet_y_n = ifelse(helmet_12m == "never", "no", "yes"))

boot <- female_no_helmet |>
  specify(response = helmet_y_n, success = "no") |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "prop")

ci <- get_ci(boot, level = 0.95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1    0.770    0.797

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

me <- 1.96 * se
me

[1] 0.01355729

We are 95% confident that the true proportion of female-identifying people who did not wear a helmet when biking in the past 12 months is between 0.771 and 0.797.

What about male-identifying people?

male_biked <- yrbss |>
  filter(gender == "male" & helmet_12m != "did not ride") # Filter out non-riders

male_no_helmet <- male_biked |>
  mutate(helmet_y_n = ifelse(helmet_12m == "never", "no", "yes"))

boot <- male_no_helmet |>
  specify(response = helmet_y_n, success = "no") |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "prop")

ci <- get_ci(boot, level = 0.95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1    0.801    0.822

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

me <- 1.96 * se
me

[1] 0.0107105

We are 95% confident that the true proportion of male-identifying people who did not wear a helmet when biking in the past 12 months is between 0.801 and 0.830.

unique(yrbss$school_night_hours_sleep)

[1] "8"   "6"   "<5"  "9"   "10+" "7"   "5"   NA

What proportion of people report typically getting at least 8 hours of sleep?

sleep_8 <- yrbss |>
  mutate(at_least_8_sleep = ifelse(school_night_hours_sleep == "8"
                                  |school_night_hours_sleep == "9"
                                  |school_night_hours_sleep == "10+",
                                  "yes", "no"))

at_least_8_sleep_prop <-sleep_8 |>
  filter(!is.na(at_least_8_sleep)) |>
  summarise(prop = mean(at_least_8_sleep == "yes")) |>
  pull()
at_least_8_sleep_prop

[1] 0.3057154

Of the people who report getting at least 8 hours of sleep, what proportion report being physically active for 60+ minutes in the last 7 days?

sleep_8_active <- sleep_8 |>
  filter(!is.na(physically_active_7d)) |>
  mutate(active_y_n = ifelse(physically_active_7d > 0, "yes", "no"))

sleep_8_active_prop <- sleep_8_active |>
  summarise(prop = mean(active_y_n == "yes")) |>
  pull()
sleep_8_active_prop

[1] 0.8368144

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

ci <- get_ci(boot, level = 0.95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1    0.831    0.843

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

[1] 0.003161668

We are 95% confident that of the population that gets at least 8 hours of sleep on a school night, the true proportion of people who were physically active for 60+ minutes in the last 7 days is between 0.830 and 0.843.

How does the proportion affect the margin of error?

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

dd <- data.frame(p = p, me = me)
ggplot(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 closer the population proportion is to 0 or 1, the lower the margin of error. Very small and very large population proportions result in a low margin of error. The margin of error is highest around the 0.50 population proportion

Success-failure condition

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.

The sampling distribution is centered on 0.010 with a mean of 0.099 and a standard error of 0.017. It is symmetrical with a range of approximately 0.10

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.

When the proportion is very small, such as when the number of successes is 9, the sampling distribution is right skewed.

When the proportion is very large, such as when the number of successes is 291, the sampling distribution is left skewed.

EXERCISE 8: Now also change n. How does n appear to affect the distribution of p^?

The smaller n is, the more skewed the sampling distribution.

More Practice

Exercise 9: Is there convincing evidence that those who sleep 10+ hours per day are more likely to strength train every day of the week? As always, write out the hypotheses for any tests you conduct and outline the status of the conditions for inference. If you find a significant difference, also quantify this difference with a confidence interval.

H₀: p_sleep 10+ = P_sleep <10

There is no relationship between sleeping at least 10 hours per day and strength training every day of the week.

H_A: p_sleep 10+ > P_sleep <10

Those who sleep at least 10 hours per day are more likely to strength train every day of the week.

sleep_strength <- yrbss |>
  mutate(sleep_10_y_n = ifelse(school_night_hours_sleep == "10+", "yes", "no")) |>
  mutate(strength_y_n = ifelse(strength_training_7d == 7, "yes", "no")) |>
  filter(!is.na(sleep_10_y_n)) |>
  filter(!is.na(strength_y_n))

# Get p_hat
sleep_10_y_prop <- sleep_strength |>
  summarise(prop = mean(sleep_10_y_n == "yes")) |>
  pull()

strength_y_prop <- sleep_strength |>
  summarise(prop = mean(strength_y_n == "yes")) |>
  pull()

# Check for S-F condition, np >= 10 and n(1-p) >= 10
nrow(sleep_strength) * sleep_10_y_prop >= 10

[1] TRUE

nrow(sleep_strength) * (1 - sleep_10_y_prop) >= 10

[1] TRUE

nrow(sleep_strength) * strength_y_prop >= 10

[1] TRUE

nrow(sleep_strength) * (1 - strength_y_prop) >= 10

[1] TRUE

# Bootsrap distributions
boot <- sleep_strength |>
  specify(sleep_10_y_n ~ strength_y_n, success = "yes") |>
  generate(reps = 1000, type = "bootstrap") |>
  calculate(stat = "diff in props", order = c("yes", "no"))

# Confidence interval
ci <- get_ci(boot, level = 0.95)
ci

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1  0.00952   0.0285

# Strandard error
se <- boot |>
  summarise(sd(stat)) |>
  pull()
se

[1] 0.004837461

Exercise 10: Let’s say there has been no difference in likeliness to strength train every day of the week for those who sleep 10+ hours. What is the probablity that you could detect a change (at a significance level of 0.05) simply by chance? Hint: Review the definition of the Type 1 error.

Exercise 11: Suppose you’re hired by the local government to estimate the proportion of residents that attend a religious service on a weekly basis. According to the guidelines, the estimate must have a margin of error no greater than 1% with 95% confidence. You have no idea what to expect for p. How many people would you have to sample to ensure that you are within the guidelines?

Hint: Refer to your plot of the relationship between p and margin of error. This question does not require using a dataset.

0.01 > (1.96)sqrt(p(1-P)/n)

0.0196 > sqrt(p(1-P)/n)

0.00038416 > p(1-p)/n

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

dd <- data.frame(p = p, me = me)
ggplot(dd, aes(x = p, y = me)) +
  geom_rect(xmin = 0, xmax = 1, ymin = 0, ymax = 0.01, alpha = 0.1, fill = "lightblue") +
  geom_line(size = 1) +
  labs(x = "Population Proportion", y = "Margin of Error") +
  geom_hline(yintercept = 0.01)

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.