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.

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
## starting httpd help server ... done

Exercise 1

What are the cases in this data set? How many cases are there in our sample?

Answer:

Each case in this dataset represents an individual high school student who participated in the CDC’s Youth Risk Behavior Surveillance System (YRBSS) survey. The dataset contains 13,583 cases, representing 13,583 student records. I confirmed the sample size by reviewing the dataset documentation and checking the total number of rows in R.

cat("I used help to view the author-provided explanation: \n")
## I used help to view the author-provided explanation:
help("yrbss")

cat("\nI reviewed the structure of the data with this:")
## 
## I reviewed the structure of the data with this:
str(yrbss)
## tibble [13,583 × 13] (S3: tbl_df/tbl/data.frame)
##  $ age                     : int [1:13583] 14 14 15 15 15 15 15 14 15 15 ...
##  $ gender                  : chr [1:13583] "female" "female" "female" "female" ...
##  $ grade                   : chr [1:13583] "9" "9" "9" "9" ...
##  $ hispanic                : chr [1:13583] "not" "not" "hispanic" "not" ...
##  $ race                    : chr [1:13583] "Black or African American" "Black or African American" "Native Hawaiian or Other Pacific Islander" "Black or African American" ...
##  $ height                  : num [1:13583] NA NA 1.73 1.6 1.5 1.57 1.65 1.88 1.75 1.37 ...
##  $ weight                  : num [1:13583] NA NA 84.4 55.8 46.7 ...
##  $ helmet_12m              : chr [1:13583] "never" "never" "never" "never" ...
##  $ text_while_driving_30d  : chr [1:13583] "0" NA "30" "0" ...
##  $ physically_active_7d    : int [1:13583] 4 2 7 0 2 1 4 4 5 0 ...
##  $ hours_tv_per_school_day : chr [1:13583] "5+" "5+" "5+" "2" ...
##  $ strength_training_7d    : int [1:13583] 0 0 0 0 1 0 2 0 3 0 ...
##  $ school_night_hours_sleep: chr [1:13583] "8" "6" "<5" "6" ...
cat("\nI counted cases by counting the number of rows in the dataset: \n")
## 
## I counted cases by counting the number of rows in the dataset:
cat("Number of cases: ", nrow(yrbss))
## Number of cases:  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

Exercise 2

How many observations are we missing weights from?

Answer:

There are 1,004 observations with missing weight values in the dataset.

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 %>%
  drop_na(weight) %>%
  mutate(physical_3plus = ifelse(physically_active_7d > 2, "yes", "no")) %>%
  drop_na(physical_3plus)

Exercise 3

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?

Answer:

The boxplots suggest that there is no strong relationship between physical activity and weight, as both groups have similar distributions with only minor differences in their medians and ranges. I did not expect a significant relationship because being physically active does not necessarily lead to weight loss. For example, students who exercise regularly or lift weights may gain muscle mass, which can increase their overall body weight.

ggplot(yrbss, aes(x = physical_3plus, y = weight)) +
  geom_boxplot() +
  labs(x = "physically active 3+ days", y = "weight", title = "Weight by physical activity levels boxplots") +
  theme_minimal()

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

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

Inference

Exercise 4

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

Answer:

To help evaluate the conditions, I can calculate group sizes, average weights, and standard deviation for the two groups.

yrbss %>%
  group_by(physical_3plus) %>%
  summarise(mean_weight = mean(weight, na.rm = TRUE),
            group_size = n(),
            sd_weight = sd(weight, na.rm = TRUE))

Conditions for Statistical Inference

1. Random Sampling: The YRBSS is a CDC survey designed to collect health-related information from young people. Its sampling design should be considered when generalizing the results.

2. Independence: The two groups consist of individuals who exercise three or more times per week and those who do not. The observations can be treated as independent if each individual contributes only one observation.

3. Large Sample Sizes: The “no” group has 4,022 individuals, and the “yes” group has 8,342. Both sample sizes are large enough for the Central Limit Theorem to apply.

4. Similar Variances: The standard deviations are 17.64 for the “no” group and 16.48 for the “yes” group. These values are relatively close, suggesting similar variances.

Conclusion: The large sample sizes and similar standard deviations support comparing the mean weights of the two groups. However, the survey design and independence assumptions should also be considered.

Exercise 5

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.

Answer:

The purpose of this hypothesis test is to determine whether the population mean weight differs between students who are physically active at least three days per week and students who are active fewer than three days per week.

The hypotheses are:

Null hypothesis (H₀): μyes − μno = 0. There is no difference in the population mean weights between the two groups.

Alternative hypothesis (H₁): μyes − μno ≠ 0. There is a difference in the population mean weights between the two groups.

Here, μyes represents the population mean weight of students who are physically active at least three days per week, while μno represents the population mean weight of students who are active fewer than three days per week.

Calculate the observed difference in means

First, remove observations with missing values for physical activity or weight, then calculate the observed difference in mean weights.

act_e3d <- yrbss %>%
  drop_na(physical_3plus, weight) %>%
  specify(weight ~ physical_3plus)

act_obs <- act_e3d %>%
  calculate(
    stat = "diff in means",
    order = c("yes", "no")
  )

act_obs

Generate the null distribution

Next, simulate 20,000 permutations under the null hypothesis of independence. This allows us to evaluate whether the observed difference in mean weights is unusual under the null hypothesis.

set.seed(123)

act_null_dist <- act_e3d %>%
  hypothesize(null = "independence") %>%
  generate(reps = 20000, type = "permute") %>%
  calculate(
    stat = "diff in means",
    order = c("yes", "no")
  )

ggplot(act_null_dist, aes(x = stat)) +
  geom_histogram(bins = 30) +
  labs(
    title = "Null Distribution of the Difference in Mean Weights",
    x = "Difference in mean weights (yes - no), in kilograms",
    y = "Count"
  ) +
  theme_minimal()

Calculate the p-value

Use the simulated null distribution to calculate the two-sided p-value.

act_pv <- act_null_dist %>%
  get_p_value(
    obs_stat = act_obs$stat,
    direction = "two-sided"
  )
## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.
act_pv

If the simulated p-value is zero, this does not mean the true p-value is exactly zero. It means that none of the 20,000 simulated permutations produced a difference as extreme as the observed difference. With 20,000 simulations, a zero count corresponds to an estimated p-value below approximately 0.00015 using the common plus-one correction.

Calculate a 95% confidence interval

A bootstrap confidence interval estimates the population difference in mean weights between the two groups.

set.seed(123)

act_ci <- act_e3d %>%
  specify(weight ~ physical_3plus) %>%
  generate(reps = 5000, type = "bootstrap") %>%
  calculate(
    stat = "diff in means",
    order = c("yes", "no")
  ) %>%
  get_ci(level = 0.95, type = "percentile")

act_ci

Calculate the margin of error:

round(
  (act_ci$upper_ci - act_ci$lower_ci) / 2,
  6
)
## [1] 0.655452

Interpretation and conclusion

The observed difference in mean weights is approximately 1.77 kg, with students who exercise at least three days per week having a higher average weight. The 95% confidence interval ranges from 1.12 to 2.43 kg. Since the interval does not include zero and the simulated p-value is below 0.05, we reject the null hypothesis. These results suggest a statistically significant association between physical activity and average weight.

Exercise 6

How many of these null permutations have a difference of at least obs_stat?

Answer: 

obs_count <- act_null_dist %>%
  filter(abs(stat) >= abs(act_obs$stat)) %>%
  summarise(count = n())

cat("\nThe number of null permutations with a difference at least as extreme as the observed statistic is:\n")
## 
## The number of null permutations with a difference at least as extreme as the observed statistic is:
obs_count

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.

act_pv_rev <- act_null_dist %>%
  get_p_value(
    obs_stat = act_obs$stat,
    direction = "two-sided"
  )
## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.
cat("The P-value is:\n")
## The P-value is:
act_pv_rev

The number of null permutations with a difference at least as extreme as the observed statistic is zero. The calculated p-value is also reported as zero, indicating strong evidence against the null hypothesis. However, a p-value of zero does not mean the true p-value is exactly zero. It means that none of the simulated permutations produced a result as extreme as the observed statistic. More simulations may provide a more precise estimate of the p-value.

Exercise 7

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.

Answer:

The 95% confidence interval for the difference in mean weights (exercise group minus non-exercise group) is 1.10 to 2.45. This suggests that individuals who exercise at least three times a week weigh, on average, between 1.10 and 2.45 units more than those who do not. Since the interval does not include zero, the difference is statistically significant. The margin of error is approximately 0.675. These results indicate a difference in average weight between the two groups but do not establish that physical activity causes the difference.

act_ci_rev <- act_e3d %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "diff in means", order = c("yes", "no")) %>%
  get_ci(level = 0.95, type = "percentile") 

cat("\nThe confidence interval is:\n")
## 
## The confidence interval is:
act_ci_rev
cat("\nThe margin of error is:\n")
## 
## The margin of error is:
round((act_ci_rev$upper_ci - act_ci_rev$lower_ci)/2, 6)
## [1] 0.674622

More Practice

Exercise 8

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

Answer:

The 95% confidence interval for the average height is 1.689 to 1.693 meters. This means we are 95% confident that the population’s average height falls within this range. The sample mean is 1.691 meters, and the median is 1.680 meters. Although individual heights range from 1.27 to 2.11 meters, the confidence interval is narrow because the large sample size provides a precise estimate of the average height.

height_ci <- yrbss %>%
  drop_na(height) %>%
  specify(response = height) %>%
  generate(reps = 5000, type = "bootstrap") %>%
  calculate(stat = "mean") %>%
  get_ci(level = 0.95, type = "percentile")

cat("The lower confidence interval is:\n")
## The lower confidence interval is:
round(height_ci$lower_ci, 6)
## [1] 1.689157
cat("\nThe upper confidence interval is\n")
## 
## The upper confidence interval is
round(height_ci$upper_ci, 6)
## [1] 1.692784
cat("The summary stats for height data in YRBSS are:\n")
## The summary stats for height data in YRBSS are:
summary(yrbss$height)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.270   1.600   1.680   1.691   1.780   2.110

Exercise 9

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.

Answer:

height_ci_90 <- yrbss %>%
  drop_na(height) %>%
  specify(response = height) %>%
  generate(reps = 5000, type = "bootstrap") %>%
  calculate(stat = "mean") %>%
  get_ci(level = 0.90, type = "percentile")

cat("The lower confidence interval is:\n")
## The lower confidence interval is:
round(height_ci_90$lower_ci, 6)
## [1] 1.689398
cat("The difference between 95% and 90% confidence lower CI is:\n")
## The difference between 95% and 90% confidence lower CI is:
 height_ci$lower_ci - height_ci_90$lower_ci 
## [1] -0.0002413256
cat("The difference between 95% and 90% confidence upper CI is:\n")
## The difference between 95% and 90% confidence upper CI is:
 height_ci$upper_ci - height_ci_90$upper_ci 
## [1] 0.0002765893
 cat("The spread changed from ", height_ci$upper_ci - height_ci$lower_ci," for 95% confidence to ",
    height_ci_90$upper_ci - height_ci_90$lower_ci, " for 90% confidence. ", "a difference of around 6 ten thousands of a percent.")
## The spread changed from  0.003626982  for 95% confidence to  0.003109067  for 90% confidence.  a difference of around 6 ten thousands of a percent.

The 90% confidence interval for the average height is approximately 1.689 to 1.692 meters. Compared with the 95% confidence interval, this interval is slightly narrower, decreasing from 0.00363 to 0.00311 meters. This is expected because a lower confidence level produces a narrower interval. Although the difference is small, both intervals estimate the average height to be approximately 1.69 meters. It would be interesting to compare this result with the average height of the broader population to determine how representative the sample is.

Exercise 10

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.

Answer:

H₀: The population mean heights are equal for students who exercise at least three times per week and those who do not.

H₁: The population mean heights are different between the two groups.

Significance level: α = 0.05.

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

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

dh3_pv <- dh3null_dist %>%
  get_p_value(obs_stat = dh3_obs, direction = "two_sided")
## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.
dh3_ci <- yrbss %>%
  drop_na(physical_3plus, height) %>%
  specify(height ~ physical_3plus) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "diff in means", order = c("yes", "no")) %>%
  get_ci(level = 0.95, type = "percentile")

cat("The P value is:\n")
## The P value is:
dh3_pv
cat("\nThe confidence interval is:\n")
## 
## The confidence interval is:
dh3_ci
round((dh3_ci$upper_ci - dh3_ci$lower_ci)/2,6)
## [1] 0.003892

The p-value is approximately zero, indicating strong evidence against the null hypothesis. Therefore, there is a statistically significant difference in average height between individuals who exercise at least three times a week and those who do not. The 95% confidence interval for the difference in mean height is approximately 0.034 to 0.042 meters. Although the difference is statistically significant, it is relatively small. These results show an association between exercise frequency and average height but do not establish a causal relationship.

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

Answer:

The hours_tv_per_school_day variable has 7 unique non-missing values in the dataset. If missing values (NA) are included, there are 8 distinct options in total.

n_distinct(yrbss$hours_tv_per_school_day, na.rm = TRUE)
## [1] 7

Exercise 12

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 α level, and conclude in context.

Answer:

Question: does the average weight of individuals differ based on how many hours they sleep on school nights?

ws_obs <- yrbss %>%
  drop_na(weight, school_night_hours_sleep) %>%
  specify(weight ~ school_night_hours_sleep) %>%
  calculate(stat = "F")

ws_null <- yrbss %>%
  drop_na(weight, school_night_hours_sleep) %>%
  specify(weight ~ school_night_hours_sleep) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "F")

ws_pv <- ws_null %>%
  get_p_value(obs_stat = ws_obs, direction = "greater")
## Warning: Please be cautious in reporting a p-value of 0. This result is an approximation
## based on the number of `reps` chosen in the `generate()` step.
## ℹ See `get_p_value()` (`?infer::get_p_value()`) for more information.
ws_ci <- yrbss %>%
  drop_na(weight, school_night_hours_sleep) %>%
  specify(weight ~ school_night_hours_sleep) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "F") %>%
  get_ci(level = 0.95, type = "percentile")

cat("The P value is:\n")
## The P value is:
ws_pv
cat("\nThe confidence interval is:\n")
## 
## The confidence interval is:
ws_ci
cat("\nThe margin of error is:\n")
## 
## The margin of error is:
round((ws_ci$upper_ci - ws_ci$lower_ci)/2,6)
## [1] 4.396171

Hypotheses

Null Hypothesis (H₀): The average weight is the same across all sleep-duration groups.

Alternative Hypothesis (H₁): At least one sleep-duration group has a different average weight.

Significance Level: α = 0.05

yrbss %>%
  drop_na(weight, school_night_hours_sleep) %>%
  group_by(school_night_hours_sleep) %>%
  summarise(
    group_size = n(),
    mean_weight = mean(weight),
    sd_weight = sd(weight),
    .groups = "drop"
  )

Assumptions

The observations should be independent, and each sleep-duration group should have sufficient observations. The YRBSS is a CDC survey, so its sampling design should be considered. Group sizes and weight distributions should also be checked to confirm that the test is appropriate.

Results and Interpretation

The simulated p-value is approximately zero, providing strong evidence against the null hypothesis. The 95% bootstrap confidence interval for the F-statistic ranges from 3.74 to 12.54. This interval describes uncertainty in the F-statistic, not the difference in average weights.

At α = 0.05, we reject the null hypothesis and conclude that average weight differs significantly across sleep-duration groups. However, this does not mean that sleep duration directly causes weight differences. Further analysis considering age and sex could provide additional insight.