Exploratory Data Analysis (EDA)

By running an EDA on weight, we can see that weight is relatively right-skewed, meaning that most youth are of lower weight relative to the range of weights, and few youth are obese in the sample.

skim(yrbss) #count number of NAs in weight - 1004; summary statistics

Variable type: character

Variable type: numeric

yrbss %>%
  ggplot(aes(x = weight)) +
    geom_density() +
    theme_bw() +
    labs(title = "Distribution of Youth Weight",
         x = "Weight", 
         y = "Density")

Relationship Between Weight and Physical Activity Level

I am someone who exercises (ie. walks) around 3 days a week, sometimes for fewer days. I am therefore interested to see whether the weight of the youth whom exercised fewer than 3 days a week significantly differed from those whom exercised more.

Firstly, I created a categorical variable: physical_3plus which takes the value yes when that individual exercises for more than 2 days a week, and no when that individual exercises for fewer than 3 days a week. All observations without an activity level are excluded from the analysis to prevent it from skewing the analysis in either direction. There is almost half of individuals as active (>2) individuals who are less active (<3).

# Create indicator for activity level
yrbss <- yrbss %>%
  mutate(physical_3plus = ifelse(physically_active_7d >= 3, "yes", "no"))

#Count number of individuals active versus not active using count
count <- yrbss %>%
  count(physical_3plus) %>%
  mutate(percentage = n/sum(n) * 100)

count

## # A tibble: 3 × 3
##   physical_3plus     n percentage
##   <chr>          <int>      <dbl>
## 1 no              4404      32.4 
## 2 yes             8906      65.6 
## 3 <NA>             273       2.01

#Count number of individuals active versus not active using group_by & summarise
group_by <- yrbss %>%
  group_by(physical_3plus) %>%
  summarise(count = n()) %>%
  mutate(percentage = count/sum(count) * 100)

group_by

## # A tibble: 3 × 3
##   physical_3plus count percentage
##   <chr>          <int>      <dbl>
## 1 no              4404      32.4 
## 2 yes             8906      65.6 
## 3 <NA>             273       2.01

I then used the physical_3plus variable to calculate the 95% confidence interval of proportion of youth who are not active. The resulting confidence interval was between 31.6% and 33.2%, meaning that we are 95% confident that the proportion of youth in the population who are not active is between 31.6% and 33.2%.

prop.test(x = 4404, # number of high schoolers who exercise fewer than 3 days a week
          n = 13583, # number of total trials
          p = 0.5, #null hypothesis being that the proportion of population of activity levels equal 
          alternative = "two.sided", # two-tailed alternative hypothesis
          conf.level = 0.95, # 95% confidence internal assuming normal distribution
          correct = F)

## 
##  1-sample proportions test without continuity correction
## 
## data:  4404 out of 13583
## X-squared = 1679, df = 1, p-value <2e-16
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.316 0.332
## sample estimates:
##     p 
## 0.324

Graphically, the medians of weight of youth don’t look too different between the two groups and the IQR either, implying that, weights do not differ depending on activity level. However, there are many outliers in both active and non-active groups, which may make the graphic analysis using medians may be misleading, and that we should use mean instead.

yrbss %>%
  filter(!is.na(physical_3plus)) %>%
  ggplot(aes(x = physical_3plus, y = weight)) +
  geom_boxplot() +
  labs(title = "Boxplot of weights for different activity levels",
       x = "Activity level",
       y = "Weight") +
  NULL

#Calculate summary statistics
ci_using_formulas <- yrbss %>%
  select(weight, physical_3plus) %>%
  filter(!is.na(physical_3plus)) %>%
  group_by(physical_3plus) %>%
  summarise(mean = mean(weight, na.rm = TRUE),
            sd = sd(weight, na.rm = TRUE),
            count = n(),
            se = sd / sqrt(count),
            t = qt(0.975, count-1),
            margin = t*se,
            l_ci = mean - margin,
            h_ci = mean + margin)

t.test(weight ~ physical_3plus,
       alternative = "two.sided",
       conf.level = 0.95,
       data = yrbss)

## 
##  Welch Two Sample t-test
## 
## data:  weight by physical_3plus
## t = -5, df = 7479, p-value = 9e-08
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -2.42 -1.12
## sample estimates:
##  mean in group no mean in group yes 
##              66.7              68.4

obs_diff <- yrbss %>%
  filter(!is.na(physical_3plus)) %>% # filter so that there are no NAs
  specify(weight ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

obs_diff #1.77 difference in means in this sample

## Response: weight (numeric)
## Explanatory: physical_3plus (factor)
## # A tibble: 1 × 1
##    stat
##   <dbl>
## 1  1.77

set.seed(1234)

null_dist <- yrbss %>%
  filter(!is.na(physical_3plus)) %>%
  # specify variables
  specify(weight ~ physical_3plus) %>%
  
  # assume independence, i.e, there is no difference
  hypothesize(null = "independence") %>%
  
  # generate 1000 reps, of type "permute"
  generate(reps = 1000, type = "permute") %>%
  
  # calculate statistic of difference, namely "diff in means"
  calculate(stat = "diff in means", order = c("yes", "no"))

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

null_dist %>% visualize() +
  shade_p_value(obs_stat = obs_diff, direction = "two-sided")

null_dist %>%
  get_p_value(obs_stat = obs_diff, direction = "two_sided")

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