Lab 7 - Numerical Inference

library(tidyverse)
library(openintro)
library(infer)

Exercise 1

The cases (rows) in this dataset are students surveyed in the YRBSS for a particular year.

There are 13583 cases in this dataset.

glimpse(yrbss)

## Rows: 13,583
## Columns: 13
## $ age                      <int> 14, 14, 15, 15, 15, 15, 15, 14, 15, 15, 15...
## $ gender                   <chr> "female", "female", "female", "female", "f...
## $ grade                    <chr> "9", "9", "9", "9", "9", "9", "9", "9", "9...
## $ hispanic                 <chr> "not", "not", "hispanic", "not", "not", "n...
## $ race                     <chr> "Black or African American", "Black or Afr...
## $ height                   <dbl> NA, NA, 1.73, 1.60, 1.50, 1.57, 1.65, 1.88...
## $ weight                   <dbl> NA, NA, 84.37, 55.79, 46.72, 67.13, 131.54...
## $ helmet_12m               <chr> "never", "never", "never", "never", "did n...
## $ text_while_driving_30d   <chr> "0", NA, "30", "0", "did not drive", "did ...
## $ physically_active_7d     <int> 4, 2, 7, 0, 2, 1, 4, 4, 5, 0, 0, 0, 4, 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, ...
## $ school_night_hours_sleep <chr> "8", "6", "<5", "6", "9", "8", "9", "6", "...

Exercise 2

There are 1004 observations that are missing weight values. This represents about 8% of all cases.

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

sum(!is.na(yrbss$weight))

## [1] 12579

sum(is.na(yrbss$weight))/sum(!is.na(yrbss$weight))

## [1] 0.07981557

yrbss_new <- filter(yrbss,!is.na(yrbss$weight),!is.na(yrbss$physically_active_7d), !is.na(height))

Creating a new variable physical_3plus that breaks students into 2 groups by whether they are physically active 3 or more days a week, we see that out of the 12,579 cases, about 2/3 of students are active 3 or more days a week and 1/3 are not.

yrbss_new <- yrbss_new %>% 
  mutate(physical_3plus = ifelse(yrbss_new$physically_active_7d >= 3, "yes", "no"))

table(yrbss_new$physical_3plus)

## 
##   no  yes 
## 4022 8342

Exercise 3

I would have expected there to be a linear relationship between weight and boolean variable +3days physically active where those with characteristic (value=1) would be more trim (weigh less) than those not physically active 3 or more days a week (value=0) would weigh more.

On initial visual comparison of side by side plot, the IQR looks fairly similar with a slightly larger mean weight for those who do exercise 3 or more days a week.

boxplot(yrbss_new$weight ~ yrbss_new$physical_3plus, xlab="3+ Days Active", ylab="Weight")

outlier_check <- yrbss_new %>%

  group_by(physical_3plus) %>% 
  summarise(n = n(),
            min = fivenum(weight)[1],
            Q1 = fivenum(weight)[2],
            median = fivenum(weight)[3],
            mean = round(mean(weight),2),
            Q3 = fivenum(weight)[4],
            max = fivenum(weight)[5],
            IQR = IQR(weight)
            )

outlier_check

## # A tibble: 2 x 9
##   physical_3plus     n   min    Q1 median  mean    Q3   max   IQR
##   <chr>          <int> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 no              4022  29.9  54.4   62.6  66.7  74.8  181.  20.4
## 2 yes             8342  33.1  56.7   65.8  68.4  77.1  160.  20.4

Exercise 4

Technically no, the 2 conditions for numerical inference are not met.

We can assume independence because the data we are observing is a random sample.

However, we may not be able to assume normality of the data however because although the number of cases (sample size) we are inferring from is greater than 30, there are 14 and 15 extreme outliers in weight respectively from both groups of students (physcially active 3 or more days, not physically active 3 or more days).

# check for mild outliers 
#is min/max more extreme than mild outlier cut-off?
MOF1 = outlier_check$Q1 - 1.5*outlier_check$IQR
MOF2 = outlier_check$Q3 + 1.5*outlier_check$IQR
MOF1

## [1] 23.815 26.085

MOF2

## [1] 105.455 107.725

outlier_check$min < MOF1

## [1] FALSE FALSE

outlier_check$max > MOF2

## [1] TRUE TRUE

# check for extreme outliers
#is max more extreme than extreme outlier cut-off?

EOF2 = outlier_check$Q3 + 3*outlier_check$IQR
EOF2

## [1] 136.07 138.34

outlier_check$max > EOF2

## [1] TRUE TRUE

#calculate how many observed weights are considered extreme outliers in the dataset
yrbss_new %>% group_by(physical_3plus) %>% 
  summarise(extreme_outliers = sum(weight > EOF2))

## # A tibble: 2 x 2
##   physical_3plus extreme_outliers
##   <chr>                     <int>
## 1 no                           15
## 2 yes                          14

Exercise 5

The sampling distribution of bootstrap approximates a normally distributed bell curve per CLT.

First, need to set seed and initialize

set.seed(74996)

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

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

Viewing the histogram, we can see that the sample mean differences congregate around 0 as center of the distribution. This makes sense assuming null hypothesis that there is no difference between the 2 groups of students and therefore no statistically significant difference in the sample means of the 2 groups.

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

Exercise 6

The null hypothesis being tested is that there is no statistical relationship between days of physical activity and weight in the population high school students. Under the null hypothesis, the mean weight of group A - students who are physically active 3 or more days a week is not statistically different than the mean weight of group B students who are NOT physically active 3 or more days a week.

To test the null hypothesis, we generate many sets of simulated data assuming null hypothesis is true. In each simulation, a sample mean is calculated from group A (physical_3plus = yes) and group B (physcial_3plus = “no”). The difference between the sample means is calculated and recorded as the test statistic.

If the groups are in fact the same, then the simulated sample mean difference generated will rarely be as or more extreme than the empirical mean difference of the 2 groups.

Out of the 1000 test statistics that are generated, we want to see how many sample mean differences were observed as being more extreme in value than the actual/empirical group mean difference.

This proportion is the p-value or the probability of obtaining a test statistic at least as extreme as the observed mean difference, assuming that the 2 groups of students have the same weight.

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

The p-value in this case came back as 0 and can be interpreted as 3 or fewer actual sample differences out of the total number generated (1000 in this case) were actually less extreme than the observed group mean difference.

With a probability so low, can safely reject the null hypothesis that there is no weight difference in students who are physically active and students who are not.

Exercise 7

The Confidence Interval generated here states that in 95% of samples taken, the actual population mean difference between the 2 groups of students will fall within -.64 kilogram and .65 kilogram.

null_dist %>%
  get_ci(level=0.95)

## # A tibble: 1 x 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1   -0.650    0.643

Exercise 8

Generating a 95% confidence interval for the average height of students in the YRBSS survey, we are confident that 95 times out of 100, the mean height of the high school student population will be within 1.6891 meters and 1.6928 meters

avg_height_dist <- yrbss_new %>%
 specify(response = height) %>%
 generate(reps = 1000, type = "bootstrap") %>%
 calculate(stat = "mean")

ci_95 <- avg_height_dist %>%
 get_ci(level = 0.95)
ci_95

## # A tibble: 1 x 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1     1.69     1.69

avg_height_dist %>%
  visualize(ci_fill="turquoise") + shade_confidence_interval(ci_95, color="turquoise")

This appears to be in line with the frequency histogram of height out of the entire students surveyed

fivenum(yrbss_new$height)

## [1] 1.27 1.60 1.68 1.78 2.11

ggplot(yrbss_new, aes(x=height)) + geom_histogram()

Exercise 9

At a 90% confidence interval, the mean height of the high school student population will be within 1.6894 meters and 1.6925 meters. This is a slightly narrower range as there is less accountability to capture the population mean at 90% than there is at 95%.

ci_90 <- avg_height_dist %>%
 get_ci(level = 0.90)
ci_90

## # A tibble: 1 x 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1     1.69     1.69

avg_height_dist %>%
  visualize() + shade_confidence_interval(ci_90, color="aquamarine")

Exercise 10

boxplot(yrbss_new$height ~ yrbss_new$physical_3plus, xlab="3+ Days Active", ylab="Height")

yrbss_new %>%
  group_by(physical_3plus) %>% 
  summarise(n = n(),
            min = fivenum(height)[1],
            Q1 = fivenum(height)[2],
            median = fivenum(height)[3],
            mean = round(mean(height),2),
            Q3 = fivenum(height)[4],
            max = fivenum(height)[5],
            IQR = IQR(height)
            )

## # A tibble: 2 x 9
##   physical_3plus     n   min    Q1 median  mean    Q3   max   IQR
##   <chr>          <int> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 no              4022  1.27  1.6    1.65  1.67  1.73  2.11 0.130
## 2 yes             8342  1.27  1.63   1.7   1.7   1.78  2.11 0.15

Actual observed sample mean difference is .03 meters in height.

obs_diff_height <- yrbss_new %>%
  specify(height ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

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

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

The probability of observing a value more extreme than observed empirical mean difference is extremely low. A p-value of 0 as calculated here indicates that only 3 or less simulations out of 1000 produced a mean difference greater than actual observed difference. We reject the null hypothesis that there is no difference in mean height between students who are physically active and students who are not.

null_dist_height %>%
  get_p_value(obs_stat = obs_diff_height, 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 x 1
##   p_value
##     <dbl>
## 1       0

…