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 work space.

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?

The cases are the observations in the dataset. There are 13583 cases (observations) 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?

1004 observations were missing weights.

#Summary of the weight variable
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
# Visualize the summary of weights
ggplot(yrbss, aes(x = weight)) +
  geom_boxplot(fill = 'blue', color = 'black') +
  labs(title = 'Distribution of Weights of Participants',
       x = 'Weight (kg)',
       y = 'Density')

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 box-plot of physical_3plus and weight. Is there a relationship between these two variables? What did you expect and why?

There does not seem to be a relationship with the two variables. I expected a negative relationship between physical activity, meaning that I expected the “yes” group might have a lower median or mean weight compared to the “no” group.

# Remove missing values
yrbss <- yrbss %>% filter(!is.na(physical_3plus))

# Plotting the box-plot
ggplot(yrbss, aes(x = physical_3plus, y = weight)) +
  geom_boxplot(fill = 'blue', color = 'black') +
  labs(title = 'Weight of Participants by Physical Activity',
       x = 'Physical Activity (3+ days a week)',
       y = 'Weight (kg)')

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: 2 × 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4

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().

All conditions necessary for inference are satisfied

Conditions for Inference Satisfied? Comments
Normality Yes The box-plot suggests the distribution is not skewed
Homoscedasticity ?? Need more test. Need the p-value
Independence of Errors Yes Observations were collected independently.
Sample Size Yes Sample size is larger than 30 
yrbss %>%
  group_by(physical_3plus) %>%
  summarise(n = n())
## # A tibble: 2 × 2
##   physical_3plus     n
##   <chr>          <int>
## 1 no              4404
## 2 yes             8906
  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 (H₀): There is no difference in the average weight between those who exercise at least three times a week and those who don’t exercise at least three times a week.

Alternative Hypothesis (H₁): There is a difference in the average weight between those who exercise at least three times a week and those who don’t exercise at least three times a week.

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, which is 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?

0 null permutations have a difference of at least obs_stat.

# Calculate how many of the null permutations have a difference of at least obs_stat
null_dist |>
  filter(stat >= obs_diff$stat) |>
  count()
## Response: weight (numeric)
## Explanatory: physical_3plus (factor)
## Null Hypothesis: independence
## # A tibble: 1 × 1
##       n
##   <int>
## 1     0
visualize(null_dist, bins = 10) +
  shade_p_value(obs_stat = obs_diff, direction = "two_sided")

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.

We are 95% confident that the true difference in average weight between those who exercise at least 3 times and those who don’t is between -0.629 and .576”

null_dist %>%
  get_ci(level = 0.95)
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   -0.658    0.619

More Practice

  1. Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context.

We can say with 95% confidence that the average height in meters is between 1.689128 and 1.692818

t.test(yrbss$height, conf.level = 0.95)
## 
##  One Sample t-test
## 
## data:  yrbss$height
## t = 1796.8, df = 12363, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.689128 1.692818
## sample estimates:
## mean of x 
##  1.690973
  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.

We can say with 90% confidence that the average height in meters is between 1.689425 and 1.692521. The 90% confidence interval will be narrower than the 95% confidence interval we calculated earlier for the same data. This is because we are willing to accept a slightly higher chance of being wrong (10% vs 5%) in exchange for a more precise range around the estimated mean height.

t.test(yrbss$height, conf.level = 0.90)
## 
##  One Sample t-test
## 
## data:  yrbss$height
## t = 1796.8, df = 12363, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 90 percent confidence interval:
##  1.689425 1.692521
## sample estimates:
## mean of x 
##  1.690973
  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.

Null Hypothesis (H₀): There is no difference in the average height between those who exercise at least three times a week and those who don’t exercise at least three times a week.

Alternative Hypothesis (H₁): There is a difference in the average height between those who exercise at least three times a week and those who don’t exercise at least three times a week.

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

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

ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram(color="darkblue", fill="lightblue") +
  geom_vline(xintercept = obs_diff_height$stat, color = "red") +
  labs(title = "Null Distribution of the Difference in Means of Height by Physical Activity",
       x = "Difference in Means",
       y = "Density")

null_dist_height %>%
  get_p_value(obs_stat = obs_diff_height, direction = "two_sided")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0
  1. 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.

There are 7 different options in the dataset for the `hours_tv_per_school_day

# Find the unique values for hours_tv_per_school_day and drop NA values
yrbss <- yrbss |>
  drop_na(hours_tv_per_school_day)

length(unique(yrbss$hours_tv_per_school_day))
## [1] 7
  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.

Question: Lets consider the possible relationship between a high schoolers weight and their sleep. Is there a difference in weight with those who sleep 8 hours or more on school nights?

First we create a new variable and code it as “yes” if they slept 8 hours or more on school nights, and “no” if not.

# Create a new variable sleep_over_8 and drop NA values
yrbss <- yrbss |>
  mutate(
    sleep_over_8 = case_when(
      school_night_hours_sleep == "<5" ~ 'no',
      school_night_hours_sleep == "10+" ~ 'yes',
      school_night_hours_sleep >= 8 ~ 'yes',
      school_night_hours_sleep < 8 ~ 'no'
    )
  ) |>
  drop_na(sleep_over_8)

Next, we will make side by side box-plot of sleep_over_8 and weight to see if there is a relationship between the two variables. Also, we will also compare the means of the distributions.

# Plotting the box-plot
ggplot(yrbss, aes(x = sleep_over_8, y = weight)) +
  geom_boxplot(fill = 'blue', color = 'black') +
  labs(title = 'Weight of Participants by Sleep',
       x = 'Sleep (8+ hours on school nights)',
       y = 'Weight (kg)')

# Compare the means of the distributions
yrbss |>
  group_by(sleep_over_8) |>
  summarise(mean_weight = mean(weight, na.rm = TRUE))
## # A tibble: 2 × 2
##   sleep_over_8 mean_weight
##   <chr>              <dbl>
## 1 no                  68.2
## 2 yes                 67.2

There is an observed difference so we will conduct a hypothesis test in order to see if it is a statistically significant difference.

Null Hypothesis (H₀): There is no difference in the average weight between those who sleep 8+ hours and those who don’t sleep 8+ hours.

Alternative Hypothesis (H₁): There is a difference in the average weight between those who sleep 8+ hours and those who don’t sleep 8+ hours.

# Initialize the test
obs_diff_sleep <- yrbss |>
  specify(weight ~ sleep_over_8) |>
  calculate(stat = "diff in means", order = c("yes", "no"))

# Simulate the test on the null distribution
null_dist_sleep <- yrbss |>
  specify(weight ~ sleep_over_8) |>
  hypothesize(null = "independence") |>
  generate(reps = 1000, type = "permute") |>
  calculate(stat = "diff in means", order = c("yes", "no"))

# Visualize the null distribution
ggplot(data = null_dist_sleep, aes(x = stat)) +
  geom_histogram(color="darkblue", fill="lightblue") +
  labs(title = "Null Distribution of the Difference in Means of Weight by Sleep",
       x = "Difference in Means",
       y = "Density")

# Calculate how many of the null permutations have a difference of at least obs_stat
null_dist_sleep |>
  filter(stat >= obs_diff_sleep$stat) |>
  count()
## Response: weight (numeric)
## Explanatory: sleep_over_8 (factor)
## Null Hypothesis: independence
## # A tibble: 1 × 1
##       n
##   <int>
## 1   998
# Calculate the p-value
null_dist_sleep |>
  get_p_value(obs_stat = obs_diff_sleep, direction = "two_sided")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.004

There is no difference in the average weight between those who sleep 8+ hours and those who don’t sleep 8+ hours. We can say with 95% confidence that the true difference in average weight between those who sleep 8+ hours and those who don’t is between -0.629 and .576.