Getting Started

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.4     ✔ tidyr     1.3.1
## ✔ purrr     1.0.4     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata
library(infer)
## Warning: package 'infer' was built under R version 4.4.3

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

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

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

This data set consists of 13,580 observation on hs. students

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

How many observations are we missing weights from?

1004 weight observation is missing.

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

Make a side-by-side boxplot of physical_3plus and weight.

ggplot(yrbss, aes(x = physical_3plus, y = weight)) +
  geom_boxplot() +
  labs(title = "Weight Distribution by Physical Activity", x = "Physically Active 3+ Days", y = "Weight (kg)")
## Warning: Removed 1004 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

Is there a relationship between these two variables? What did you expect and why? According to the boxplot, students who work out three or more days a week often weigh a little more. This makes sense because those who are physically engaged may have larger muscles.

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

Are all conditions necessary for inference satisfied?

Both groups’ sample sizes are sizable enough, and the data seems to be independent.

Comment on each. You can compute the group sizes with the summarize command above by defining a new variable with the definition n().

Hypothesis test setup

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"))
## Warning: Removed 946 rows containing missing values.

null distribution

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"))
## Warning: Removed 946 rows containing missing values.

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()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

How many of these null permutations have a difference of at least obs_stat? The rarity of the observed statistic under the null hypothesis may be inferred from the histogram and p-value calculation below.

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")
## 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.
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

This the standard workflow for performing hypothesis tests. # 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.

null_dist %>%
  get_confidence_interval(level = 0.95, point_estimate = obs_diff)
## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   -0.722    0.620
yrbss %>%
  summarise(
    lower = mean(height, na.rm = TRUE) - qt(0.975, df = n()-1) * sd(height, na.rm = TRUE)/sqrt(n()),
    upper = mean(height, na.rm = TRUE) + qt(0.975, df = n()-1) * sd(height, na.rm = TRUE)/sqrt(n())
  )
## # A tibble: 1 × 2
##   lower upper
##   <dbl> <dbl>
## 1  1.69  1.69

The range that the actual difference in means is most likely to fall within is provided by the confidence interval.

90% Confidence interval for height

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

Because we are allowing for greater uncertainty, a 90% confidence interval is smaller than a 95% interval

yrbss %>%
  summarise(
    lower = mean(height, na.rm = TRUE) - qt(0.95, df = n()-1) * sd(height, na.rm = TRUE)/sqrt(n()),
    upper = mean(height, na.rm = TRUE) + qt(0.95, df = n()-1) * sd(height, na.rm = TRUE)/sqrt(n())
  )
## # A tibble: 1 × 2
##   lower upper
##   <dbl> <dbl>
## 1  1.69  1.69

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.

yrbss %>%
  summarise(
    lower = mean(height, na.rm = TRUE) - qt(0.90, df = n()-1) * sd(height, na.rm = TRUE)/sqrt(n()),
    upper = mean(height, na.rm = TRUE) + qt(0.90, df = n()-1) * sd(height, na.rm = TRUE)/sqrt(n())
  )
## # A tibble: 1 × 2
##   lower upper
##   <dbl> <dbl>
## 1  1.69  1.69

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.

The infer package may be used to do a comparable hypothesis test for height.

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.

length(unique(yrbss$hours_tv_per_school_day))
## [1] 8

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.

  labs(title = "Relationship Between Sleep and Weight", x = "Hours of Sleep on School Nights", y = "Weight (kg)")
## $x
## [1] "Hours of Sleep on School Nights"
## 
## $y
## [1] "Weight (kg)"
## 
## $title
## [1] "Relationship Between Sleep and Weight"
## 
## attr(,"class")
## [1] "labels"

We look into whether weight and sleep duration are related. Hypothesis tests comparing the mean weight of various sleep groups might be used to investigate this.