Inference for numerical data

Getting Started

Load packages

In this lab, we will explore and visualize the data using the tidyverse suite of packages, and perform statistical inference using infer. The data can be found in the companion package for OpenIntro resources, openintro.

Let’s load the packages.

library(tidyverse)
library(openintro)
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.

Load the yrbss data set into your workspace.

data(yrbss)

There are observations on 13 different variables, some categorical and some numerical. The meaning of each variable can be found by bringing up the help file:

?yrbss

## starting httpd help server ... done

Exercise 1

What are the cases in this data set? How many cases are there in our sample?

glimpse(yrbss)

They are 13,583 cases and case represents a high schooler (9th through 12th grade)

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?

We are missing 1004 observations from weights.

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?

# Create a subset for physical_3plus and weight

yrbss_phys_weight <- yrbss %>%
  filter(physical_3plus == "yes", weight != "NA")

yrbss_no_phys_weight <- yrbss %>%
  filter(physical_3plus == "no", weight != "NA")

boxplot(yrbss_phys_weight$weight, yrbss_no_phys_weight$weight,
        names = c("Weight for phys. act", "Weight for no phys. act"))

As the box plots show, there is not much of the difference between the median of the two variables (I also did the summary below to compare the mean but it is not significantly different neither).

There is not very clear relationship the two variables. My expectation was to see youth physically active to have a low average weight compared to no active ones but no a very significant difference because there are many factors that play an important role in someone’s weight.

summary(yrbss_no_phys_weight$weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   29.94   54.43   62.60   66.67   74.84  180.99

summary(yrbss_phys_weight$weight)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   33.11   56.70   65.77   68.45   77.11  160.12

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

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 3 x 2
##   physical_3plus mean_weight
##   <chr>                <dbl>
## 1 no                    66.7
## 2 yes                   68.4
## 3 <NA>                  69.9

There is an observed difference, but is this difference statistically significant? In order to answer this question we will conduct a hypothesis test.

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

yrbss %>%
  group_by(physical_3plus) %>%
  summarise(n_weight = n())

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 3 x 2
##   physical_3plus n_weight
##   <chr>             <int>
## 1 no                 4404
## 2 yes                8906
## 3 <NA>                273

Independence: The samples observations are high schoolers and they are independent each others. As we can see the count of “yes”, “no”, and “NA” give us the total.

Random samples: As stated in the “data”, it is a random sample of observations.

Approximately normal: more tan 30 samples, CLT.

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 for those who exercise at least times a week is the same to those who don’t.

H1: The average weights for those who exercise at least times a week differ to those who don’t.

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

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(1412)
null_dist <- yrbss %>%
  specify(weight ~ physical_3plus) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("yes", "no"))

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 is 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?

count_null <- null_dist %>%
  filter(stat >= obs_diff)

count_null

## # A tibble: 0 x 2
## # ... with 2 variables: replicate <int>, stat <dbl>

None of these null permutations have a difference of at least obs_stat

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

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.

# Find a point estimate 

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

## Warning: Removed 1219 rows containing missing values.

# Generate null distribution

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

# Confidence interval

ci_null_dist %>%
  
  get_confidence_interval(point_estimate = point_estimate,
                          
                          level = 0.95,
                          
                          type = "se")

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

We are 95% confident that the difference between the weights of those who exercise at least three times a week and those who don’t falls in (1.13, 2.41)

---

More Practice

Exercise 8

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

height_dist1 <- yrbss %>%
  specify(response = height) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean") %>%
  get_ci(level = 0.95)

## Warning: Removed 1004 rows containing missing values.

We are 95% confident that the average heights in meter falls in (1.689633, 1.693176)

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.

height_dist2 <- yrbss %>%
  specify(response = height) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "mean") %>%
  get_ci(level = 0.90)

## Warning: Removed 1004 rows containing missing values.

90% confident that the average heights in meter falls in (1.689685, 1.692724)

Although the difference of both interval is in order of micro, 90% confidence is still narrower (very small difference and can barely tell) than the one obtained in the previous exercise.

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.

Hypotheses

H0: The average height for those who exercise at least times a week is the same to those who don’t.

H1: The average height for those who exercise at least times a week differ to those who don’t.

# Find a point estimate 

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

## Warning: Removed 1219 rows containing missing values.

# Generate null distribution

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

# p_value

null_dist3 %>%
  
  get_p_value(obs_stat = point_estimate3,
              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

  # Visualize


null_dist3 %>%
  visualise() 

We reject the null hypothesis. The data provides evidence for the alternative hypothesis.

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 = n())

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 8 x 2
##   hours_tv_per_school_day     n
##   <chr>                   <int>
## 1 <1                       2168
## 2 1                        1750
## 3 2                        2705
## 4 3                        2139
## 5 4                        1048
## 6 5+                       1595
## 7 do not watch             1840
## 8 <NA>                      338

There are 7 options excluding “NA”

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 \(\alpha\) level, and conclude in context.

Conduct a hypothesis test evaluating whether the average height is different for those who sleep 10+ hours and the rest of the high schoolers.

Assumptions

Independence: The samples observations are high schoolers and they are independent each others.

Random samples: As stated in the “data”, it is a random sample of observations.

Approximately normal:more tan 30 samples, CLT.

Hypotheses

H0: The average height for those who sleep 10+ hours is the same to the rest of the high schoolers.

H1: The average height for those who sleep 10+ hours differ to the rest of the high schoolers.

Test

\(\alpha = 0.05\)

# Create a new var

yrbss <- yrbss %>%
  mutate(sleeper = ifelse(yrbss$school_night_hours_sleep == "10+", "yes", "no"))

# Check the mean for both cases

yrbss %>%
  group_by(sleeper) %>%
  summarise(mean_height = mean(height, na.rm = TRUE))

## `summarise()` ungrouping output (override with `.groups` argument)

## # A tibble: 3 x 2
##   sleeper mean_height
##   <chr>         <dbl>
## 1 no             1.69
## 2 yes            1.68
## 3 <NA>           1.70

# Find a point estimate 

point_estimate <- yrbss %>%
  specify(height ~ sleeper) %>%
  calculate(stat = "diff in means" , order = c("yes", "no")) 

## Warning: Removed 2102 rows containing missing values.

# Generate null distribution

set.seed(1412)

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

## Warning: Removed 2102 rows containing missing values.

# p_value

null_dist %>%
  
  get_p_value(obs_stat = point_estimate,
              direction = "two-sided")

## # A tibble: 1 x 1
##   p_value
##     <dbl>
## 1   0.086

# Confidence interval

null_dist %>%
  
  get_confidence_interval(point_estimate = point_estimate,
                          
                          level = 0.95,
                          
                          type = "se")

## # A tibble: 1 x 2
##   lower_ci upper_ci
##      <dbl>    <dbl>
## 1  -0.0249  0.00177

# Visualize


null_dist3 %>%
  visualise() 

Results

From the hypthesis test, we can see that the p_value is greater than 0.05. Thus we fail to reject the null hypothesis, the data don’t provide convincing evidence that the average height for those who sleep 10+ hours differ to the rest of the high schoolers.

Looking at the confidence interval, We are 95% confident that the diffrence of the average height for those who sleep 10+ hours and the rest of the high schoolers falls in (-0.0249, 0.0018).

Also, the difference of mean (mean_diff of yes and no is 1.68 - 1.69 = -0.01) is include to our 95% confidence interval founded.

Same conclusion from the hypothesis and the confidence interval.

---

This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.