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

Every two years, the Centers for Disease Control and Prevention conduct the Youth Risk Behavior Surveillance System (YRBSS) survey, where it takes data from high schoolers (9th through 12th grade), to analyze health patterns. You will work with a selected group of variables from a random sample of observations during one of the years the YRBSS was conducted.

Load the yrbss data set into your workspace.

data('yrbss', package='openintro')

There are observations on 13 different variables, some categorical and some numerical. The meaning of each variable can be found by bringing up the help file:

?yrbss
  1. What are the cases in this data set? How many cases are there in our sample?
nrow(yrbss)
## [1] 13583

Remember that you can answer this question by viewing the data in the data viewer or by using the following command:

glimpse(yrbss)
## Rows: 13,583
## Columns: 13
## $ age                      <int> 14, 14, 15, 15, 15, 15, 15, 14, 15, 15, 15, 1…
## $ gender                   <chr> "female", "female", "female", "female", "fema…
## $ grade                    <chr> "9", "9", "9", "9", "9", "9", "9", "9", "9", …
## $ hispanic                 <chr> "not", "not", "hispanic", "not", "not", "not"…
## $ race                     <chr> "Black or African American", "Black or Africa…
## $ height                   <dbl> NA, NA, 1.73, 1.60, 1.50, 1.57, 1.65, 1.88, 1…
## $ weight                   <dbl> NA, NA, 84.37, 55.79, 46.72, 67.13, 131.54, 7…
## $ helmet_12m               <chr> "never", "never", "never", "never", "did not …
## $ text_while_driving_30d   <chr> "0", NA, "30", "0", "did not drive", "did not…
## $ physically_active_7d     <int> 4, 2, 7, 0, 2, 1, 4, 4, 5, 0, 0, 0, 4, 7, 7, …
## $ hours_tv_per_school_day  <chr> "5+", "5+", "5+", "2", "3", "5+", "5+", "5+",…
## $ strength_training_7d     <int> 0, 0, 0, 0, 1, 0, 2, 0, 3, 0, 3, 0, 0, 7, 7, …
## $ school_night_hours_sleep <chr> "8", "6", "<5", "6", "9", "8", "9", "6", "<5"…

Exploratory data analysis

You will first start with analyzing the weight of the participants in kilograms: weight.

Using visualization and summary statistics, describe the distribution of weights. The summary function can be useful.

summary(yrbss$weight)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   29.94   56.25   64.41   67.91   76.20  180.99    1004
  1. How many observations are we missing weights from?
sum(is.na(yrbss$weight))
## [1] 1004

Next, consider the possible relationship between a high schooler’s weight and their physical activity. Plotting the data is a useful first step because it helps us quickly visualize trends, identify strong associations, and develop research questions.

First, let’s create a new variable physical_3plus, which will be coded as either “yes” if they are physically active for at least 3 days a week, and “no” if not.

yrbss <- yrbss %>% 
  mutate(physical_3plus = ifelse(yrbss$physically_active_7d > 2, "yes", "no"))
  1. Make a side-by-side boxplot of physical_3plus and weight. Is there a relationship between these two variables? What did you expect and why?
ggplot(yrbss, aes(x = physical_3plus, y = weight)) +
  geom_boxplot() +
  labs(
    title = "Boxplot of Weight by Physical Activity Level",
    x = "Physically Active at Least 3 Days a Week",
    y = "Weight (kg)"
  ) +
  theme_minimal()

It would be reasonable to expect that participants who are physically active for at least three days a week might have lower or more consistent weights compared to those who are less active. Physical activity is generally associated with maintaining a healthy weight and reducing weight variability.

The box plots show how the medians of the two distributions compare, but we can also compare the means of the distributions using the following to first group the data by the physical_3plus variable, and then calculate the mean weight in these groups using the mean function while ignoring missing values by setting the na.rm argument to TRUE.

yrbss %>%
  group_by(physical_3plus) %>%
  summarise(mean_weight = mean(weight, na.rm = TRUE))
## # A tibble: 3 × 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4
## 3 <NA>                  69.9

There is an observed difference, but is this difference statistically significant? In order to answer this question we will conduct a hypothesis test.

Inference

  1. Are all conditions necessary for inference satisfied? Comment on each. You can compute the group sizes with the summarize command above by defining a new variable with the definition n().

Independence of observation: Assuming data was collected randomly and independently by CDC, this condition is likely satisfied

Sample Size:

group_sizes <- yrbss %>%
  group_by(physical_3plus) %>%
  summarize(group_size = n())

group_sizes
## # A tibble: 3 × 2
##   physical_3plus group_size
##   <chr>               <int>
## 1 no                   4404
## 2 yes                  8906
## 3 <NA>                  273

Normality of Distribution:

ggplot(yrbss, aes(x = weight)) +
  geom_histogram(binwidth = 5) +
  facet_wrap(~ physical_3plus) +
  labs(
    title = "Distribution of Weight by Physical Activity Level",
    x = "Weight (kg)",
    y = "Count"
  ) +
  theme_minimal()

  1. Write the hypotheses for testing if the average weights are different for those who exercise at least times a week and those who don’t.

Null Hypothesis: The average weight of high school students who are physically active at least 3 times a week is equal to the average weight of those who are not.

Next, we will introduce a new function, hypothesize, that falls into the infer workflow. You will use this method for conducting hypothesis tests.

But first, we need to initialize the test, which we will save as obs_diff.

obs_diff <- yrbss %>%
  drop_na(physical_3plus) %>%
  specify(weight ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

Notice how you can use the functions specify and calculate again like you did for calculating confidence intervals. Here, though, the statistic you are searching for is the difference in means, with the order being yes - no != 0.

After you have initialized the test, you need to simulate the test on the null distribution, which we will save as null.

null_dist <- yrbss %>%
  drop_na(physical_3plus) %>%
  specify(weight ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

Here, hypothesize is used to set the null hypothesis as a test for independence. In one sample cases, the null argument can be set to “point” to test a hypothesis relative to a point estimate.

Also, note that the type argument within generate is set to permute, whichis the argument when generating a null distribution for a hypothesis test.

We can visualize this null distribution with the following code:

ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram()

  1. How many of these null permutations have a difference of at least obs_stat? Convert physical_3plus to a Factor
# Step 1: Create the 'physical_3plus' variable and convert it to a factor
yrbss <- yrbss %>%
  mutate(
    physical_3plus = ifelse(physically_active_7d > 2, "yes", "no"),
    physical_3plus = factor(physical_3plus, levels = c("yes", "no"))
  )

# Step 2: Remove rows with missing values in 'weight' or 'physical_3plus'
yrbss_filtered <- yrbss %>%
  drop_na(weight, physical_3plus)

# Step 3: Calculate the observed difference in means
obs_diff <- yrbss_filtered %>%
  specify(weight ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

# Step 4: Generate the null distribution using permutations
null_dist <- yrbss_filtered %>%
  specify(weight ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

# Step 5: Calculate how many null permutations have a difference at least as extreme as obs_diff
extreme_count <- null_dist %>%
  filter(stat >= obs_diff$stat) %>%
  nrow()

# Step 6: Visualize the null distribution
ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram(binwidth = 0.5) +
  geom_vline(xintercept = obs_diff$stat, color = "red", linetype = "dashed") +
  labs(
    title = "Null Distribution of Difference in Means",
    x = "Difference in Means",
    y = "Count"
  ) +
  theme_minimal()

# Display the observed difference and the count of extreme null statistics
obs_diff
## Response: weight (numeric)
## Explanatory: physical_3plus (factor)
## # A tibble: 1 × 1
##    stat
##   <dbl>
## 1  1.77
extreme_count
## [1] 0

extreme_count is 0, this indicates that none of the null permutations were as extreme as the observed difference, suggesting that the observed difference is statistically significant. Now that the test is initialized and the null distribution formed, you can calculate the p-value for your hypothesis test using the function get_p_value.

null_dist %>%
  get_p_value(obs_stat = obs_diff, direction = "two_sided")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

This the standard workflow for performing hypothesis tests.

  1. Construct and record a confidence interval for the difference between the weights of those who exercise at least three times a week and those who don’t, and interpret this interval in context of the data.
# Step 1: Generate the bootstrap distribution
boot_dist <- yrbss_filtered %>%
  specify(weight ~ physical_3plus) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

# Step 2: Construct the 95% confidence interval
ci <- boot_dist %>%
  get_confidence_interval(level = 0.95, type = "percentile")

# Display the confidence interval
ci
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1     1.10     2.43

The confidence interval for the difference in weights between those who exercise at least three times a week and those who do not is from 1.12 kg to 2.41 kg.

More Practice

  1. Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context.
# Step 1: Generate the bootstrap distribution for height
boot_dist_height <- yrbss_filtered %>%
  specify(response = height) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean")

# Step 2: Construct the 95% confidence interval for height
ci_height <- boot_dist_height %>%
  get_confidence_interval(level = 0.95, type = "percentile")

# Display the confidence interval
ci_height
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1     1.69     1.69

The data truly has very little variability, then this confidence interval accurately reflects the consistency of students’ heights in the dataset.

  1. Calculate a new confidence interval for the same parameter at the 90% confidence level. Comment on the width of this interval versus the one obtained in the previous exercise.
# Construct the 90% confidence interval for height
ci_height_90 <- boot_dist_height %>%
  get_confidence_interval(level = 0.90, type = "percentile")

# Display the 90% confidence interval
ci_height_90
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1     1.69     1.69
  1. Conduct a hypothesis test evaluating whether the average height is different for those who exercise at least three times a week and those who don’t.
# Step 1: Calculate the observed difference in means for height
obs_diff_height <- yrbss_filtered %>%
  specify(height ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

# Step 2: Generate the null distribution for height
null_dist_height <- yrbss_filtered %>%
  specify(height ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

# Step 3: Calculate the p-value
p_value_height <- null_dist_height %>%
  get_p_value(obs_stat = obs_diff_height, direction = "two_sided")

# Display the p-value
p_value_height
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

It is highly unlikely that the difference in average height between the two groups is due to random variation 11. Now, a non-inference task: Determine the number of different options there are in the dataset for the hours_tv_per_school_day there are.

# Count the number of different options for hours_tv_per_school_day
num_tv_options <- yrbss_filtered %>%
  summarize(num_options = n_distinct(hours_tv_per_school_day))

num_tv_options
## # A tibble: 1 × 1
##   num_options
##         <int>
## 1           8
  1. Come up with a research question evaluating the relationship between height or weight and sleep. Formulate the question in a way that it can be answered using a hypothesis test and/or a confidence interval. Report the statistical results, and also provide an explanation in plain language. Be sure to check all assumptions, state your \(\alpha\) level, and conclude in context.

Create a New Variable for Sleep Categories First, create a new categorical variable (sleep_8plus) indicating whether the student sleeps at least 8 hours per school night.

# Step 1: Convert 'school_night_hours_sleep' to numeric
yrbss_filtered <- yrbss_filtered %>%
  mutate(
    sleep_numeric = case_when(
      school_night_hours_sleep == "<5" ~ 4,
      school_night_hours_sleep == "10+" ~ 10,
      TRUE ~ as.numeric(school_night_hours_sleep)
    )
  )

# Step 2: Create 'sleep_8plus' variable
yrbss_filtered <- yrbss_filtered %>%
  mutate(sleep_8plus = ifelse(sleep_numeric >= 8, "8+ hours", "<8 hours"))

# Step 3: Remove missing values and ensure 'sleep_8plus' is a binary factor
yrbss_filtered <- yrbss_filtered %>%
  drop_na(weight, sleep_8plus) %>%
  mutate(sleep_8plus = factor(sleep_8plus, levels = c("8+ hours", "<8 hours")))

# Step 4: Calculate the observed difference in means for weight
obs_diff_weight_sleep <- yrbss_filtered %>%
  specify(weight ~ sleep_8plus) %>%
  calculate(stat = "diff in means", order = c("8+ hours", "<8 hours"))

# Step 5: Generate the null distribution for weight
null_dist_weight_sleep <- yrbss_filtered %>%
  specify(weight ~ sleep_8plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("8+ hours", "<8 hours"))

# Step 6: Calculate the p-value
p_value_weight_sleep <- null_dist_weight_sleep %>%
  get_p_value(obs_stat = obs_diff_weight_sleep, direction = "two_sided")

# Step 7: Construct the 95% confidence interval for weight difference
boot_dist_weight_sleep <- yrbss_filtered %>%
  specify(weight ~ sleep_8plus) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "diff in means", order = c("8+ hours", "<8 hours"))

ci_weight_sleep <- boot_dist_weight_sleep %>%
  get_confidence_interval(level = 0.95, type = "percentile")

# Display results
p_value_weight_sleep
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.004
ci_weight_sleep
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1    -1.67   -0.321

The confidence interval for the difference in average weight between students who sleep at least 8 hours per school night and those who sleep less is from -1.67 kg to -0.377 kg.

Interpretation: Since the confidence interval ranges from -1.67 to -0.377, this suggests that, on average, students who sleep at least 8 hours per night weigh less than those who sleep less than 8 hours. The fact that the entire interval is negative indicates that there is a significant difference in weight between the two groups, and specifically that students who sleep more tend to weigh less on average. We are 95% confident that the average weight of students who sleep at least 8 hours per night is between 1.67 kg and 0.377 kg lower than the average weight of those who sleep less. * * *