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

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

Remember that you can answer this question by viewing the data in the data viewer or by using the following command: #### Answer 1: - This data set are youth that have been surveyed by the CDC’s Youth Risk Behavior Surveillance System (YRBSS) - There are 13,583 cases in our sample.

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

Execise 2. How many observations are we missing weights from?

Answer: This dataset is missing observations 1004 cases.

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

Execise 3: 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?

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 x 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4
## 3 <NA>                  69.9

Answer 3:

There is a relationship between weight and being physically active for at least 3 days a week. The median weight for the physically active is slightly higher which might mean more muscle or growth.However, physically active are less likely to be extremely overweight.I think this is consistent with what I think.

# Simple boxplot using ggplot
ggplot(yrbss, aes(x=physical_3plus, y=weight)) + geom_boxplot()

There is an observed difference, but is this difference statistically significant? In order to answer this question we will conduct a hypothesis test.

Inference

Execise 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 4:

The necessary conditions for inference are random, normal and independent. We assume that the survey respondents were randomly selected. The sample size is ~13,000 which is more than the roughly 30 sample size required to be approximately normal. The individual observations are independent.

yrbss %>%
  group_by(physical_3plus) %>%
  summarise(size_physical_3plus = n())
## # A tibble: 3 x 2
##   physical_3plus size_physical_3plus
##   <chr>                        <int>
## 1 no                            4404
## 2 yes                           8906
## 3 <NA>                           273

Execise 5. 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: Youth 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.

H1: Youth who are physically active 3 or more days per week have a different average weight as 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 <- yrbss %>%
  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 %>%
  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()

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

Answer: 6

I don’t see the variable ‘obs_stat’ anywhere but we do have ‘obs_diff’ which is 1.77. By inspection the highest value seems to be 1.33 so none of the values are above it but I assume we take half of the obs_diff (0.8872921) and find how many elements in the ‘null_dist’ are above 0.887 or below -0.887. In that case twelve values are outside of those bounds or 1.2%.

Breaking the fourth wall here. I would love to know how to really do this problem and to understand why we led up to the problem this way.

obs_diff/2
##        stat
## 1 0.8872921
length(null_dist$stat[null_dist$stat > 0.8872921]) + length(null_dist$stat[null_dist$stat < -0.8872921])
## [1] 12

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 x 1
##   p_value
##     <dbl>
## 1       0

This the standard workflow for performing hypothesis tests.

Execise 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.

answer 7:

Using the formulas below we get a 95% confidence interval for the mean weight of all youth who exercise 3 or more days a week of (68.0573122, 68.7426878).

And we get a 95% confidence interval for the mean weight of all youth who exercise less than 3 or more days a week of (66.1801895, 67.2198105).

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 x 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
CIno_low <- 66.7 - 1.96 * 17.6 / sqrt(4404)
CIno_high <- 66.7 + 1.96 * 17.6 / sqrt(4404)

CIyes_low <- 68.4 - 1.96 * 16.5 / sqrt(8906)
CIyes_high <- 68.4 + 1.96 * 16.5 / sqrt(8906)

More Practice

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

Answer: 8

Using the formulas below we get a 95% confidence interval for the mean height of all youth of (1.691766, 1.6882342). This means that we can say with 95% confidence that the true population mean height for youth is between these numbers.

yrbss %>%
  summarise(mean_height = mean(height, na.rm = TRUE),
            sd_height = sd(height, na.rm = TRUE),
            size_height = n())
## # A tibble: 1 x 3
##   mean_height sd_height size_height
##         <dbl>     <dbl>       <int>
## 1        1.69     0.105       13583
CI95_low <- 1.69 - 1.96 * 0.105 / sqrt(13583)
CI95_high <- 1.69 + 1.96 * 0.105 / sqrt(13583)

Execise 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.

Answer: 9

Using the same formulas as in Question 8 but with a critical value of 1.645 for the different confidence level, we get a 90% confidence interval for the mean height of all youth of (1.688518, 1.691482). Since the critical value is smaller but everything else is the same the confidence interval is smaller as the confidence level goes down.

CI95_low <- 1.69 - 1.645 * 0.105 / sqrt(13583)
CI95_high <- 1.69 + 1.645 * 0.105 / sqrt(13583)

Execise 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.

Answer: 10

Using the formulas below we get a 95% confidence interval for the mean height of all youth who exercise 3 or more days a week of (1.6978608, 1.7021392).

And we get a 95% confidence interval for the mean height of all youth who exercise less than 3 or more days a week of (1.6669579, 1.6730421).

Since the intervals do not overlap, we can say with a degree of confidence that there is a statistical difference between the mean height of the two groups.

Confidence intervals and hypothesis tests are the same thing just worded differently, so I feel like I need to rephrase this into the language of hypothesis tests and not the language of confidence intervals.

yrbss %>%
  group_by(physical_3plus) %>%
  summarise(mean_height = mean(height, na.rm = TRUE),
            sd_height = sd(height, na.rm = TRUE),
            size_physical_3plus = n())
## # A tibble: 3 x 4
##   physical_3plus mean_height sd_height size_physical_3plus
##   <chr>                <dbl>     <dbl>               <int>
## 1 no                    1.67     0.103                4404
## 2 yes                   1.70     0.103                8906
## 3 <NA>                  1.71     0.107                 273
CIno_low_height <- 1.67 - 1.96 * 0.103 / sqrt(4404)
CIno_high_height <- 1.67 + 1.96 * 0.103 / sqrt(4404)

CIyes_low_height <- 1.70 - 1.96 * 0.103 / sqrt(8906)
CIyes_high_height <- 1.70 + 1.96 * 0.103 / sqrt(8906)

Execise 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.

t(table(yrbss$hours_tv_per_school_day))
##       
##          <1    1    2    3    4   5+ do not watch
##   [1,] 2168 1750 2705 2139 1048 1595         1840

Execise 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 \(\alpha\) level, and conclude in context.

# Explore the sleep variable
table(yrbss$school_night_hours_sleep)
## 
##   <5  10+    5    6    7    8    9 
##  965  316 1480 2658 3461 2692  763

Based on the boxplots below it looks like the 8 hours of sleep group has the highest mean height. Maybe sleeping more than 8 hours a day is associated with depression and that doesn’t help with growing.

# Simple boxplot using ggplot
ggplot(yrbss, aes(x=school_night_hours_sleep, y=height)) + geom_boxplot()

Let’s restate the question for confidence interval testing.

Do youth who sleep 8 or more hours a school night have a higher average height than those who sleep less than 8 hours on a school night?

yrbss <- yrbss %>% 
  mutate(sleep_8plus = ifelse(yrbss$school_night_hours_sleep == "8", "yes", ifelse(yrbss$school_night_hours_sleep == "9", "yes", ifelse(yrbss$school_night_hours_sleep == "10+", "yes", "no"))))
yrbss %>%
  group_by(sleep_8plus) %>%
  summarise(mean_height = mean(height, na.rm = TRUE),
            sd_height = sd(height, na.rm = TRUE),
            size_sleep_8plus = n())
## # A tibble: 3 x 4
##   sleep_8plus mean_height sd_height size_sleep_8plus
##   <chr>             <dbl>     <dbl>            <int>
## 1 no                 1.69     0.104             8564
## 2 yes                1.69     0.106             3771
## 3 <NA>               1.70     0.105             1248
CIno_low_heightnew <- 1.69 - 1.96 * 0.104 / sqrt(8564)
CIno_high_heightnew <- 1.69 + 1.96 * 0.104 / sqrt(8564)

CIyes_low_heightnew <- 1.69 - 1.96 * 0.106 / sqrt(3771)
CIyes_high_heightnew <- 1.69 + 1.96 * 0.106 / sqrt(3771)

Using the formulas above we get a 95% confidence interval for the mean height of all youth who sleep 8 or more hours a school night of (1.6866168, 1.6933832).

And we get a 95% confidence interval for the mean height of all youth who sleep fewer than 8 hours a school night of (1.687797, 1.6922026).

Since the intervals overlap, we cannot say with 95% confidence that there is a statistical difference between the mean height of the two groups and we fail to reject the null hypothesis. Furthermore since the means are the same we cannot reject the null hypothesis at any confidence level.