Lab 7

library(tidyverse)

## Warning: package 'tidyverse' was built under R version 4.2.3

## Warning: package 'ggplot2' was built under R version 4.2.3

## Warning: package 'tibble' was built under R version 4.2.3

## Warning: package 'tidyr' was built under R version 4.2.3

## Warning: package 'readr' was built under R version 4.2.3

## Warning: package 'purrr' was built under R version 4.2.3

## Warning: package 'dplyr' was built under R version 4.2.3

## Warning: package 'stringr' was built under R version 4.2.3

## Warning: package 'forcats' was built under R version 4.2.3

## Warning: package 'lubridate' was built under R version 4.2.3

library(openintro)

## Warning: package 'openintro' was built under R version 4.2.3

library(infer)

## Warning: package 'infer' was built under R version 4.2.3

The data

data('yrbss', package='openintro')

?yrbss

## starting httpd help server ... done

Exercise 1

What are the cases in this data set? How many cases are there in our sample? There are 13583 cases in this dataset.

ls(yrbss)

##  [1] "age"                      "gender"                  
##  [3] "grade"                    "height"                  
##  [5] "helmet_12m"               "hispanic"                
##  [7] "hours_tv_per_school_day"  "physically_active_7d"    
##  [9] "race"                     "school_night_hours_sleep"
## [11] "strength_training_7d"     "text_while_driving_30d"  
## [13] "weight"

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

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

Exercise 2

How many observations are we missing weights from?

sum(is.na(yrbss))

## [1] 9476

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

Exercise 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? Yes. We observe that the weight measurements are quite alike, although there is a more pronounced clustering of weight measures among individuals who engage in exercise compared to those who do not.

yrbss_plot <- yrbss %>% 
  mutate(physical_3plus = ifelse(yrbss$physically_active_7d > 2, "yes", "no")) %>%
  na.exclude()
ggplot(yrbss_plot, 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: 3 × 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4
## 3 <NA>                  69.9

Exercise 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(). The data appears to meet these conditions as it represents a diverse sample of students across national, state, and tribal territories, and all students appear to be independent of one another.

yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(freq = table(weight)) %>%
  summarise(n = sum(freq))

## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `summarise()` has grouped output by 'physical_3plus'. You can override using
## the `.groups` argument.

## # A tibble: 3 × 2
##   physical_3plus     n
##   <chr>          <int>
## 1 no              4022
## 2 yes             8342
## 3 <NA>             215

Exercise 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): The average weights of individuals who exercise at least three times a week and those who don’t are not significantly different.

(Ha): The average weights of individuals who exercise at least three times a week are significantly different from those who don’t exercise as frequently.

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

## Warning: Removed 946 rows containing missing values.

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

## Warning: Removed 946 rows containing missing values.

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

Exercise 6

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

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

## 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       0

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

The standard deviation is 17.638 for those who do are not physically active at least 3 days per week and 16.478 for those who are.

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

Sample size

yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(freq = table(weight)) %>%
  summarise(n = sum(freq))

## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## `summarise()` has grouped output by 'physical_3plus'. You can override using
## the `.groups` argument.

## # A tibble: 3 × 2
##   physical_3plus     n
##   <chr>          <int>
## 1 no              4022
## 2 yes             8342
## 3 <NA>             215

xnot3 <- 66.67389
nnot3 <- 4022
snot3 <- 17.63805
x3 <- 68.44847
n3 <- 8342
s3 <- 16.47832

z = 1.96

uci_not <- xnot3 + z*(snot3/sqrt(nnot3))
lci_not <- xnot3 - z*(snot3/sqrt(nnot3))
uci_not

## [1] 67.219

lci_not

## [1] 66.12878

u_ci <- x3 + z*(s3/sqrt(n3))
l_ci <- x3 - z*(s3/sqrt(n3))

u_ci

## [1] 68.80209

l_ci

## [1] 68.09485

We can state with 95% confidence that students who engage in exercise at least three times a week have an average weight ranging from 68.09 kg to 68.8 kg. Similarly, students who do not meet the minimum exercise frequency of three times a week exhibit an average weight falling between 66.13 kg and 67.22 kg with 95% confidence.

Exercise 8

Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context. We can state with 95% confidence that the mean height of the students in this population falls within the range of 1.689 meters to 1.693 meters.

tb <- as.data.frame(table(yrbss$height))
freq <- sum(tb$Freq)

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

## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

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

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

x <- abs(height_lower_2 - height_lower)
y <- abs(height_upper_2 - height_upper)

c(x,y)

## $n
## [1] 0.0002940511
## 
## $n
## [1] 0.0002940511

The lower limit of the confidence interval is higher, but the upper limit is lower at a 90% confidence level.

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

group1 <- yrbss$height[yrbss$physical_3plus == "yes"]
group2 <- yrbss$height[yrbss$physical_3plus == "no"]

# Perform a two-sample t-test to compare the means of the two groups
t_test_result <- t.test(group1, group2)

# Print the t-test result
t_test_result

## 
##  Welch Two Sample t-test
## 
## data:  group1 and group2
## t = 19.029, df = 7973.3, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.03374994 0.04150183
## sample estimates:
## mean of x mean of y 
##  1.703213  1.665587

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

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

Exercise 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 α level, and conclude in context.

Is there a difference in the sleep patterns between students shorter than the mean height and those taller than the mean height?

Based on the results, we cannot conclude that there is a statistically significant difference in school night sleep duration between students who are shorter than the mean height and those who are taller.

yrbss <- yrbss %>%
  mutate(sleep_less = ifelse(yrbss$school_night_hours_sleep < 6, "yes", "no"))

height_less <- yrbss %>% 
  select(height, sleep_less) %>% 
  filter(sleep_less == "no") %>%
  na.omit()

height_more <- yrbss %>% 
  select(height, sleep_less) %>% 
  filter(sleep_less == "yes") %>%
  na.omit()

boxplot(height_less$height, height_more$height,
        names = c("less sleep", "more sleep"))

# less sleep
less_sleep_mean <- mean(height_less$height)
less_sleep_mean

## [1] 1.692256

less_sleep_sd <- sd(height_less$height)
less_sleep_sd

## [1] 0.1042161

max <- max(height_less$height)
max

## [1] 2.11

# more sleep
more_sleep_mean <- mean(height_more$height)
more_sleep_mean

## [1] 1.685185

more_sleep_sd <- sd(height_more$height)
more_sleep_sd

## [1] 0.1059036

max_2 <- max(height_more$height)
max_2

## [1] 2.11

#Difference
diff_mean <- more_sleep_mean - less_sleep_mean
diff_mean

## [1] -0.0070715

diff_sd <- sqrt(((more_sleep_mean^2) / nrow(height_more)) + ((less_sleep_mean^2) / nrow(height_less)))
diff_sd

## [1] 0.03818596

sleep_df <- 2492-1
t_sleep <- qt(.05/2, sleep_df, lower.tail = FALSE)

# confidence interval
upper_sleep<- diff_mean + t_sleep * diff_sd
upper_sleep

## [1] 0.06780798

lower_sleep<- diff_mean - t_sleep * diff_sd
lower_sleep

## [1] -0.08195098

# p-value
p_value_sleep <- 2 * pt(t_sleep, sleep_df, lower.tail = FALSE)
p_value_sleep

## [1] 0.05

…