Lab 1: Intro to R

library(tidyverse)

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

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

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

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

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

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

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

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

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

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

library(openintro)

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

## Warning: package 'airports' was built under R version 4.3.3

## Warning: package 'cherryblossom' was built under R version 4.3.3

## Warning: package 'usdata' was built under R version 4.3.3

library(infer)

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

library(skimr)

## Warning: package 'skimr' was built under R version 4.3.3

data("yrbss")

Exercise 1

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

#Each case in this dataset represents a single student, the entire set being comprised of 13,583 cases.

Exercise 2

yrbss %>% 
  skim()

Data summary

Name	Piped data
Number of rows	13583
Number of columns	13
_______________________
Column type frequency:
character	8
numeric	5
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
gender	12	1.00	4	6	0	2	0
grade	79	0.99	1	5	0	5	0
hispanic	231	0.98	3	8	0	2	0
race	2805	0.79	5	41	0	5	0
helmet_12m	311	0.98	5	12	0	6	0
text_while_driving_30d	918	0.93	1	13	0	8	0
hours_tv_per_school_day	338	0.98	1	12	0	7	0
school_night_hours_sleep	1248	0.91	1	3	0	7	0

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
age	77	0.99	16.16	1.26	12.00	15.00	16.00	17.00	18.00	▁▂▅▅▇
height	1004	0.93	1.69	0.10	1.27	1.60	1.68	1.78	2.11	▁▅▇▃▁
weight	1004	0.93	67.91	16.90	29.94	56.25	64.41	76.20	180.99	▆▇▂▁▁
physically_active_7d	273	0.98	3.90	2.56	0.00	2.00	4.00	7.00	7.00	▆▂▅▃▇
strength_training_7d	1176	0.91	2.95	2.58	0.00	0.00	3.00	5.00	7.00	▇▂▅▂▅

#We are missing weights from 1,219 observations.

yrbss1 <- yrbss %>% 
  mutate(physical_3plus = if_else(physically_active_7d > 2, "yes", "no"))|>
    filter(!is.na(weight) & !is.na(physical_3plus))

Exercise 3

ggplot(data = yrbss1, aes (physical_3plus, weight)) + geom_violin(position = "dodge")

#Yes, there is a relationship between the variables "physical_3plus" and "weight". I expected there to be a greater difference between the plots of the factors since one's physical activity greatly impacts their weight, but there are countless other factors that do so too, so the results are not too shocking.

yrbss1 %>%
  group_by(physical_3plus) %>%
  summarise(mean_weight = mean(weight, na.rm = TRUE))

## # A tibble: 2 × 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4

Exercise 4

yrbss1 %>%
  group_by(physical_3plus) %>%
  summarize(n = n())

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

yrbss1 %>%
  group_by(weight) %>%
  summarize(n = n())

## # A tibble: 239 × 2
##    weight     n
##     <dbl> <int>
##  1   29.9     1
##  2   31.8     1
##  3   33.1     1
##  4   34.0     3
##  5   34.9     1
##  6   35.4     3
##  7   35.8     1
##  8   36.3     7
##  9   36.7     1
## 10   37.2     1
## # ℹ 229 more rows

#The conditions necessary for inference are satisfied: there are more than 10 success and failures for "physical_3plus" and more than 30 numerical entries for "weight".

Exercise 5

#Ho: There is no difference in the average weights of adolescents who exercise at least 3 times a week and those who do not. Ho : mu1 = mu2
# Ha: There is a difference in the average weights of adolescents who exercise at least 3 times a week and those who do not. Ha : mu1 =/ mu2

obs_diff <- yrbss1 %>%
  specify(weight ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

null_dist <- yrbss1 %>%
  specify(weight ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

ggplot(data = null_dist, aes(x = stat)) +
  geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Exercise 6

ggplot(data = null_dist, aes(x = stat)) +
   geom_histogram()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

geom_vline(xintercept = 1.774584, color="red")

## mapping: xintercept = ~xintercept 
## geom_vline: na.rm = FALSE
## stat_identity: na.rm = FALSE
## position_identity

#The line does not fall on the distribution.

Exercise 7

#None of the "null_dist" permutations have a difference equal to or larger than the "obs_diff" 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()` (`?infer::get_p_value()`) for more information.

## # A tibble: 1 × 1
##   p_value
##     <dbl>
## 1       0

Exercise 8

# I got a warning message to be cautious about reporting the p-value that the system procured as it is "0". I believe I got this warning because such a p-value is interpreted as there being absolutely no way that the null hypothesis can fail to be rejected, and such broad statements hardly hold true.

Exercise 9

percentile_ci <- null_dist |>
  get_confidence_interval(level = 0.95)
percentile_ci

## # A tibble: 1 × 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   -0.611    0.625

#Given the inclusion of 0 in the confidence interval, we are 95% sure that there is no statistically significant difference between the average weights of adolescents who exercise at least 3 times a week and those who do not.