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:
?yrbss- 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:
Ans: There are 13583 cases in the data set. There are 13583 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- How many observations are we missing weights from?
Ans: There are 1004 missing weights.
sum(is.na(yrbss$weight))## [1] 1004Next, 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")) %>%
na.exclude()- Make a side-by-side boxplot of
physical_3plusandweight. Is there a relationship between these two variables? What did you expect and why?
No, there is relationship between physical activity and weight. They have a similar mean a long with similar outliars. Alaso the Interquartile range is very similar.
ggplot(yrbss, aes(x=weight, y=physical_3plus)) + geom_boxplot() + theme_bw()
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: 2 x 2
## physical_3plus mean_weight
## <chr> <dbl>
## 1 no 67.1
## 2 yes 68.7There is an observed difference, but is this difference statistically significant? In order to answer this question we will conduct a hypothesis test.
Inference
- 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().
Ans: Yes, we have all conditions for inference as the data is normally distributed and the samples are independent of each other. We know the data is normal based on the box plot.
- 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.
Ans: H0: Students who are physically active 3 or more days/week have the same average weight as those who are not physically active.
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.
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()
- How many of these
nullpermutations have a difference of at leastobs_stat?
visualize(null_dist) +
shade_p_value(obs_stat = obs_diff, direction = "two_sided")
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.002This the standard workflow for performing hypothesis tests.
- 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, na.rm = TRUE)) ## # A tibble: 2 x 2
## physical_3plus `sd(weight, na.rm = TRUE)`
## <chr> <dbl>
## 1 no 18.0
## 2 yes 16.4yrbss %>%
group_by(physical_3plus) %>%
summarise(mean(weight, na.rm = TRUE))## # A tibble: 2 x 2
## physical_3plus `mean(weight, na.rm = TRUE)`
## <chr> <dbl>
## 1 no 67.1
## 2 yes 68.7yrbss %>%
group_by(physical_3plus) %>%
summarise(freq = table(weight)) %>%
summarise(sum(freq))## # A tibble: 2 x 2
## physical_3plus `sum(freq)`
## <chr> <int>
## 1 no 2656
## 2 yes 5695PA_mean=68.67
PA_sd=16.42
PA_n=5695
PNA_mean= 67.15
PNA_sd=17.99
PNA_n=2656
z=1.96
UCI_PA = PA_mean+z*(PA_sd/sqrt(PA_n))
LCI_PA = PA_mean-z*(PA_sd/sqrt(PA_n))
UCI_PNA = PNA_mean+z*(PNA_sd/sqrt(PNA_n))
LCI_PNA = PNA_mean-z*(PNA_sd/sqrt(PNA_n))
c("Physically Active CI: ", LCI_PA, UCI_PA )## [1] "Physically Active CI: " "68.2435360158172" "69.0964639841828"c("Physically Not Active CI: ", LCI_PNA, UCI_PNA )## [1] "Physically Not Active CI: " "66.4658155500846"
## [3] "67.8341844499154"More Practice
- Calculate a 95% confidence interval for the average height in meters (
height) and interpret it in context.
Ans. 95% CI is between 1.694m and 1.699m, we can say with 95% certainty the average height of the population is with in the CI.
mean_height <- mean(yrbss$height, na.rm = TRUE)
sd_height <- sd(yrbss$height, na.rm = TRUE)
n_height <- yrbss %>%
summarise(freq = table(height)) %>%
summarise(n = sum(freq, na.rm = TRUE))
H_UCI <- mean_height + z*(sd_height/sqrt(n_height))
H_LCI <- mean_height - z*(sd_height/sqrt(n_height))
c(H_LCI,H_UCI)## $n
## [1] 1.694811
##
## $n
## [1] 1.699298- 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.
The CI is 1.695 and 1.699. The interval with is smaller as the variability decreases thus the lower confidence interval width.
z_2= 1.645
H_UCI2 <- mean_height + z_2*(sd_height/sqrt(n_height))
H_LCI2 <- mean_height - z_2*(sd_height/sqrt(n_height))
c(H_LCI2,H_UCI2)## $n
## [1] 1.695171
##
## $n
## [1] 1.698937- 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.
HO: There is no difference in average height of those who are physically active at least 3 days per week. HA: There is a difference.
With 95% confident that the average height of students who are physically active at least 3 days per week is between 1.706 and 1.711 and the average height of students who are not physically active at least 3 days per week is between 1.667 and 1.673.
height_exercise <- yrbss %>%
filter(physical_3plus == "yes") %>%
select(height) %>%
na.omit()
height_noexercise <- yrbss %>%
filter(physical_3plus == "no") %>%
select(height) %>%
na.omit()
n__height <- height_exercise %>%
summarise(freq = table(height)) %>%
summarise(n = sum(freq, na.rm = TRUE))
n_no_height <- height_noexercise %>%
summarise(freq = table(height)) %>%
summarise(n = sum(freq, na.rm = TRUE))
mean_nh1 <- mean(height_noexercise$height)
sd_nh1<- sd(height_noexercise$height)
mean_h1 <- mean(height_exercise$height)
sd_h1 <- sd(height_exercise$height)
mean_h1## [1] 1.709147sd_h1## [1] 0.1035233n__height## # A tibble: 1 x 1
## n
## <int>
## 1 5695mean_nh1## [1] 1.671126sd_nh1## [1] 0.1022101n_no_height## # A tibble: 1 x 1
## n
## <int>
## 1 2656z_3 = 1.96
H_LCI3 <- mean_h1 - z_3*(sd_h1/sqrt(n__height))
H_UCI3 <- mean_h1 + z_3*(sd_h1/sqrt(n__height))
NH_LCI3 <- mean_nh1 - z_3*(sd_nh1/sqrt(n_no_height))
NH_UCI3 <- mean_nh1 + z_3*(sd_nh1/sqrt(n_no_height))
c(H_LCI3,H_UCI3)## $n
## [1] 1.706458
##
## $n
## [1] 1.711835c(NH_LCI3,NH_UCI3)## $n
## [1] 1.667239
##
## $n
## [1] 1.675013- 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 7 options available.
yrbss %>%group_by(hours_tv_per_school_day)%>% summarise(n())## # A tibble: 7 x 2
## hours_tv_per_school_day `n()`
## <chr> <int>
## 1 <1 1407
## 2 1 1172
## 3 2 1738
## 4 3 1309
## 5 4 627
## 6 5+ 966
## 7 do not watch 1132- 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.
Are students havier or lighter based on sleep.
HO: There is a relationship between weight and sleep
HA: There is no relationship between weight and sleep
95% confident level
Conditions:
-Independent sample-yes - Normality - yes
Because our P-value is equal to our alpha, 0.05, we cannot reject the null hypothesis. There is no relationship between the variables
yrbss <- yrbss %>%
mutate(sleep_less = ifelse(yrbss$school_night_hours_sleep < 6, "yes", "no"))
weight_less <- yrbss %>%
select(weight, sleep_less) %>%
filter(sleep_less == "yes") %>%
na.omit()
weight_more <- yrbss %>%
select(weight, sleep_less) %>%
filter(sleep_less == "no") %>%
na.omit()
n__wight <- weight_less %>%
summarise(freq = table(weight)) %>%
summarise(n = sum(freq, na.rm = TRUE))
n_more_weight <- weight_more %>%
summarise(freq = table(weight)) %>%
summarise(n = sum(freq, na.rm = TRUE))
mean_nw1 <- mean(weight_less$weight)
sd_nw1<- sd(weight_less$weight)
mean_w1 <- mean(weight_more$weight)
sd_w1 <- sd(weight_more$weight)
mean_w1## [1] 67.81277sd_w1## [1] 16.45589n__wight## # A tibble: 1 x 1
## n
## <int>
## 1 1794mean_nw1## [1] 69.57611sd_nw1## [1] 18.56974n_more_weight## # A tibble: 1 x 1
## n
## <int>
## 1 6557z_3 = 1.96
W_LCI3 <- mean_w1 - z_3*(sd_w1/sqrt(n__wight))
W_UCI3 <- mean_w1 + z_3*(sd_w1/sqrt(n__wight))
NW_LCI3 <- mean_nw1 - z_3*(sd_nw1/sqrt(n_more_weight))
NW_UCI3 <- mean_nw1 + z_3*(sd_nw1/sqrt(n_more_weight))
c(W_LCI3,W_UCI3)## $n
## [1] 67.05128
##
## $n
## [1] 68.57426c(NW_LCI3,NW_UCI3)## $n
## [1] 69.12664
##
## $n
## [1] 70.02559df <- 6557-1
t <- qt(.05/2, df, lower.tail = FALSE)
p_value <- 2*pt(t,df, lower.tail = FALSE)
p_value## [1] 0.05