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?

There are 13 case types for a total of 13,583 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
  1. How many observations are we missing weights from?

1004 weight observations are 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"))
  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?
yrbss2 <- yrbss %>% 
  mutate(physical_3plus = ifelse(yrbss$physically_active_7d > 2, "yes", "no")) %>%
  na.exclude()
ggplot(yrbss2, aes(x=weight, y=physical_3plus)) + geom_boxplot() + theme_bw()

The results are pretty similar. Based on this data, I wouldn’t say with certainty that there is some sort of relationship between weight and physical activity for 3 days per week. However, it does look like those who partake in physical activity weight more, as based on the box plot of the “yes” answer. Q1 and Q3 are greater than they are for the “no” boxplot, and so is the median.

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

Yes, seeing as this is a randomized sample, that includes a large enough number of individuals (at least 30).

yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(freq = table(weight)) %>%
  summarise(n = sum(freq))
## # A tibble: 3 × 2
##   physical_3plus     n
##   <chr>          <int>
## 1 no              4022
## 2 yes             8342
## 3 <NA>             215
  1. 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.

H0: There is no difference in average weights for those who exercise at least 3 times a week and those who don’t. Ha: There is a difference in average weights for those who exercise at least 3 times a week and 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(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")

It looks like 0 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.
yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(sd_weight = sd(weight, na.rm = TRUE))
## # A tibble: 3 × 2
##   physical_3plus sd_weight
##   <chr>              <dbl>
## 1 no                  17.6
## 2 yes                 16.5
## 3 <NA>                17.6
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
yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(freq = table(weight)) %>%
  summarise(n = sum(freq))
## # A tibble: 3 × 2
##   physical_3plus     n
##   <chr>          <int>
## 1 no              4022
## 2 yes             8342
## 3 <NA>             215

It looks like the mean and the standard deviation are similar for both groups. Those who exercise at least three times a week have a mean of about 68.45 and a standard deviation of about 16.48. Those who do not exercise at least three times a week have a mean of about 66.67 and a standard deviation of about 17.64. These measurements are taken from 8342 individuals who do not exercise three times a week, and 4022 who do not.

mean_no_phys <- round(66.67389,2)
n_no_phys <- 4022
sd_no_phys <- 17.63805

mean_phys<- round(68.44847,2)
n_phys <- 8342
sd_phys <- 16.47832

z = 1.96
upper_ci_no_phys <- mean_no_phys + z*(sd_no_phys/sqrt(n_no_phys))
lower_ci_no_phys <- mean_no_phys - z*(sd_no_phys/sqrt(n_no_phys))

upper_ci_phys <- mean_phys + z*(sd_phys/sqrt(n_phys))
lower_ci_phys <- mean_phys - z*(sd_phys/sqrt(n_phys))

c("CI for not active:", lower_ci_no_phys, upper_ci_no_phys)
## [1] "CI for not active:" "66.1248881694363"   "67.2151118305637"
c("CI for active:", lower_ci_phys, upper_ci_phys)
## [1] "CI for active:"   "68.0963823684916" "68.8036176315084"

The confidence interval is also similar for both groups. The confidence interval for the group that does not partake in physical activity is about 66.12 - 67.22, and the confidence interval for the group that does partake in physical activity is about 68.10 - 68.80. This means that the actual range for the weight of the population of the two groups is petty similar. * * *

More Practice

  1. Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context.
mean_height <- mean(yrbss$height, na.rm = TRUE)
sd_height <- sd(yrbss$height, na.rm = TRUE)
sample_height <- yrbss %>% 
  summarise(freq = table(height)) %>%
  summarise(n = sum(freq, na.rm = TRUE))

height_upper <- mean_height + z*(sd_height/sqrt(sample_height))
height_lower <- mean_height - z*(sd_height/sqrt(sample_height))
c(height_lower,height_upper)
## $n
## [1] 1.689411
## 
## $n
## [1] 1.693071

The 95% confidence interval is between 1.689411 and 1.693071. This means that the actual height for the population of the sample is between 1.689411 and 1.693071 meters.

  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.
z2 <- 1.645
height_upper_2 <- mean_height + z2*(sd_height/sqrt(sample_height))
height_lower_2 <- mean_height - z2*(sd_height/sqrt(sample_height))
c(height_lower_2 ,height_upper_2)
## $n
## [1] 1.689705
## 
## $n
## [1] 1.692777

The 90% confidence interval is between 1.689705 and 1.692777, making the interval width 0.003072. The width for the 95% confidence interval is 0.00366. The widths are similar.

  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 height is not different for those who exercise at least three times a week and those who don’t. Ha: The average height is different for those who exercise at least three times a week and those who don’t.

height_exercise <- yrbss %>% 
  filter(physical_3plus == "yes") %>% 
  select(height) %>% 
  na.omit()
height_noexercise <- yrbss %>% 
  filter(physical_3plus == "no") %>% 
  select(height) %>% 
  na.omit()
boxplot(height_exercise$height, height_noexercise$height,
        names = c("exercise", "no_exercise"))

#no exercise
mean_height_no <- mean(height_noexercise$height)
print(mean_height_no)
## [1] 1.665587
sd_height_no<- sd(height_noexercise$height)
print(sd_height_no)
## [1] 0.1028581
max_height_no <- max(height_noexercise$height)
print(max_height_no)
## [1] 2.11
#exercise
mean_height_yes <- mean(height_exercise$height)
print(mean_height_yes)
## [1] 1.703213
sd_height_yes <- sd(height_exercise$height)
print(sd_height_yes)
## [1] 0.1032956
max_height_yes <- max(height_exercise$height)
print(max_height_yes)
## [1] 2.11
meandiff <- mean_height_yes - mean_height_no
stderror <-sqrt(((mean_height_yes^2) / nrow(height_exercise)) + ((mean_height_no^2) / nrow(height_noexercise)))

degfreedom_height <- 4022-1
tvalue_height <- qt(.05/2, degfreedom_height, lower.tail = FALSE)
rightintervalht <- meandiff + tvalue_height * stderror
leftintervalht <- meandiff - tvalue_height * stderror
pvalue_height <- 2*pt(tvalue_height,degfreedom_height, lower.tail = FALSE) #p-value

print(round(leftintervalht,2))
## [1] -0.03
print(round(rightintervalht,2))
## [1] 0.1
print(pvalue_height)
## [1] 0.05

We can see that based on the box plots, the median, Q1, and Q3 values for those who exercise are all higher than those who do not. The mean height of those who do not exercise about 1.67 meters (with a standard deviation of about .10) as compared to 1.70 meters (with a standard deviation of about .10) for those who do exercise at least 3 times a week. For both groups, the maximal height is the same: 2.11 meters.The calculated p-value is .05 meaning that we can reject the null hypothesis and accept the alternate hypothesis.

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

There are 7 options: 1, 2, 3, 4, 5+, <1, do not watch.

  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.

Research Question: Is there an association between students who sleep at least 6 hours and their mean height?

H0: There is no association between sleeping at least 6 hours and average height. Ha: There is an association between sleeping at least 6 hours and average height.

yrbss <- yrbss %>%
  mutate(at_least_6_hrs = ifelse(yrbss$school_night_hours_sleep > 6, "yes", "no"))
less_than_6_hrs_height <- yrbss %>% 
  select(height, at_least_6_hrs) %>% 
  filter(at_least_6_hrs == "no") %>%
  na.omit()
at_least_6_hrs_height <- yrbss %>% 
  select(height, at_least_6_hrs) %>% 
  filter(at_least_6_hrs == "yes") %>%
  na.omit()

boxplot(less_than_6_hrs_height$height, at_least_6_hrs_height$height,
        names = c("less_sleep", "more_sleep"))

#less than 6 hours
mean_height_5 <- mean(less_than_6_hrs_height$height)
print(round(mean_height_5,2))
## [1] 1.69
sd_height_5 <- sd(less_than_6_hrs_height$height)
print(round(sd_height_5,2))
## [1] 0.1
max_height_5 <- max(less_than_6_hrs_height$height)
print(round(max_height_5,2))
## [1] 2.11
#at least 6 hours
mean_height_6<- mean(at_least_6_hrs_height$height)
print(round(mean_height_6,2))
## [1] 1.69
sd_height_6 <- sd(at_least_6_hrs_height$height)
print(round(sd_height_6,2))
## [1] 0.1
max_height_6 <- max(at_least_6_hrs_height$height)
print(round(max_height_6,2))
## [1] 2.11
meandiff_sleep <- mean_height_6 - mean_height_5
stderror_sleep <- sqrt(((mean_height_6^2) / nrow(at_least_6_hrs_height)) +((mean_height_5^2) / nrow(less_than_6_hrs_height)))
degfreedom_sleep <- 2492-1 #t-value
tvaluesleep <- qt(.05/2, degfreedom_sleep, lower.tail = FALSE)
rightinterval_sleep <- meandiff_sleep + tvaluesleep * stderror_sleep
leftinterval_sleep <- meandiff_sleep - tvaluesleep * stderror_sleep

pvalue_sleep <- 2*pt(tvaluesleep,degfreedom_sleep, lower.tail = FALSE) #p-value

print(leftinterval_sleep)
## [1] -0.05742148
print(rightinterval_sleep)
## [1] 0.06737374
print(pvalue_sleep)
## [1] 0.05

We can see that based on the box plots, the median, Q1, and Q3 values for those who sleep an average of at least 6 hours daily are all higher than those who do not. The mean height of those who do not sleep at least 6 hours daily is about 1.69 meters (with a standard deviation of about 0.1) as compared to 1.69 meters (with a standard deviation of about 0.1) for those who do sleep at least 6 hours daily. For both groups, the maximal height is the same: 2.11 meters. Even though the plots for these groups looks different, their calculated values are similar. The calculated p-value is .05 meaning that we can reject the null hypothesis and accept the alternate hypothesis. * * *