infer.Rmd

# Statistical inference with infer

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  out.width = "100%"
)
suppressPackageStartupMessages(library(dplyr))
suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(readr))
suppressPackageStartupMessages(library(infer))
county_complete <- read_rds("data/county_complete.rds")
```

## Case study: Comparing work travel times

The following is a worked example using the county complete dataset from Section \@ref(representing-distributions) that demonstrates how to use the `infer` package to test if the difference in mean values of two distributions is statistically significant.
To follow along, load the following libraries and dataset into your R session:

```r
# Packages
library(dplyr)
library(ggplot2)
library(readr)
library(infer)
# Dataset
county_complete <- read_rds(
  url("http://data.cds101.com/county_complete.rds")
)
```

### Work travel times in Iowa and Nebraska

In Section \@ref(representing-distributions), we looked at the distribution of average work travel times in Iowa and Nebraska measured on the county level.

```{r ia-ne-travel-times-pmf, echo = FALSE}
county_complete %>%
  filter(state == "Iowa" | state == "Nebraska") %>%
  ggplot() +
  geom_histogram(
    mapping = aes(
      x = mean_work_travel,
      y = ..density..,
      fill = state,
      color = state
    ),
    position = "identity",
    alpha = 0.5,
    binwidth = 1
  ) +
  labs(
    title = "Average work travel times in Iowa and Nebraska",
    x = "Average work travel time (min)", y = "PMF"
  )
```

In our visual inspection of the probability mass function, we concluded that the distribution of work travel times is shifted slightly towards shorter commutes in Nebraska versus Iowa.
The summary statistics for the two distributions seem to support this conclusion:

```{r ia-ne-summary-statistics, results = "hide"}
# Filter dataset to only include counties in Iowa and Nebraska
ia_ne_county <- county_complete %>%
  filter(state == "Iowa" | state == "Nebraska")

# Compute summary statistics of mean_work_travel column for Iowa and Nebraska
ia_ne_summary_stats <- ia_ne_county %>%
  group_by(state) %>%
  summarize(
    mean = mean(mean_work_travel),
    median = median(mean_work_travel),
    sd = sd(mean_work_travel),
    iqr = IQR(mean_work_travel)
  )
```

```{r ia-ne-summary-statistics-table, echo = FALSE}
ia_ne_summary_stats %>%
  knitr::kable() %>% #format = "html") %>%
  kableExtra::kable_styling(
    bootstrap_options = combine("striped", "responsive"),
    full_width = FALSE
  )
```

Let's compute the difference in means directly using `ia_ne_summary_stats`.
First, we need to grab the `mean` column:

```{r ia-ne-get-means}
ia_ne_means <- ia_ne_summary_stats %>%
  pull(mean)
```

`ia_ne_means` is a simple vector that stores two numbers.

```{r ia-ne-means-output}
ia_ne_means
```

The first number of the vector, which we can access with `ia_ne_means[1]`, is the mean work travel times in Iowa.
The second number of the vector, which we can access with `ia_ne_means[2]`, is the mean work travel times in Nebraska.
Knowing this, we easily compute the difference in means,

```{r ia-ne-diff-in-means}
ia_ne_diff_in_means <- ia_ne_means[1] - ia_ne_means[2]
```

obtaining the following value of `ia_ne_diff_in_means`:

```{r ia-ne-diff-in-means-output, echo = FALSE, comment = NA}
cat(
  glue::glue(
    "# (mean work travel time in Iowa) - ",
    "(mean work travel time in Nebraska)\n",
    "ia_ne_diff_in_means = {round(ia_ne_diff_in_means, 2)} minutes\n"
  )
)
```

However, only checking the numerical difference between the mean values is not enough if we want to claim there is a **statistically significant** difference in commute times between Iowa and Nebraska.
There is significant overlap between the two distributions, with the standard deviation (sd) of each distribution being 3--4 times *larger* than the difference in means.
Before we can conclude that the difference in means is statistically significant, we need to estimate how likely it is that this difference could have arisen from chance alone.

### Defining the hypothesis test

To determine whether or not there is a meaningful difference in means between the two distributions, we will conduct a two-sided hypothesis test using the following null and alternative hypotheses:

*   **Null hypothesis**: There is no difference between the mean work travel time in Iowa and the mean work travel time in Nebraska.

    $$\text{H}_{0}:\mu_{\text{IA}}-\mu_{\text{NE}}=0$$

*   **Alternative hypothesis**: There is a difference between the mean work travel time in Iowa and the mean work travel time in Nebraska.
    
    $$\text{H}_{\text{A}}:\mu_{\text{IA}}-\mu_{\text{NE}}\neq{}0$$

We set our significance level to $\alpha=0.05$, which we will use as a cutoff for determining whether or not a result is statistically significant.

### Building the null distribution

To conduct our hypothesis test, we first must generate our null distribution using the `infer` package:

```{r ia-ne-mean-work-travel-null}
ia_ne_mean_work_travel_null <- ia_ne_county %>%
  specify(formula = mean_work_travel ~ state) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 10000, type = "permute") %>%
  calculate(stat = "diff in means", order = combine("Iowa", "Nebraska"))
```

Let's explain what each of these lines is doing.
We start by piping our filtered dataset stored in `ia_ne_county` into `specify`:

```r
specify(formula = mean_work_travel ~ state)
```

Here we *specify* the **response** variable and the **explanatory** variable in the `formula =` input keyword.
The formula syntax is given as follows:

```
response ~ explanatory
```

Therefore, in our formula `mean_work_travel ~ state`, `mean_work_travel` is the response variable and `state` is the explanatory variable.
How do we know to use `mean_work_travel` as the response variable?
It's because of how we formulated our null and alternative hypotheses, such that we are testing whether a person's mean work travel time is affected by the state he or she lives in, **not** whether the state that a person lives in is affected by his or her mean work travel time.
Finally, because the **response variable** is numerical, we don't use the `success =` keyword in `specify`.

The next line that we pipe into is:

```r
hypothesize(null = "independence")
```

Here we indicate that we are conducting a hypothesis test to check whether or not the variables provided in `specify` are independent of one another.
A basic rule of thumb is, if you specify **both** an explanatory and a response variable in the formula of `specify`, then you will use `null = "independence"`.

Next we pipe into this line:

```r
generate(reps = 10000, type = "permute")
```

Here we indicate that we want to *permute* the `mean_work_travel` and `state` columns 10,000 times to build up our null distribution.
In effect, we are randomly shuffling the rows in `mean_work_travel` into a different order, and then randomly shuffling the rows in `state` into yet another order, which should remove any connections between the `mean_work_travel` and `state` columns (this is what we mean by null distribution).
We know to use `type = "permute"` because we have specified both a response and an explanatory variable and we want to see what happens when we assume the two columns are independent.
As for `reps = 10000`, this controls the overall accuracy of our hypothesis tests.
There is no hard and fast rule for the value to set, but if you want more accurate results, increase `reps`.
However, keep in mind that the larger `reps` is, the *longer* it will take to complete the calculation.

We finish up by piping into the last line:

```r
calculate(stat = "diff in means", order = combine("Iowa", "Nebraska"))
```

Here we indicate that we are comparing the difference in means between our two distributions (the mean work travel times in Iowa and Nebraska).
The input `order = combine("Iowa", "Nebraska")` defines the subtraction order, so we will find the difference in means by subtracting the mean time in Iowa from the mean time in Nebraska, which is the same order we used for `ia_ne_diff_in_means <- ia_ne_means[1] - ia_ne_means[2]`.
As for the `stat` input, we know to use `stat = "diff in means"` because of how we formulated our null and alternative hypotheses, i.e. we are comparing the **difference** in **mean** work travel times between Iowa and Nebraska.
The way we know that we **can't** use something like `stat = "diff in props"` is because at least one of our variables in `specify` is numerical (`mean_work_travel`).
The input `stat = "diff in props"` only applies if **both** the explanatory and response variables are categorical.

The null distribution that we generated looks as follows:

```{r ia-ne-diff-in-means-null-distribution}
ia_ne_mean_work_travel_null %>%
  visualize() +
  labs(
    x = "difference in means",
    title = "Difference in mean work travel times null distribution"
  )
```

### Computing the two-sided p-value

After we've generated our null distribution, we can compute the *p*-value using `get_p_value`.
To compute a *p*-value, we need an **observed statistic**, which is simply the difference in means we computed using the dataset itself.
This is what we computed using `ia_ne_diff_in_means <- ia_ne_means[1] - ia_ne_means[2]`.
However, `infer` also provides a shortcut for computing the observed statistic so that you don't have to use `summarize()`,

```{r obs-stat-using-infer}
ia_ne_diff_in_means <- ia_ne_county %>%
  specify(formula = mean_work_travel ~ state) %>%
  calculate(stat = "diff in means", order = combine("Iowa", "Nebraska"))
```

Note that all we've done is taken the code we used to generate the null distribution and removed the `hypothesize()` and `generate()` lines.
Now that we have our observed statistic, we can compute the *p*-value for our **two-sided hypothesis test**,

```{r ia-ne-p-value-two-sided}
ia_ne_p_value_two_sided <- ia_ne_mean_work_travel_null %>%
  get_p_value(obs_stat = ia_ne_diff_in_means, direction = "both")
```

```{r ia-ne-p-value-two-sided-print, echo = FALSE}
ia_ne_p_value_two_sided %>%
  knitr::kable() %>% #format = "html") %>%
  kableExtra::kable_styling(
    bootstrap_options = combine("striped", "responsive"),
    full_width = FALSE
  )
```

```{r summarize-two-sided-test-results, echo = FALSE, results = "asis"}
cat(
  glue::glue(
    "Our two-sided p-value is {ia_ne_p_value_two_sided}, which is ",
    "{if_else(ia_ne_p_value_two_sided < 0.05, 'below', 'above')} ",
    "our significance level of $\\alpha=0.05$. We therefore ",
    "{if_else(ia_ne_p_value_two_sided < 0.05, 'reject', 'are unable to reject')} ",
    "the null hypothesis in favor of the alternative hypothesis. ",
    "The difference between the mean work travel times in Iowa and Nebraska ",
    "{if_else(ia_ne_p_value_two_sided < 0.05, 'is', 'is not')} ",
    "statistically significant."
  )
)
```

We visualize the meaning of the *p*-value in the following way:

```{r ia-ne-diff-in-means-null-distribution-p-values}
ia_ne_mean_work_travel_null %>%
  visualize() +
  shade_p_value(obs_stat = ia_ne_diff_in_means, direction = "both") +
  labs(
    x = "difference in means",
    title = "Difference in mean work travel times null distribution"
  )
```

The portions of the null distribution that are in the shaded red area correspond to results that are **more extreme** than the observed value `r round(ia_ne_diff_in_means, 2)` in a two-sided hypothesis test.
The more that a null distribution lies within the red portions of the visualization, the more likely it becomes that we will fail to reject the null distribution and the difference in means that we observed in the dataset arose due to random chance alone.

### Computing the 95% confidence interval

In addition to the two-sided hypothesis test, it is also good practice to compute the 95% confidence interval for the difference in mean work travel times between Iowa and Nebraska.
To compute this, we need to generate the **bootstrap distribution** for the difference in means:

```{r ia-ne-mean-work-travel-bootstrap}
ia_ne_mean_work_travel_bootstrap <- ia_ne_county %>%
  specify(formula = mean_work_travel ~ state) %>%
  generate(reps = 10000, type = "bootstrap") %>%
  calculate(stat = "diff in means", order = combine("Iowa", "Nebraska"))
```

You'll notice that the code to generate the bootstrap distribution is pretty similar to the code used to generate the null distribution.
All we had to do was remove the `hypothesize(null = "independence")` line and change the `type =` keyword in `generate` from `"permute"` to `"bootstrap"`, and that's it!

You might be wondering, how are the methods for generating the null distribution and bootstrap distribution different?
For this example, they differ in two important ways:

1.  Unlike when we generated the null distribution, the `mean_work_travel` and `state` columns are **not** randomly shuffled

2.  Generating a bootstrap distribution requires *sampling with replacement*, meaning we grab different rows in `ia_ne_county` at random (we might even grab the same row more than), which we do not do if we permute the columns

The bootstrap distribution we just generated looks as follows:

```{r ia-ne-mean-work-travel-bootstrap-visualization}
ia_ne_mean_work_travel_bootstrap %>%
  visualize() +
  labs(
    x = "difference in means",
    title = "Difference in mean work travel times bootstrap distribution"
  )

```

This shows what could happen if we went out and collected another set of samples of the data in this dataset.
Due to random variations in the people we survey, sometimes we might end up with a dataset where the difference between the mean work travel times in Iowa and Nebraska is closer to 2 minutes instead of 1 minute, or it might be closer to zero.
Most of the time, it will be close to 1 minute.

The 95% confidence interval is defined as the fraction of the data in the bootstrap distribution that lies between the 2.5th and 97.5th percentiles:

```{r ia-ne-mean-work-travel-ci}
ia_ne_ci95 <- ia_ne_mean_work_travel_bootstrap %>%
  get_confidence_interval()
```

```{r ia-ne-mean-work-travel-ci-print, echo = FALSE}
ia_ne_ci95 %>%
  knitr::kable() %>% #format = "html") %>%
  kableExtra::kable_styling(
    bootstrap_options = combine("striped", "responsive"),
    full_width = FALSE
  )
```

This means that, if we keep collecting new samples of mean work travel times in Iowa and Nebraska, we expect that 95% of our samples will have a difference in mean work travel times that lies between `r ia_ne_ci95[1]` and `r ia_ne_ci95[2]`.
Also note that the lower bound of the confidence interval `r if_else(ia_ne_ci95[1] < 0, "intersects", "does not intersect")` with a difference of means equal to zero, lending further support to the conclusion we reached using our two-sided hypothesis test.

We close this example by visualizing the 95% confidence interval by shading the middle 95% of the data in the bootstrap distribution:

```{r ia-ne-mean-work-travel-bootstrap-visualization-ci95}
ia_ne_mean_work_travel_bootstrap %>%
  visualize() +
  shade_confidence_interval(endpoints = ia_ne_ci95) +
  labs(
    x = "difference in means",
    title = "Difference in mean work travel times 95% confidence interval"
  )

```