HAZAL GUNDUZ

Inference for numerical data

Getting Started

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.0     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(openintro)
## Loading required package: airports
## Loading required package: cherryblossom
## Loading required package: usdata
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.

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

Exercise 1. What are the cases in this data set? How many cases are there 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"…

=> There at 13,583 rows and also the number of cases in this sample. There are 13 cases in this data set.

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

Exercise 2. How many observations are we missing weights from?

sum(is.na(yrbss))
## [1] 9476

=> There are 9,476 missing.

The observations of weights we have 1004 missing values.

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

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?

yrbss2 <- yrbss %>%
  mutate(physical_3plus = ifelse(yrbss$physically_active_7d > 2, "yes", "no")) %>%
  na.exclude()
ggplot(yrbss2, aes(x = weight, y = physical_3plus)) + geom_boxplot() + theme_bw()

=> There is a relationship between physical activity and weight. We see that the weights are pretty similar, but there is a higher concentration of weight measures clustered together for those who exercise than those who don’t. For those that don’t exercise, we see more outliers in weight. The data is more normally distributed for those who exercise than for those who don’t, but both sets of data appears to be normally distributed.

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

Inference

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

=> There two conditions - independence and normality. According to the data, representative samples across national, state, tribal systems and is independent. The sample size and distribution of the boxplots would help us verify normality. We can assume with a sample size well over 1000 and no particularly outliers that the condition is satisfied.

yrbss %>%
  group_by(physical_3plus) %>%
  summarise(freq = table(weight)) %>%
  summarise(n = sum(freq))
## `summarise()` has grouped output by 'physical_3plus'. You can override using the `.groups` argument.

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: Who are physically active 3 or more days per week have the same average weight as those who are not physically active 3 or more days per week.

=> HA: Who are physically active 3 or more days per week have a different average weight when compared to those who are not physically active 3 or more days per week.

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"))
## Warning: Removed 1219 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.

set.seed(999)
null_dist <- yrbss %>%
  specify(weight ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
## Warning: Removed 1219 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, which 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")
## Warning: F usually corresponds to right-tailed tests. Proceed with caution.

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.

=> This the standard workflow for performing hypothesis tests.

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.

Standard Deviation

yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(sd_weight = sd(weight, na.rm = TRUE))

Mean

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

Sample size: N

yrbss %>% 
  group_by(physical_3plus) %>% 
  summarise(freq = table(weight)) %>%
  summarise(n = sum(freq))
## `summarise()` has grouped output by 'physical_3plus'. You can override using the `.groups` argument.
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))

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

uci_not
## [1] 67.219
lci_not
## [1] 66.12878
u_ci
## [1] 68.80209
l_ci
## [1] 68.09485

More Practice

Exercise 8. Calculate a 95% confidence interval for the average height in meters (height) and interpret it in context.

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

x_h  <- mean(yrbss$height, na.rm = TRUE)
sd_h <- sd(yrbss$height, na.rm = TRUE)
n_h  <- yrbss %>% 
  summarise(freq = table(height)) %>%
  summarise(n = sum(freq, na.rm = TRUE))

u_h <- x_h + z*(sd_h/sqrt(n_h))
l_h <- x_h - z*(sd_h/sqrt(n_h))
u_h
l_h

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.

t1 <- 1.645
uh1 <- x_h + t1*(sd_h/sqrt(n_h))
lh1 <- x_h - t1*(sd_h/sqrt(n_h))
uh1
lh1

=> The new confidence interval is 1.689705 to 1.692777. There is a slight difference in these two confidence intervals.

r_95 <- (u_h - l_h)
r_90 <- (uh1 - lh1)
r_95
r_90

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.

=> HO: There is no difference in the average height of those who are physically active at least 3 days per week and those who are not.

=> HA: There is a difference.

obs_diff_hgt <- yrbss %>%
  specify(height ~ physical_3plus) %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
## Warning: Removed 1219 rows containing missing values.
set.seed(45698)
null_dist_hgt <- yrbss %>%
  specify(height ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))
## Warning: Removed 1219 rows containing missing values.
visualize(null_dist_hgt) + 
  shade_p_value(obs_stat = obs_diff_hgt, direction = "two_sided")
## Warning: F usually corresponds to right-tailed tests. Proceed with caution.

null_dist_hgt %>%
  get_p_value(obs_stat = obs_diff_hgt, 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.
x_t <- 1.6665
n_t <- 4022
s_t <- 0.1029
x_yt <- 1.7032
n_yt <- 8342
s_yt <- 0.1033

z = 1.96

ut <- x_t + z*(s_t/sqrt(n_t))
lt <- x_t - z*(s_t/sqrt(n_t))
ut
## [1] 1.66968
lt
## [1] 1.66332
uyt <- x_yt + z*(s_yt/sqrt(n_yt))
lyt <- x_yt - z*(s_yt/sqrt(n_yt))
uyt
## [1] 1.705417
lyt
## [1] 1.700983

With 95% confident that the average height of students who are physically active at least 3 days per week is between 1.705 and 1.701 and the average height of students who are not physically active at least 3 days per week is between 1.670 and 1.663.

11.Exercise. 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())

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.

  1. Are there evidence that students who are heavier than the mean weight sleep more than students who weight less than the mean weight?

=> There is a relationship between weight and sleep, and no relationship between weight and sleep.

yrbss <- yrbss %>%
  mutate(sleep_less = ifelse(yrbss$school_night_hours_sleep < 7, "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()
boxplot(weight_less$weight, weight_more$weight,
        names = c("less_sleep", "more_sleep"))

mn <- mean(weight_less$weight)
sd <- sd(weight_less$weight)
max <- max(weight_less$weight)
max
## [1] 180.99
mn1 <- mean(weight_more$weight)
sd2 <- sd(weight_more$weight)
max2 <- max(weight_more$weight)
meandiff <- mn1 - mn
sd <- 
  sqrt(
  ((mn1^2) / nrow(weight_more)) +
  ((mn^2) / nrow(weight_less))
  )
df <- 2492-1
t <- qt(.05/2, df, lower.tail = FALSE)

upper_ci <- meandiff + t * sd
lower_ci <- meandiff - t * sd

c(lower_ci ,upper_ci)
## [1] -4.018176  1.009936
p_value <- 2*pt(t,df, lower.tail = FALSE)
p_value
## [1] 0.05

=> The null hypothesis can be rejected because p-value is 0.05..