Getting Started

Load packages

In this lab, we will explore and visualize the data using the tidyverse suite of packages, and perform statistical inference using infer. The data can be found in the companion package for OpenIntro resources, openintro.

Let’s load the packages.

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

The data

You will be analyzing the same dataset as in the previous lab, where you delved into a sample from the Youth Risk Behavior Surveillance System (YRBSS) survey, which uses data from high schoolers to help discover health patterns. The dataset is called yrbss.

  1. What are the counts within each category for the amount of days these students have texted while driving within the past 30 days?
head(yrbss)
## # A tibble: 6 × 13
##     age gender grade hispa…¹ race  height weight helme…² text_…³ physi…⁴ hours…⁵
##   <int> <chr>  <chr> <chr>   <chr>  <dbl>  <dbl> <chr>   <chr>     <int> <chr>  
## 1    14 female 9     not     Blac…  NA      NA   never   0             4 5+     
## 2    14 female 9     not     Blac…  NA      NA   never   <NA>          2 5+     
## 3    15 female 9     hispan… Nati…   1.73   84.4 never   30            7 5+     
## 4    15 female 9     not     Blac…   1.6    55.8 never   0             0 2      
## 5    15 female 9     not     Blac…   1.5    46.7 did no… did no…       2 3      
## 6    15 female 9     not     Blac…   1.57   67.1 did no… did no…       1 5+     
## # … with 2 more variables: strength_training_7d <int>,
## #   school_night_hours_sleep <chr>, and abbreviated variable names ¹​hispanic,
## #   ²​helmet_12m, ³​text_while_driving_30d, ⁴​physically_active_7d,
## #   ⁵​hours_tv_per_school_day
yrbss %>%
  count(text_while_driving_30d, sort=TRUE) 
## # A tibble: 9 × 2
##   text_while_driving_30d     n
##   <chr>                  <int>
## 1 0                       4792
## 2 did not drive           4646
## 3 1-2                      925
## 4 <NA>                     918
## 5 30                       827
## 6 3-5                      493
## 7 10-19                    373
## 8 6-9                      311
## 9 20-29                    298
  1. What is the proportion of people who have texted while driving every day in the past 30 days and never wear helmets?
no_helmet_everyday_driving <- yrbss %>%
  select(text_while_driving_30d, helmet_12m) %>%
  filter(!is.na(text_while_driving_30d)) %>%
  filter(helmet_12m == "never") %>%
  mutate(text_everyday_never_helmet = ifelse(text_while_driving_30d == "30", "Yes", "No")) %>%
  group_by(text_everyday_never_helmet) %>%
  summarise(count=n())%>%
  mutate(props= round(count/sum(count), 2))
no_helmet_everyday_driving
## # A tibble: 2 × 3
##   text_everyday_never_helmet count props
##   <chr>                      <int> <dbl>
## 1 No                          6040  0.93
## 2 Yes                          463  0.07

There are 7% of people who have texted while driving everyday in the past 30 days and never wear helmets.

Remember that you can use filter to limit the dataset to just non-helmet wearers. Here, we will name the dataset no_helmet.

data('yrbss', package='openintro')
no_helmet <- yrbss %>%
  filter(helmet_12m == "never")

Also, it may be easier to calculate the proportion if you create a new variable that specifies whether the individual has texted every day while driving over the past 30 days or not. We will call this variable text_ind.

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

Inference on proportions

When summarizing the YRBSS, the Centers for Disease Control and Prevention seeks insight into the population parameters. To do this, you can answer the question, “What proportion of people in your sample reported that they have texted while driving each day for the past 30 days?” with a statistic; while the question “What proportion of people on earth have texted while driving each day for the past 30 days?” is answered with an estimate of the parameter.

The inferential tools for estimating population proportion are analogous to those used for means in the last chapter: the confidence interval and the hypothesis test.

no_helmet <- no_helmet %>%
  filter(!is.na(text_ind))
set.seed(1234)
no_helmet %>%
  specify(response = text_ind, success = "yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_ci(level = 0.95)
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   0.0652   0.0777

Note that since the goal is to construct an interval estimate for a proportion, it’s necessary to both include the success argument within specify, which accounts for the proportion of non-helmet wearers than have consistently texted while driving the past 30 days, in this example, and that stat within calculate is here “prop”, signaling that you are trying to do some sort of inference on a proportion.

  1. 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 <- no_helmet %>%
  drop_na(text_ind)

n <- count(no_helmet)
p <- count(no_helmet %>% filter(text_ind == "yes"))/n
SE <- sqrt((p*(1-p))/n)

ME <- (1.96*SE)*100

print(ME)
##           n
## 1 0.6250207

There’s a margin of error of 0.62% of non-helmet wearers that have texted while driving each day for the past 30 days.

  1. 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.
yrbss %>%
  count(hours_tv_per_school_day, sort=TRUE)
## # A tibble: 8 × 2
##   hours_tv_per_school_day     n
##   <chr>                   <int>
## 1 2                        2705
## 2 <1                       2168
## 3 3                        2139
## 4 do not watch             1840
## 5 1                        1750
## 6 5+                       1595
## 7 4                        1048
## 8 <NA>                      338
no_watch<- yrbss %>%
  select(hours_tv_per_school_day, school_night_hours_sleep) %>%
  filter(!is.na(hours_tv_per_school_day)) %>%
  filter(!is.na(school_night_hours_sleep)) %>%
  filter(hours_tv_per_school_day == "do not watch") %>%
  mutate(not_watch_tv = ifelse(school_night_hours_sleep >= "8", "Yes", "No")) %>%
  group_by(not_watch_tv) %>%
  summarise(count=n())%>%
  mutate(props= round(count/sum(count), 2))
no_watch
## # A tibble: 2 × 3
##   not_watch_tv count props
##   <chr>        <int> <dbl>
## 1 No            1265  0.76
## 2 Yes            410  0.24

Approximately 24% of students slept 8 or more hours a day and did not watch tv.

I wanted to examine how many people slept 8 or more hours that did not watch tv:

set.seed(1234)
not_watch_tv <- yrbss %>%
  mutate(no_watch_sleep = ifelse(
    hours_tv_per_school_day == "do not watch" &
    school_night_hours_sleep >= "8", "Yes", "No"))
no_tv <- not_watch_tv %>%
  filter(!is.na(no_watch_sleep))
no_watch_ci <- no_tv %>%
  specify(response = no_watch_sleep, success = "Yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_ci(level = 0.95)
no_watch_ci
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   0.0283   0.0342
no_tv <- no_tv %>%
  drop_na(no_watch_sleep)

n <- count(no_tv)
p <- count(no_tv %>% filter(no_watch_sleep == "Yes"))/n
SE <- sqrt((p*(1-p))/n)

ME <- (1.96*SE)*100

print(ME)
##           n
## 1 0.2971934

The margin of error is 0.29%

Another aspect of the data I wanted to see was the proportion of females who were physically active at least 3 days a week:

set.seed(1234)
female_phys_active <- yrbss %>%
  mutate(f_active = ifelse(
    gender == "female" &
    physically_active_7d >= "3", "Yes", "No"))
days_active <- female_phys_active %>%
  filter(!is.na(f_active))
active_ci <- days_active %>%
  specify(response = f_active, success = "Yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_ci(level = 0.95)
active_ci
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1    0.270    0.285
days_active<- days_active %>%
  drop_na(f_active)

n <- count(days_active)
p <- count(days_active %>% filter(f_active == "Yes"))/n
SE <- sqrt((p*(1-p))/n)

ME <- (1.96*SE)*100

print(ME)
##         n
## 1 0.75571

There is a 0.75% margin of error of females who are physically active 3 or more days in a week.

How does the proportion affect the margin of error?

Imagine you’ve set out to survey 1000 people on two questions: are you at least 6-feet tall? and are you left-handed? Since both of these sample proportions were calculated from the same sample size, they should have the same margin of error, right? Wrong! While the margin of error does change with sample size, it is also affected by the proportion.

Think back to the formula for the standard error: \(SE = \sqrt{p(1-p)/n}\). This is then used in the formula for the margin of error for a 95% confidence interval:

\[ ME = 1.96\times SE = 1.96\times\sqrt{p(1-p)/n} \,. \] Since the population proportion \(p\) is in this \(ME\) formula, it should make sense that the margin of error is in some way dependent on the population proportion. We can visualize this relationship by creating a plot of \(ME\) vs. \(p\).

Since sample size is irrelevant to this discussion, let’s just set it to some value (\(n = 1000\)) and use this value in the following calculations:

n <- 1000

The first step is to make a variable p that is a sequence from 0 to 1 with each number incremented by 0.01. You can then create a variable of the margin of error (me) associated with each of these values of p using the familiar approximate formula (\(ME = 2 \times SE\)).

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

Lastly, you can plot the two variables against each other to reveal their relationship. To do so, we need to first put these variables in a data frame that you can call in the ggplot function.

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

  1. 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?
ggplot(data = dd, aes(x = p, y = me)) + 
  geom_line() +
  labs(x = "Population Proportion", y = "Margin of Error")

It appears that the margin of error is at its highest when the population proportion p is at 0.5. The margin of error increases as p values increase from 0 to 0.5, and it decreases when p values are from 0.51 to 1.0.

Success-failure condition

We have emphasized that you must always check conditions before making inference. For inference on proportions, the sample proportion can be assumed to be nearly normal if it is based upon a random sample of independent observations and if both \(np \geq 10\) and \(n(1 - p) \geq 10\). This rule of thumb is easy enough to follow, but it makes you wonder: what’s so special about the number 10?

The short answer is: nothing. You could argue that you would be fine with 9 or that you really should be using 11. What is the “best” value for such a rule of thumb is, at least to some degree, arbitrary. However, when \(np\) and \(n(1-p)\) reaches 10 the sampling distribution is sufficiently normal to use confidence intervals and hypothesis tests that are based on that approximation.

You can investigate the interplay between \(n\) and \(p\) and the shape of the sampling distribution by using simulations. Play around with the following app to investigate how the shape, center, and spread of the distribution of \(\hat{p}\) changes as \(n\) and \(p\) changes.

  1. 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 distributions appear to show a normal distribution, with its center around 0.1, which is the population proportion of this sample, and its spread between 0.05 and 0.15.

  1. 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 increasing the population proportion and keeping the min and max axis the same, the shape keeps its normal distribution, its center is around the same as the population proportion, and the spread appearing a little more wider. As the population proportion gets closer to the min and max values, the distribution becomes skewed, to the left when it’s closer to the max value and to the right when closer to the min value. The shape begins to change when the population proportion is within approximately 0.05 of the min and max values, respectively, without going over the max value or under the min value. The width of the distribution appears to depend on the min and max values in relation to the population proportion. The closer the population proportion is to the center between the min and max values, the wider the shape. However, it doesn’t appear to be a significant increase in the spread when the population proportion is in between the min and max values.

  1. Now also change \(n\). How does \(n\) appear to affect the distribution of \(\hat{p}\)?

When increasing the sample size n, it appears that there is still a normal distribution. However, the spread becomes more narrow towards the size of the population proportion.

More Practice

For some of the exercises below, you will conduct inference comparing two proportions. In such cases, you have a response variable that is categorical, and an explanatory variable that is also categorical, and you are comparing the proportions of success of the response variable across the levels of the explanatory variable. This means that when using infer, you need to include both variables within specify.

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

Ho: Null hypothesis: There is no evidence that those who sleep 10+ hours per day are more likely to strength train every day of the week.

Ha: Alternative hypothesis: There is evidence that those who sleep 10+ hours per day are more likely to strength train every day of the week.

sleep_train <- yrbss %>%
  select(strength_training_7d, school_night_hours_sleep) %>%
  filter(!is.na(strength_training_7d)) %>%
  filter(!is.na(school_night_hours_sleep)) %>%
  mutate(s_training = strength_training_7d == "7", "Yes", "No") %>%
  mutate(sleep_10 = ifelse(school_night_hours_sleep == "10+", "Yes", "No")) %>%
  group_by(sleep_10) %>%
  summarise(count=n())%>%
  mutate(props=round(count/sum(count), 2))
sleep_train
## # A tibble: 2 × 3
##   sleep_10 count props
##   <chr>    <int> <dbl>
## 1 No       11907  0.97
## 2 Yes        312  0.03

3% of people slept 10+ hours per day and strength trained 7 days a week.

set.seed(1234)
school_sleep_train <- yrbss %>%
  mutate(sleepy = ifelse(school_night_hours_sleep == "10+", "yes", "no")) %>%
  mutate(s_train = ifelse(strength_training_7d == "7", "yes", "no")) %>%
  filter(!is.na(sleepy) & !is.na(s_train)) 
school_sleep_train
## # A tibble: 12,219 × 15
##      age gender grade hispanic race        height weight helme…¹ text_…² physi…³
##    <int> <chr>  <chr> <chr>    <chr>        <dbl>  <dbl> <chr>   <chr>     <int>
##  1    14 female 9     not      Black or A…  NA      NA   never   0             4
##  2    14 female 9     not      Black or A…  NA      NA   never   <NA>          2
##  3    15 female 9     hispanic Native Haw…   1.73   84.4 never   30            7
##  4    15 female 9     not      Black or A…   1.6    55.8 never   0             0
##  5    15 female 9     not      Black or A…   1.5    46.7 did no… did no…       2
##  6    15 female 9     not      Black or A…   1.57   67.1 did no… did no…       1
##  7    15 female 9     not      Black or A…   1.65  132.  did no… <NA>          4
##  8    14 male   9     not      Black or A…   1.88   71.2 never   <NA>          4
##  9    15 male   9     not      Black or A…   1.75   63.5 never   <NA>          5
## 10    15 male   10    not      Black or A…   1.37   97.1 did no… <NA>          0
## # … with 12,209 more rows, 5 more variables: hours_tv_per_school_day <chr>,
## #   strength_training_7d <int>, school_night_hours_sleep <chr>, sleepy <chr>,
## #   s_train <chr>, and abbreviated variable names ¹​helmet_12m,
## #   ²​text_while_driving_30d, ³​physically_active_7d
set.seed(1234)
school_sleep_train_ci <- school_sleep_train %>%
  specify(response = sleepy, success = "yes") %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "prop") %>%
  get_ci(level = 0.95)
school_sleep_train_ci
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   0.0228   0.0282

The confidence interval is 0.023 to 0.028

n <- count(school_sleep_train)
p <- count(school_sleep_train %>% filter(sleepy == "yes" & s_train == "yes"))/n
SE <- sqrt((p*(1-p))/n)

ME <- 1.96*SE

print(ME)
##             n
## 1 0.001465083

The margin of error is 0.0014.

We are 95% confident that the population proportion of people who sleep 10+ hours a day and strength train 7 days a week will be between 0.023 and 0.028, with a margin of error of 0.0014. Because the population proportion of this subset is 3%, we fail to reject the null hypothesis that there is no evidence that people who sleep 10+ hours per day strength train 7 days a week, even with a margin of error of 0.0014, since the percentage of students who sleep 10+ hours a day and strength train 7 days a week is outside of the confidence interval and margin of error.

  1. 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 probability 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.

You can detect a change at a significance level of 0.05 when you reject the null hypothesis when it is actually true. In this scenario, a smaller significance level would be more appropriate to use.

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

Using the previously made plot, the population proportion with a margin of error of 0.01 is approximately 0.05. With a confidence level of 95%, we know that the z-score will be approximately 1.96. Using this information, I plugged in the numbers into the formula for finding the appropriate sample size for an unknown population:

ME <- 0.01
Z <- 1.96
p <- 0.05

round((Z**2)*p*(1-p)/ME**2)
## [1] 1825

Based on the results, the appropriate amount of people to sample with a 95% confidence level and 1% margin of error is approximately 1,825. * * *