DATA606wk7Lab
PK O’Flaherty
2022-03-21
Inference for numerical data
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:
?yrbssQuestion 1
- What are the cases in this data set? How many cases are there in our sample?
The cases (or records or rows) in this data set are youth that have been surveyed by the CDC’s Youth Risk Behavior Surveillance System (YRBSS). There are 13,583 cases or records.
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 1004Question 2
- How many observations are we missing weights from?
We are missing observations from 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"))Question 3
- Make a side-by-side boxplot of
physical_3plusandweight. Is there a relationship between these two variables? What did you expect and why?
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 but there are fewer high weight outliers meaning the physically active are less likely to be extremely overweight. This is consistent with what I would expect.
# Simple boxplot using ggplot
ggplot(yrbss, aes(x=physical_3plus, y=weight)) + geom_boxplot()## Warning: Removed 1004 rows containing non-finite values (stat_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.9There is an observed difference, but is this difference statistically significant? In order to answer this question we will conduct a hypothesis test.
Inference
Question 4
- Are all conditions necessary for inference satisfied? Comment on
each. You can compute the group sizes with the
summarizecommand above by defining a new variable with the definitionn().
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: While there is not resampling, the total sample size is less than 10% of the USA teenage population.
yrbss %>%
group_by(physical_3plus) %>%
summarise(size_physical_3plus = n())## # A tibble: 3 × 2
## physical_3plus size_physical_3plus
## <chr> <int>
## 1 no 4404
## 2 yes 8906
## 3 <NA> 273Question 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.
set.seed(4591)
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()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Question 6
- How many of these
nullpermutations have a difference of at leastobs_stat?
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.8872921length(null_dist$stat[null_dist$stat > 0.8872921]) + length(null_dist$stat[null_dist$stat < -0.8872921])## [1] 12Now 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 0This the standard workflow for performing hypothesis tests.
Question 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),
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 273CIno_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)Using the formulas above 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).
Since the intervals do not overlap, we can say with a degree of confidence that there is a statistical difference between the mean weight of the two groups.
I would love to know how to really do this problem. It seems like there should be another way to do it as one confidence interval.
More Practice
Question 8
- Calculate a 95% confidence interval for the average height in meters
(
height) and interpret it in context.
yrbss %>%
summarise(mean_height = mean(height, na.rm = TRUE),
sd_height = sd(height, na.rm = TRUE),
size_height = n())## # A tibble: 1 × 3
## mean_height sd_height size_height
## <dbl> <dbl> <int>
## 1 1.69 0.105 13583CI95_low <- 1.69 - 1.96 * 0.105 / sqrt(13583)
CI95_high <- 1.69 + 1.96 * 0.105 / sqrt(13583)Using the formulas above we get a 95% confidence interval for the mean height of all youth of (1.6882342, 1.6917658). This means that we can say with 95% confidence that the true population mean height for youth is between these numbers.
Question 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.
CI90_low <- 1.69 - 1.645 * 0.105 / sqrt(13583)
CI90_high <- 1.69 + 1.645 * 0.105 / sqrt(13583)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.
Question 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.
H0: Youth who are physically active 3 or more days per week have the same average height 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 height as those who are not physically active 3 or more days per week.
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 × 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 273CIno_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)Using the formulas above 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.
Question 11
- Now, a non-inference task: Determine the number of different options
there are in the dataset for the
hours_tv_per_school_daythere are.
There are six options: none, <1, 1, 2, 3, 4, and 5+ hours of TV a day.
t(table(yrbss$hours_tv_per_school_day))##
## <1 1 2 3 4 5+ do not watch
## [1,] 2168 1750 2705 2139 1048 1595 1840Question 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 763The 2nd tallest people in the world are Montenegrins and they have a reputation in the Balkans for being the sleepiest. Is there a relationship between height and sleep?
# Simple boxplot using ggplot
ggplot(yrbss, aes(x=school_night_hours_sleep, y=height)) + geom_boxplot()## Warning: Removed 1004 rows containing non-finite values (stat_boxplot).
Based on the boxplots above 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.
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"))))I wish I could have used an OR instead of nested loops but I couldn’t figure out how.
H0: Youth who sleep 8 or more hours a school night have the same average height as those who do not sleep 8 or more hours a school night.
H1: Youth who sleep 8 or more hours a school night do not have the same average height as those who do not sleep 8 or more hours a school night.
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 × 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 1248CIno_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.6669579, 1.6730421).
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.