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

print('The cases are the individuals in the dataset.' )
## [1] "The cases are the individuals in the dataset."
print(paste0('The total number of cases in our sample are: ', nrow(yrbss),' cases.'))
## [1] "The total number of cases in our sample are: 13583 cases."

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?

print(paste0('We are missing weights from ', sum(is.na(yrbss$weight)),' observations.'))
## [1] "We are missing weights from 1004 observations."

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

Exercise 3

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

Answer

physical_3plus_weight <- yrbss %>% drop_na()

physical_3plus_weight %>% ggplot(aes(x = physical_3plus, y = weight)) + geom_boxplot()

Is there a relationship between these two variables? What did you expect and why?

The box plot above shows that there is a relationship between the physical_3plus and weight variables. The relations shows that the weights are similiar however those who are physically active at least three times per week weigh more than those who are not physically active at least three times per week. I was not expecting this outcome. I expected that those who are physically active at least three times would have a significantly lower weight because this is the normal outcome in the real world.

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

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

All conditions necessary for inference is satisfied. The conditions are: ##### independence - The data presented depicts school age students across the nationwide school systems and is considered independent ##### normality - The sample size (male/female school age student with age range from 14 to 18) and distribtion verifies thate the distribution is normal 

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

Exercise 5

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

Answer

H_null = There’s no difference between The average weights for those who exercise at least 3 times per week and Those who don’t.

H_alt = There is a difference between The average weights for those who exercise at least 3 times per week and those who don’t. Students who exercise at least 3 times per week have a different average weight than those who don’t.

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

Exercise 6

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

Answer

print(paste0('The total null permutations that have a difference of at least obs_stat is: ', nrow(null_dist %>% filter(stat >= obs_diff)),' cases.'))
## [1] "The total null permutations that have a difference of at least obs_stat is: 0 cases."

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()` for more information.
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

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.

yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(mean_weight = mean(weight, na.rm = TRUE), sd_weight = sd(weight, na.rm = TRUE), count = n())
## # A tibble: 3 × 4
##   physical_3plus mean_weight sd_weight count
##   <chr>                <dbl>     <dbl> <int>
## 1 no                    66.7      17.6  4404
## 2 yes                   68.4      16.5  8906
## 3 <NA>                  69.9      17.6   273
calc_weight <- yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(mean_weight = mean(weight, na.rm = TRUE), sd_weight = sd(weight, na.rm = TRUE), count = n())
t <- 1.96

#Confidence Interval for those who exercise at least three times a week
uci_exercise3plus <- calc_weight$mean_weight[2] + t*(calc_weight$sd_weight[2]/sqrt(calc_weight$count[2]))

lci_exercise3plus <- calc_weight$mean_weight[2] - t*(calc_weight$sd_weight[2]/sqrt(calc_weight$count[2]))

#Confidence Interval for who do not exercise at least three times a week
uci_noexercise <- calc_weight$mean_weight[1] + t*(calc_weight$sd_weight[1]/sqrt(calc_weight$count[1]))

lci_noexercise <- calc_weight$mean_weight[1] - t*(calc_weight$sd_weight[1]/sqrt(calc_weight$count[1]))

print(paste0('The 95% confidence interval of students who exercise at least three times a week have a range between ', round(lci_exercise3plus,4),' and Lower CI ', round(uci_exercise3plus,4)))
## [1] "The 95% confidence interval of students who exercise at least three times a week have a range between 68.1062 and Lower CI 68.7907"
print(paste0('The 95% confidence interval of students who do not exercise at least three times a week have a range between ', round(lci_noexercise,4),' and Lower CI ', round(uci_noexercise,4)))
## [1] "The 95% confidence interval of students who do not exercise at least three times a week have a range between 66.153 and Lower CI 67.1948"

More Practice

Exercise 8

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

data_tb <- as.data.frame(table(yrbss$height))
freqncy <- sum(data_tb$Freq)
m_height  <- mean(yrbss$height, na.rm = TRUE)
sd_height <- sd(yrbss$height, na.rm = TRUE)
mx_height <- max(yrbss$height, na.rm = TRUE)

n_height <- yrbss %>% 
  summarise(freqncy = table(height)) %>%
  summarise(n = sum(freqncy, na.rm = TRUE))

u_height <- m_height  + t * (sd_height/sqrt(n_height))
l_height <- m_height  - t * (sd_height/sqrt(n_height))

print(paste0('The 95% confidence interval have a range between Upper CI ', round(u_height,4),' and Lower CI ', round(l_height,4)))
## [1] "The 95% confidence interval have a range between Upper CI 1.6931 and Lower CI 1.6894"

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.

There is a very small difference difference with the Upper and Lower CI between the 90% and 95% confidence interval.

nu_height <- m_height + 1.645 * (sd_height/sqrt(n_height))
nl_height <- m_height - 1.645 * (sd_height/sqrt(n_height))


print(paste0("'The 90% confidence interval have a range between Upper CI  ", round(nu_height,4), " and Lower CI ", round(nl_height,4)))
## [1] "'The 90% confidence interval have a range between Upper CI  1.6928 and Lower CI 1.6897"

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.

exercise3plus <- yrbss %>% 
  filter(physical_3plus == "yes") %>% 
  select(height) %>% 
  na.omit()

noexercise <- yrbss %>% 
  filter(physical_3plus == "no") %>% 
  select(height) %>% 
  na.omit() 

mean_exercise3plus <- mean(exercise3plus$height)
mean_noexercise <- mean(noexercise$height)

print(paste0("'The average height for those who exercise at least three times a week is ", round(mean_exercise3plus,4), " and the average height for those who don’t is ", round(mean_noexercise,4)))
## [1] "'The average height for those who exercise at least three times a week is 1.7032 and the average height for those who don’t is 1.6656"
boxplot(exercise3plus$height, noexercise$height,
        names = c("exercise", "no_exercise"))

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.

yrbss %>% group_by(hours_tv_per_school_day) %>% summarise(n())
## # A tibble: 8 × 2
##   hours_tv_per_school_day `n()`
##   <chr>                   <int>
## 1 <1                       2168
## 2 1                        1750
## 3 2                        2705
## 4 3                        2139
## 5 4                        1048
## 6 5+                       1595
## 7 do not watch             1840
## 8 <NA>                      338

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.

Research Question

Is there a evidence that there is a relationship between high schoolers weight and the amount of sleeps they get on a school night?

Hypothesis test: Using a 95% confidence interval, conduct a one-tailed correlation test to determine whether there is a significant positive relationship between weight and sleep quality.

Null hypothesis: There is no significant relationship between weight and sleep quality.

Alternative hypothesis: There is a significant positive relationship between weight and sleep quality for high schoolers who get at least 8 hours sleep.

Assumptions: The data should be normally distributed, and the relationship between weight and sleep quality should be linear.

α level: Let’s set α at 0.05.

The Data:

hrs_sleep <- yrbss %>%
  mutate(sleep_schlnght = ifelse(yrbss$school_night_hours_sleep >= 7, "yes", "no"))

weight_sleep <- hrs_sleep %>% 
  select(weight, sleep_schlnght) %>% 
  filter(sleep_schlnght == "yes") %>%
  na.omit()

weight_nosleep <- hrs_sleep %>% 
  select(weight, sleep_schlnght) %>% 
  filter(sleep_schlnght == "no") %>%
  na.omit()

mean_sleep <- mean(weight_sleep$weight)

mean_nosleep <- mean(weight_nosleep$weight)

print(paste0("'The average weight for those who get at least 8 hours sleep is ", round(mean_sleep,4), " and the average weight for those who do not get sleep is ", round(mean_nosleep,4)))
## [1] "'The average weight for those who get at least 8 hours sleep is 67.2397 and the average weight for those who do not get sleep is 68.7439"
boxplot(weight_nosleep$weight, weight_sleep$weight,
        names = c("less_sleep", "more_sleep"))

Explanation: The study found that there is no significant positive relationship between weight and sleep quality among high schoolers. This suggests that weight may not be a significant factor that affects sleep quality.

Conclusion: Based on the statistical results, we can conclude that there is insufficient evidence to support a significant positive relationship between weight and sleep quality among high schoolers.