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?

The cases in this dataset are students from the YRBSS survey. There are 13,583 cases or rows in this sample.

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?

There are 1004 NA’s which means that there are 1004 observations that we are missing weights from.

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?

There is a relationship between the two variables but certain aspects are slightly different from what I expected but it does make sense. The students who are active 3 times per week have higher weights than those who are not active. This is not what I expected but when I analyze and think about it this is correct. The reason being that the students who exercise more would have more muscle mass, which weighs more than body fat, than the students who don’t work out at least 3 times per week. The students who don’t work out though do have more outliers which means more students who are overwieght or obese when compared to those who actively work out 3 times per week.

ggplot(yrbss, aes(x=physical_3plus, y=weight)) + geom_boxplot()

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

I would say yes, all of the conditions necessary for inference have been met. The necessary conditions for inference are randomization, independence, normality and equal variance. The sample size is way above the amount required for normality which is 30 and the observations are independent.

  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.

H0: The average weights of students who exercise at least 3 times a week and those who don’t exercise at least 3 times a week ARE THE SAME

H1: The average weights of students who exercise at least 3 times a week and those who don’t exercise at least 3 times a week ARE DIFFERENT

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?
visualize(null_dist) +  
  shade_p_value(obs_stat = obs_diff, direction = "two_sided")

 sum(null_dist$stat >= obs_diff$stat | null_dist$stat <= -obs_diff$stat)
## [1] 0

The red line in this distribution represents obs_stat which equals 1.774584. None of the null permutations have a difference of at least obs_stat.

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.
set.seed(1234)
yrbss %>% group_by(physical_3plus) %>% summarise(mean_weight = mean(weight, na.rm = TRUE),
sd_weight = sd(weight, na.rm = TRUE), size_physical_3plus = n())
## # A tibble: 3 × 4
##   physical_3plus mean_weight sd_weight size_physical_3plus
##   <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
active_mean <- 68.4 # Mean weight for active students
active_sd <- 16.5 # Standard D. of weight for active students
active_n <- 8906 # Sample size of active students

not_active_mean <- 66.7 # Mean weight for not active students
not_active_sd <- 17.6 # Standard deviation of weight for not active students
not_active_n <- 4404 # Sample size of not active students

z <- qnorm(0.975)  # Z-score for a two-tailed test at 95% confidence level

# Calculate the confidence interval for those who exercise at least three times a week
upper_active <- active_mean + z * (active_sd / sqrt(active_n))
lower_active <- active_mean - z * (active_sd / sqrt(active_n))
confidence_interval_active <- c(lower_active, upper_active)

# Calculate the confidence interval for those who don't exercise at least three times a week
upper_not_active <- not_active_mean + z * (not_active_sd / sqrt(not_active_n))
lower_not_active <- not_active_mean - z * (not_active_sd / sqrt(not_active_n))
confidence_interval_not_active <- c(lower_not_active, upper_not_active)

confidence_interval_active
## [1] 68.05732 68.74268
confidence_interval_not_active
## [1] 66.1802 67.2198

Using a 95% confidence interval, the students who exercise at least 3 times per week is 68.06kg and 68.74kg. The students who are not active at least 3 times per week had an average weight of 66.18kg and 67.22kg.

More Practice

  1. Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context.
set.seed(1235)
height_mean <- mean(yrbss$height, na.rm = TRUE) # Mean for height
height_sd <- sd(yrbss$height, na.rm = TRUE) # Standard D. for height
height_n <- sum(!is.na(yrbss$height))

z_height <- qnorm(0.975)  # For a two-tailed test at 95% confidence level

# Calculate the standard error for the mean height
se_height <- height_sd / sqrt(height_n)

# Calculate the margin of error
margin_of_error_height <- z_height * se_height

# Calculate the confidence interval for height
upper_height <- height_mean + margin_of_error_height
lower_height <- height_mean - margin_of_error_height

c(lower_height, upper_height)
## [1] 1.689411 1.693071

Using a 95% confidence interval we get an average height in meters between 1.68941 and 1.69307. These results indicate that we can be 95% confident that the average height of the students fall between these numbers.

  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.
z_height_90 <- qnorm(0.95)  # For a two-tailed test at 90% confidence level

# Calculate the margin of error for the new confidence level of 90%
margin_of_error_height_90 <- z_height_90 * se_height

# Calculate the confidence interval for height at the 90% confidence level
upper_height_90 <- height_mean + margin_of_error_height_90
lower_height_90 <- height_mean - margin_of_error_height_90

c(lower_height_90, upper_height_90)
## [1] 1.689705 1.692776

Using a 90% confidence interval we get an average height in meters between 1.68971 and 1.69278. These results indicate that we can be 90% confident that the average height of the students fall between these numbers.

  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.

H0: The average heights of students who exercise at least 3 times a week and those who don’t exercise at least 3 times a week ARE THE SAME

H1: The average heights of students who exercise at least 3 times a week and those who don’t exercise at least 3 times a week ARE DIFFERENT

# Perform two-sample t-test
t.test(height ~ physical_3plus, data = yrbss)
## 
##  Welch Two Sample t-test
## 
## data:  height by physical_3plus
## t = -19.029, df = 7973.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.04150183 -0.03374994
## sample estimates:
##  mean in group no mean in group yes 
##          1.665587          1.703213

For this hypothesis test I utilized a t-test in which we concluded that the average heights of students who exercise at least 3 times a week and those who don’t exercise at least 3 times a week ARE DIFFERENT. This is the alternative hypothesis which is H1. In such test we use the p-value to determine whether to reject or fail to reject the null hypothesis. This is determined on whether the p-value is greater or less than the significance level which is 1-0.95=0.05. In this case the p-value is less than 0.05.

  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.
unique(yrbss$hours_tv_per_school_day)
## [1] "5+"           "2"            "3"            "do not watch" "<1"          
## [6] "4"            "1"            NA
length(unique(yrbss$hours_tv_per_school_day))
## [1] 8

There are 8 total different options when it comes to the hours_tv_per_school _day dataset. 7 are time intervals and there is 1 option for N/A in the dataset

  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.

Do students who sleep at least 8 hours per night develop more height than students who sleep less than 8 hours per night?

H0: The average heights of students who sleep at least 8 hours per night and those who don’t sleep at least 8 hours per night ARE THE SAME

H1: The average heights of students who sleep at least 8 hours per night and those who don’t sleep at least 8 hours per night ARE DIFFERENT