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?

Answer:

According to glimpse() function there are in total 13583 cases in this data set and each case has 13 attributes which are characteristics of 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?

Answer:

missing_weight = sum(is.na(yrbss$weight))
missing_weight
## [1] 1004

The data set has in total 1004 missing values in weight column. It can be confirmed from the summary command above the question too.

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?

Answer:

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

we can see that there are some missing values and we were expecting that because in Question 2 we saw that weight coulmn has 1004 missing values. So we can exclude that from the graph by na.exclude()

yrbss2 <- na.exclude(yrbss)
ggplot(yrbss2, aes(x=weight, y= physical_3plus))+geom_boxplot()+theme_bw()

Relationship:

Well, I was expecting students who were less physically active to weight more than those who are physically active but data depicts an anti climax. According to the graph, the students who were physically active for more than three days weigh slightly more than those students who were not. But the difference seems not be that significant and requires further investigation.

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

Answer:

The conditions for inference are: * To check if our data is random enough * The variable we are measuring has normal distribution
* The population standard deviation is known.

From the overview of data it looks random and normally distributed.

Computing the group sizes:

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 times a week and those who don’t.

Answer:

Null Hypothesis (Ho) = Students who are physically active 3 or more days per week have the same average weight as those who are not physically active 3 or more days per week

Alternate Hypothesis (Ha) = Students who are physically active 3 or more days per week have a different average weight when compared to those who are not physically active 3 or more days per week.

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 <- yrbss2 %>%
  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.

set.seed(1122)
null_dist <- yrbss2 %>%
  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()

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

Answer:

We will use visualize() from infer package to plot the data in null_dist.

visualize(null_dist) +
  shade_p_value(obs_stat = obs_diff, direction = "two_sided")

The red line on the graph being our indicator of the obs_stat it does appear to be far from the data

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.

Answer:

We know that the confidence interval = sample mean +- margin of error. Similarly, margin of error is,

M_E = Z * sqrt((sd*sd)/n)

Mow in order to find out Confidence interval all we have find is the standard deviation for physically active and non active students. Since we have already found sample mean and n for both in previous questions

# standard deviation 
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

Now we have all the values and we can plug it into the formula to find confidence intervals

For active students:

a_sd <- 16.5
a_n <- 8342
a_mean <- 68.4
z <- 1.96

a_upper <- a_mean + z * sqrt((a_sd*a_sd)/a_n)
a_lower <- a_mean - z * sqrt((a_sd*a_sd)/a_n)

print(a_upper)
## [1] 68.75408
print(a_lower)
## [1] 68.04592

For non-active students:

n_sd <- 17.6
n_n <- 4022
n_mean <- 66.7
z <- 1.96

n_upper <- n_mean + z * sqrt((n_sd*n_sd)/n_n)
n_lower <- n_mean - z * sqrt((n_sd*n_sd)/n_n)

print(n_upper)
## [1] 67.24394
print(n_lower)
## [1] 66.15606

We are 95% confident that students who are active, have their between 68.045 and 68.754 while those who were non-active have their weight between 66.156 and 67.24.

More Practice

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

Answer:

First we have to the standard deviation (sd), mean and sample size (n)

h_mean <- mean(yrbss$height, na.rm = TRUE)
h_sd <- sd(yrbss$height, na.rm = TRUE)
h_n <- yrbss %>% 
  summarise(freq = table(height)) %>% 
  summarise(n = sum(freq, na.rm = TRUE))
# z value is already defined
h_upper<- h_mean + z * (h_sd / sqrt(h_n))
print(h_upper)
##          n
## 1 1.693071
h_lower <- h_mean - z * (h_sd / sqrt(h_n))
print(h_lower)
##          n
## 1 1.689411

We 95% confident that the height in meters is between 1.689 and 1.693

  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.

Answer:

h_mean <- mean(yrbss$height, na.rm = TRUE)
h_sd <- sd(yrbss$height, na.rm = TRUE)
h_n <- yrbss %>% 
  summarise(freq = table(height)) %>% 
  summarise(n = sum(freq, na.rm = TRUE))
z_h = 1.645 # z value for 90% confidence interval
h_upper<- h_mean + z_h * (h_sd / sqrt(h_n))
print(h_upper)
##          n
## 1 1.692777
h_lower <- h_mean - z_h * (h_sd / sqrt(h_n))
print(h_lower)
##          n
## 1 1.689705

The interval narrows down to 1.6897 and 1.692 as the confidence interval reduces to 90%

  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.

Answer:

Null Hypothesis (Ho) = There is no difference in the average height of those who are physically active 3 or more days per week Alternate hypothesis (Ha) = There is difference in the average height of those who are physically active 3 or more days per week

obs_diff_hgt <- yrbss2 %>% # NA Problem
  specify(height ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
set.seed(1123)
null_dist_hgt <- yrbss2 %>% # NA Problem
  specify(height ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
visualize(null_dist_hgt) + 
  shade_p_value(obs_stat = obs_diff_hgt, direction = "two_sided")

Lets find out the p-value

null_dist_hgt %>%
  get_p_value(obs_stat = obs_diff_hgt, direction = "two_sided")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

Since the p-values is below 0.05, we reject the null hypothesis. There is a difference in the average height of the students who are physically active and those who are not.

  1. Now, a non-inference task: Determine the number of different options there are in the data set for the hours_tv_per_school_day there are.

Answer:

yrbss%>%
  filter(!is.na(hours_tv_per_school_day))%>%
  select(hours_tv_per_school_day)%>%
  unique()
## # A tibble: 7 × 1
##   hours_tv_per_school_day
##   <chr>                  
## 1 5+                     
## 2 2                      
## 3 3                      
## 4 do not watch           
## 5 <1                     
## 6 4                      
## 7 1

There are 7 different options for the data set hours_tv_per_school_day and 1 option for NA.

  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.

Answer:

Question: Does students who are taller have more weight?

Null Hypothesis (Ho): There is no relationship between height and weight of the student Alternate Hypothesis (Ha): There is a relationship between height and weight of the student

yrbss <- yrbss %>%
  mutate(sleep_less = ifelse(yrbss$school_night_hours_sleep < 6, "yes", "no"))|>
  filter(!is.na(sleep_less))
obs_diff_hw <- yrbss %>%
  specify(height ~ sleep_less) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
set.seed(1123)
null_dist_hw <- yrbss %>%
  specify(height ~ sleep_less) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
null_dist_hw %>%
  get_p_value(obs_stat = obs_diff_hw, direction = "two_sided")
## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1   0.002

Since p-value is less than .05 so we can reject the null hypothesis and choose alternate hypothesis over null