Inference for numerical data

Author

Amanda Rose Knudsen

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', package='openintro')

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

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

There are 13,583 cases in the YRBSS dataset. A case is an observation here, and the help file found with running ?yrbss specifies that it is “A data frame with 13583 observations on the following 13 variables.” Additionally when running glimpse(yrbss) we see in the output that the number of rows is specified: 13,583. Each row (each case) represents an individual student who participated in the study. So there are 13,583 students represented in the YRBSS sample.

Remember that you can answer this question by viewing the data in the data viewer or by using the following command:

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

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.

yrbss_weight_summary <- summary(yrbss$weight)

yrbss_weight_summary

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  29.94   56.25   64.41   67.91   76.20  180.99    1004

ggplot(yrbss, aes(x = weight, na.rm = TRUE)) +
  geom_histogram(binwidth = 10,fill = "navy", color = "darkgrey") +
  labs(title = "Distribution of weight in kilograms in YRBSS")

weight_summary <- yrbss |> 
  summarise(
    mean_weight = mean(weight, na.rm = TRUE),
    median_weight = median(weight, na.rm = TRUE),
    sd_weight = sd(weight, na.rm = TRUE)
  )

weight_summary

# A tibble: 1 × 3
  mean_weight median_weight sd_weight
        <dbl>         <dbl>     <dbl>
1        67.9          64.4      16.9

In the weight_summary we can see that the minimum value is 29.94, the median is 64.41, the mean is 67.92, and there are 1004 NA’s (missing values). As we can see in the plot and with these summary statistics, the distribution is right skewed: the values are concentrated on the left wide of the plot.

How many observations are we missing weights from?

We can see that there are 1004 observations where we are missing 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"))

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?

yrbss_dropna_weight_physical3plus <- yrbss |> 
  drop_na(weight) |> 
  drop_na(physical_3plus)

ggplot(yrbss_dropna_weight_physical3plus, aes(x = physical_3plus, y = weight)) +
  geom_boxplot() +
  labs(x = "physically active 3+ days", y = "weight", 
       title = "Weight by physical activity levels boxplots") +
  theme_minimal()

The side-by-side boxplots show the medians are very close, as are the IQRs. We can see that the students who are physically active for 3+ days per week (“yes”) have a slightly higher median weight compared with those who are not physically active for 3+ days per week (“no”). At first this surprised me because I expected there to be lower weight overall for those who are physically fit, but then I realized that ‘working out’ or being physically active leads to a higher weight due to muscle mass increase. The ‘higher weights’ from not being physically fit are not related to muscle but rather to fat, and so this doesn’t surprise me that we don’t see too much of a difference in the median weights and IQRs, since “weight” alone does not tell the whole story of “fitness” or “health”.

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

# A tibble: 3 × 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

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

When we use the ‘original’ sample dataset YRBSS…

yrbss |> 
  group_by(physical_3plus) |> 
  summarise(
    n = n(),
    mean_weight = mean(weight, na.rm = TRUE),
    sd_weight = sd(weight, na.rm = TRUE)
    )

# A tibble: 3 × 4
  physical_3plus     n mean_weight sd_weight
  <chr>          <int>       <dbl>     <dbl>
1 no              4404        66.7      17.6
2 yes             8906        68.4      16.5
3 <NA>             273        69.9      17.6

When we use the ‘modified’ sample dataset YRBSS which has the NA values removed …

yrbss_dropna_weight_physical3plus |> 
  group_by(physical_3plus) |> 
  summarise(
    n = n(),
    mean_weight = mean(weight, na.rm = TRUE),
    sd_weight = sd(weight, na.rm = TRUE)
    )

# A tibble: 2 × 4
  physical_3plus     n mean_weight sd_weight
  <chr>          <int>       <dbl>     <dbl>
1 no              4022        66.7      17.6
2 yes             8342        68.4      16.5

Yes, conditions for inference are satisfied. - We know this is from a simple random sample from the CDC (condition that the data is from a random sample) - We know that the sample size (over 13k) is more than the n greater than or equal to 30 sample size required to be considered to be approximately normal (condition for normality) - We know that it would not be possible for an individual to exist in both groups (grouping is based on distinct individuals) and this sample size represents less than 10% of the population (condition that the observations should be independent within each group and between each group)

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: Highschoolers who are physically active 3+ days per week have the same average weight as those who are not physically active 3+ days per week. Ha: Highschoolers who are physically active 3+ days per week have a different average weight than those who are not physically active 3+ 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 %>%
  drop_na(physical_3plus) %>%
  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.

null_dist <- yrbss %>%
  drop_na(physical_3plus) %>%
  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, whichis 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()

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

obs_diff

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

min(null_dist$stat)

[1] -0.8859871

max(null_dist$stat)

[1] 1.02176

It appears that none (0) of the null permutations have a difference of at least obs_stat. The obs_stat (obs_diff) is 1.775, visualized with the red line in the plot above. None of the null permutations have a difference of at least this value. So I think this means that with an `obs_stat (obs_diff) is 1.775 we can reject the null in favor of the alternative hpyothesis.

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

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

This the standard workflow for performing hypothesis tests.

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.

I’m going to construct a confidence interval for the difference between the weights of those who exercise at least three times a week and those who don’t using two different approaches: firts, using infer, followed by the methods outlined in our textbook.

set.seed(58)
ci_dist <- yrbss  |> 
  drop_na(physical_3plus) |> 
  specify(weight ~ physical_3plus)  |> 
  generate(reps = 1000, type = "bootstrap")  |> 
  calculate(stat = "diff in means", order = c("yes", "no")) |> 
  get_confidence_interval(point_estimate = obs_diff, level = .95, type = "se")

ci_dist

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1     1.09     2.46

yrbss |> 
  drop_na(physical_3plus, weight) |> 
  group_by(physical_3plus)  |> 
  summarise(mean_weight = mean(weight, na.rm = TRUE),
            sd_weight = sd(weight, na.rm = TRUE),
            size_physical_3plus = n())

# A tibble: 2 × 4
  physical_3plus mean_weight sd_weight size_physical_3plus
  <chr>                <dbl>     <dbl>               <int>
1 no                    66.7      17.6                4022
2 yes                   68.4      16.5                8342

no_exercise3x_mean <- 66.7
no_exercise3x_sd <- 17.6
no_exercise3x_n <- 4022

yes_exercise3x_mean <- 68.4
yes_exercise3x_sd <- 16.5
yes_exercise3x_n <- 8342

point_estimate_difference <- yes_exercise3x_mean - no_exercise3x_mean
point_estimate_difference

[1] 1.7

se_difference <- sqrt( 
                      ((yes_exercise3x_sd^2) / yes_exercise3x_n) 
                       + ((no_exercise3x_sd^2) / no_exercise3x_n)
                      )

se_difference

[1] 0.3311381

Based on guidance in the textbook (page 267-268)… we can use the smaller of n1-1 and n2-2 for degrees of freedom if using a t-table to find tail areas. So I will use the smaller of n which is no_exercise3x_n equal to 4022. So the df is 4022-1.

t_critical_4021 <- qt(0.975, df = 4021)
t_critical_4021

[1] 1.960554

ci_lower_difference <- 
  point_estimate_difference - t_critical_4021 * se_difference
ci_upper_difference <- 
  point_estimate_difference + t_critical_4021 * se_difference

ci_difference <- c(ci_lower_difference, ci_upper_difference)
ci_difference

[1] 1.050786 2.349214

The 95% confidence interval using infer and bootstrapping is from 1.09 to 2.46. In the context of the data, we can say that we are 95% confidence that the true difference in means weights between highschoolers who do/don’t exercise 3x per week is between 1.10 and 2.45.

The 95% confidence interval using the method from the book is from 1.05 to 2.35. In the context of the data, we can say that we are 95% confident that the true difference in mean weights between highschoolers who do/don’t exercise 3x per week is between 1.05 and 2.35.

More Practice

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

yrbss  |> 
  drop_na(height) |> 
  summarise(mean_height = mean(height, na.rm = TRUE),
            sd_height = sd(height, na.rm = TRUE),
            n_height = n(),
            se_height = (sd_height) / (sqrt(n_height)),
            t_critical_height = qt(0.975, df = (n_height - 1))
            )

# A tibble: 1 × 5
  mean_height sd_height n_height se_height t_critical_height
        <dbl>     <dbl>    <int>     <dbl>             <dbl>
1        1.69     0.105    12579  0.000933              1.96

lower_ci_height <- 1.69 - 1.96 * 0.000933
upper_ci_height <- 1.69 + 1.96 * 0.000933

height_ci95 <- c(lower_ci_height, upper_ci_height)

height_ci95

[1] 1.688171 1.691829

We are 95% confident that the true population mean for the average height of highschoolers is between 1.688171 meters and 1.691829 meters.

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.

yrbss  |> 
  drop_na(height) |> 
  summarise(mean_height = mean(height, na.rm = TRUE),
            sd_height = sd(height, na.rm = TRUE),
            n_height = n(),
            se_height = (sd_height) / (sqrt(n_height)),
            t_critical_height_90 = qt(0.95, df = (n_height - 1))
            )

# A tibble: 1 × 5
  mean_height sd_height n_height se_height t_critical_height_90
        <dbl>     <dbl>    <int>     <dbl>                <dbl>
1        1.69     0.105    12579  0.000933                 1.64

lower_ci_height_90 <- 1.69 - 1.64 * 0.000933
upper_ci_height_90 <- 1.69 + 1.64 * 0.000933

height_ci90 <- c(lower_ci_height_90, upper_ci_height_90)

height_ci90

[1] 1.68847 1.69153

Using a different t critical value of 1.64 for a 90% confidence interval, we get a slightly different confidence interval (with otherwise the same inputs in the ‘finding a confidence interval for the mean’, with the formula again with guidance from page 256 of the textbook). This time, for a 90% confidence interval we see the values from 1.68847 to 1.69153. We are 90% confident that the true population mean for the average height of highschoolers is between 1.68847 and 1.69153, representing a narrower interval than the 95% confidence interval we obtained in the previous question. It makes sense that we have a narrower band with a lower confidence level – more precision but potentially lower accuracy. This goes along with our general understanding of how confidence levels work.

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.

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

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

diffheight3daysactive_pvalue <- diffheight3daysactive_nulldist  |> 
  get_p_value(obs_stat = diffheight3daysactive_obs, direction = "two_sided")

set.seed(58)
ci_dist_height <- yrbss  |> 
  drop_na(physical_3plus, height)|> 
  specify(height ~ physical_3plus)  |> 
  generate(reps = 1000, type = "bootstrap")  |> 
  calculate(stat = "diff in means", order = c("yes", "no")) |> 
  get_confidence_interval(point_estimate = diffheight3daysactive_obs, level = .95, type = "se")

ci_dist_height

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1   0.0336   0.0416

The confidence interval based on the above bootstrapped sampling distribution is from 0.0336 to 0.0416

yrbss |> 
  drop_na(physical_3plus, height) |> 
  group_by(physical_3plus)  |> 
  summarise(mean_height = mean(height, na.rm = TRUE),
            sd_height = sd(height, na.rm = TRUE),
            size_physical_3plus = n())

# A tibble: 2 × 4
  physical_3plus mean_height sd_height size_physical_3plus
  <chr>                <dbl>     <dbl>               <int>
1 no                    1.67     0.103                4022
2 yes                   1.70     0.103                8342

no_exercise3x_mean_height <- 1.67
no_exercise3x_sd_height <- 0.103
no_exercise3x_n_height <- 4022

yes_exercise3x_mean_height <- 1.70
yes_exercise3x_sd_height <- 0.103
yes_exercise3x_n_height <- 8342

point_estimate_difference_height <- 
  yes_exercise3x_mean_height - no_exercise3x_mean_height
point_estimate_difference_height

[1] 0.03

se_difference_height <- sqrt( 
                      ((yes_exercise3x_sd_height^2) / yes_exercise3x_n_height) 
                       + ((no_exercise3x_sd_height^2) / no_exercise3x_n_height)
                      )

se_difference_height

[1] 0.001977246

Like with the similar question above, using guidance in the textbook (page 267-268)… we can use the smaller of n1-1 and n2-2 for degrees of freedom if using a t-table to find tail areas. So I will use the smaller of n which is no_exercise3x_n equal to 4022. So the df is 4022-1.

We can use the same defined variable as above, t_critical_4021

ci_lower_difference_height <- 
  point_estimate_difference_height - t_critical_4021 * se_difference_height
ci_upper_difference_height <- 
  point_estimate_difference_height + t_critical_4021 * se_difference_height

ci_difference_height <- c(ci_lower_difference_height, ci_upper_difference_height)
ci_difference_height

[1] 0.0261235 0.0338765

Using methods outlined in the textbook, the 95% confidence interval is from 0.0261235 to 0.0338765. In the context of the data, we can say that we are 95% confident that the true difference in mean heights between highschoolers who do/don’t exercise 3x per week is between 0.026 meters and 0.034 meters.

The confidence interval based on the above bootstrapped sampling distribution is from 0.0336 to 0.0416. In the context of the data, we can say that we are 95% confident that the true difference in mean hights between highschoolers who do/don’t exercise 3x per week is between 0.033 meters and 0.0416 meters.

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

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

There are 7 different non-null options for the tv_hours_per_school_day variable. If we count NA (null, no answer) as an option, there are 8 different ‘options’.

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.

yrbss |> 
  group_by(school_night_hours_sleep) |> 
  summarise(n = n())

# A tibble: 8 × 2
  school_night_hours_sleep     n
  <chr>                    <int>
1 10+                        316
2 5                         1480
3 6                         2658
4 7                         3461
5 8                         2692
6 9                          763
7 <5                         965
8 <NA>                      1248

yrbss <- yrbss |> 
  mutate(
    sleep_ind = if_else(school_night_hours_sleep %in% c("10+", "9", "8"), "yes", "no")
  )
yrbss

# A tibble: 13,583 × 15
     age gender grade hispanic race                     height weight helmet_12m
   <int> <chr>  <chr> <chr>    <chr>                     <dbl>  <dbl> <chr>     
 1    14 female 9     not      Black or African Americ…  NA      NA   never     
 2    14 female 9     not      Black or African Americ…  NA      NA   never     
 3    15 female 9     hispanic Native Hawaiian or Othe…   1.73   84.4 never     
 4    15 female 9     not      Black or African Americ…   1.6    55.8 never     
 5    15 female 9     not      Black or African Americ…   1.5    46.7 did not r…
 6    15 female 9     not      Black or African Americ…   1.57   67.1 did not r…
 7    15 female 9     not      Black or African Americ…   1.65  132.  did not r…
 8    14 male   9     not      Black or African Americ…   1.88   71.2 never     
 9    15 male   9     not      Black or African Americ…   1.75   63.5 never     
10    15 male   10    not      Black or African Americ…   1.37   97.1 did not r…
# ℹ 13,573 more rows
# ℹ 7 more variables: text_while_driving_30d <chr>, physically_active_7d <int>,
#   hours_tv_per_school_day <chr>, strength_training_7d <int>,
#   school_night_hours_sleep <chr>, physical_3plus <chr>, sleep_ind <chr>

yrbss |> 
  group_by(sleep_ind) |> 
  summarise(n = n())

# A tibble: 2 × 2
  sleep_ind     n
  <chr>     <int>
1 no         9812
2 yes        3771

Question: Is there a relationship between the average weight of highschoolers and if they get 8+ hours of sleep per night on school nights? We will see if there is a relationship between the average weight of highschoolers and if they get 8+ hours of sleep per night on school nights. We will see if there is a difference in the average weights of highschoolers for those who do get 8+ hours per sleep on school nights and those who do not get 8+ hours per sleep on school nights.

Conditions for inference are satisfied. - We know this is from a simple random sample from the CDC (condition that the data is from a random sample) - We know that the sample size (over 13k) is more than the n greater than or equal to 30 sample size required to be considered to be approximately normal (condition for normality) - We know that it would not be possible for an individual to exist in both groups (grouping is based on distinct individuals) and this sample size represents less than 10% of the population (condition that the observations should be independent within each group and between each group)

Our null hypothesis was a p-value of 0.006 which is lower than our alpha of 0.05 but also is not zero like we’d seen previously in our other tests of no difference.

weightsleep_obs <- yrbss  |> 
  drop_na(weight, school_night_hours_sleep) |> 
  specify(weight ~ sleep_ind)  |> 
  calculate(stat = "diff in means", order = c("yes", "no"))

weightsleep_null <- yrbss  |> 
  drop_na(weight, school_night_hours_sleep) |> 
  specify(weight ~ sleep_ind) |> 
  hypothesize(null = "independence")  |> 
  generate(reps = 1000, type = "permute")  |> 
  calculate(stat =  "diff in means", order = c("yes", "no"))

weightsleep_pvalue <- weightsleep_null  |> 
  get_p_value(obs_stat = weightsleep_obs, direction = "two_sided")

visualize(weightsleep_null) +  
  shade_p_value(obs_stat = obs_diff, direction = "two_sided")

set.seed(58)
ci_dist_sleep <- yrbss  |> 
  drop_na(weight, school_night_hours_sleep) |> 
  specify(weight ~ sleep_ind)  |> 
  generate(reps = 1000, type = "bootstrap")  |> 
  calculate(stat = "diff in means", order = c("yes", "no")) |> 
  get_confidence_interval(
    point_estimate = weightsleep_obs, level = .95, type = "se")

ci_dist_sleep

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1    -1.66   -0.306

Our 95% confidence interval indicates that we can be 95% confident that the true difference in weight between those who get 8+ hours of sleep per night and those who don’t is a difference between -1.66 and -0.306. Our p-value for one side is 0.06 (doubled per guidance in the textbook would be 0.12 which is greater than our alpha of 0.05, so we would not reject the null.)

# weightsleep_pvalue
min(weightsleep_null$stat)

[1] -1.067403

max(weightsleep_null$stat)

[1] 1.268649

weightsleep_obs

Response: weight (numeric)
Explanatory: sleep_ind (factor)
# A tibble: 1 × 1
    stat
   <dbl>
1 -0.984