---
title: "Tidy t-Tests with infer"
description: "Conducting t-Tests on tidy data with infer."
output: rmarkdown::html_vignette
vignette: |
  %\VignetteIndexEntry{t-Tests}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r settings, include=FALSE}
knitr::opts_chunk$set(fig.width = 6, fig.height = 4.5) 
options(digits = 4)
```

```{r load-packages, echo = FALSE, message = FALSE, warning = FALSE}
library(ggplot2)
library(dplyr)
devtools::load_all()
```

### Introduction

In this vignette, we'll walk through conducting t-tests using `infer`. We'll start out with a 1-sample t-test, which compares a sample mean to a hypothesized true mean value. Then, we'll discuss paired t-tests, which are a special use case of 1-sample t-tests, and evaluate whether differences in paired values (e.g. some measure taken of a person before and after an experiment) differ from 0. Finally, we'll wrap up with 2-sample t-tests, testing the difference in means of two populations using a sample of data drawn from them.

Throughout this vignette, we'll make use of the `gss` dataset supplied by `infer`, which contains a sample of data from the General Social Survey. See `?gss` for more information on the variables included and their source. Note that this data (and our examples on it) are for demonstration purposes only, and will not necessarily provide accurate estimates unless weighted properly. For these examples, let's suppose that this dataset is a representative sample of a population we want to learn about: American adults. The data looks like this:

```{r glimpse-gss-actual, warning = FALSE, message = FALSE}
dplyr::glimpse(gss)
```

### 1-Sample t-Test

The 1-sample t-test can be used to test whether a sample of continuous data could have plausibly come from a population with a specified mean. As an example, we'll test whether the average American adult works 40 hours a week using data from the `gss`. To do so, we make use of the `hours` variable, giving the number of hours that respondents reported having worked in the previous week. The distribution of `hours` in the observed data looks like this:

```{r plot-1-sample, echo = FALSE}
gss %>%
  ggplot2::ggplot() +
  ggplot2::aes(x = hours) +
  ggplot2::geom_histogram(bins = 20) +
  ggplot2::labs(x = "hours: Number of Hours Worked",
                y = "Number of Responses") +
  ggplot2::scale_x_continuous(breaks = seq(0, 90, 10))
```

Note the warning about missing values---many respondents' entries are missing. If we were actually carrying out this hypothesis test, we might look further into how this data was collected; it's possible that some of the missing values should actually be 0 hours.

In general, though, it looks like most respondents reported having worked 40 hours, but there's quite a bit of variability. Let's test whether we have evidence that the true mean number of hours that Americans work per week is 40.

First, to calculate the observed statistic, we can use `specify()` and `calculate()`.

```{r calc-obs-stat-1-sample, warning = FALSE, message = FALSE}
# calculate the observed statistic
observed_statistic <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  calculate(stat = "t")
```

The observed statistic is `r observed_statistic`. Now, we want to compare this statistic to a null distribution, generated under the assumption that the mean was actually 40, to get a sense of how likely it would be for us to see this observed statistic if the true number of hours worked per week in the population was really 40.

We can `generate` the null distribution using the bootstrap. In the bootstrap, for each replicate, a sample of size equal to the input sample size is drawn (with replacement) from the input sample data. This allows us to get a sense of how much variability we'd expect to see in the entire population so that we can then understand how unlikely our sample mean would be.

```{r generate-null-1-sample, warning = FALSE, message = FALSE}
# generate the null distribution
null_distribution_1_sample <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  generate(reps = 1000, type = "bootstrap") %>%
  calculate(stat = "t")
```

To get a sense for what these distributions look like, and where our observed statistic falls, we can use `visualize()`:

```{r visualize-1-sample, warning = FALSE, message = FALSE}
# visualize the null distribution and test statistic!
null_distribution_1_sample %>%
  visualize() + 
  shade_p_value(observed_statistic,
                direction = "two-sided")
```

It looks like our observed statistic of `r observed_statistic` would be relatively unlikely if the true mean was actually 40 hours a week. More exactly, we can calculate the p-value:

```{r p-value-1-sample, warning = FALSE, message = FALSE}
# calculate the p value from the test statistic and null distribution
p_value_1_sample <- null_distribution_1_sample %>%
  get_p_value(obs_stat = observed_statistic,
              direction = "two-sided")

p_value_1_sample
```

Thus, if the true mean number of hours worked per week was really 40, the probability that we would see a test statistic as or more extreme than `r observed_statistic` is approximately `r p_value_1_sample`.

Note that, equivalently to the steps shown above, the package supplies a wrapper function, `t_test`, to carry out 1-sample t-Tests on tidy data. The syntax looks like this:

```{r t-test-wrapper, message = FALSE, warning = FALSE}
t_test(gss, response = hours, mu = 40)
```

Also, to just calculate the observed statistic, the package supplies the wrapper `t_stat()`.

```{r t-stat-wrapper, message = FALSE, warning = FALSE}
t_stat(gss, response = hours, mu = 40)
```

### Paired t-Test

You might be interested in running a paired t-Test. Paired t-tests can be used in situations when there is a natural pairing between values in distributions---a common example would be two columns, `before` and `after`, say, that contain measurements from a patient before and after some treatment. To compare these two distributions, then, we're not necessarily interested in how the two distributions look different altogether, but how these two measurements from each individual change across time. (Pairings don't necessarily have to be over time; another common usage is measurements from two married people, for example.) Thus, we can create a new column (see `mutate()` from the `dplyr` package if you're not sure how to do this) that is the difference between the two: `difference = after - before`, and then examine _this_ distribution to see how each individuals' measurements changed over time.

Once we've `mutate()`d that new `difference` column, we can run a 1-sample t-Test on it, where our null hypothesis is that `mu = 0` (i.e. the difference between these measurements before and after treatment is, on average, 0). To do so, we'd use the procedure outlined in the above section.


### 2-Sample t-Test

2-Sample t-Tests evaluate the difference in mean values of two populations using data randomly-sampled from the population that approximately follows a normal distribution. As an example, we'll test if Americans work the same number of hours a week regardless of whether they have a college degree or not using data from the `gss`. The `college` and `hours` variables allow us to do so:

```{r plot-2-sample, echo = FALSE}
gss %>%
  ggplot2::ggplot() +
  ggplot2::aes(x = college, y = hours) +
  ggplot2::geom_boxplot() +
  ggplot2::labs(x = "college: Whether the Respondent has a College Degree",
                y = "hours: Number of Hours Worked")
```

It looks like both of these distributions are centered near 40 hours a week, but the distribution for those with a degree is slightly right skewed.

Again, note the warning about missing values---many respondents' values are missing. If we were actually carrying out this hypothesis test, we might look further into how this data was collected; it's possible that whether or not a value in either of these columns is missing is related to what that value would be. 

Now, to calculate the observed $t$ statistic, we can use `specify()` and `calculate()`.

```{r calc-obs-stat-2-sample, warning = FALSE, message = FALSE}
# calculate the observed statistic
observed_statistic <- gss %>%
  specify(hours ~ college) %>%
  calculate(stat = "t", order = c("degree", "no degree"))

observed_statistic
```

Note that, in the line `specify(hours ~ college)`, we could have swapped this out with the syntax `specify(response = hours, explanatory = college)`!

The `order` argument in that `calculate` line gives the order to subtract the mean values in: in our case, we're taking the mean number of hours worked by those with a degree minus the mean number of hours worked by those without a degree; a positive difference, then, would mean that people with degrees worked more than those without a degree.

Now, we want to compare this $t$-statistic to a null distribution, generated under the assumption that the number of hours worked a week has no relationship with whether or not one has a college degree, to get a sense of how likely it would be for us to see this $t$-statistic if there were really no relationship between these two variables.

We can `generate` the null distribution using permutation, where, for each replicate, each value of degree status will be randomly reassigned (without replacement) to a new number of hours worked per week in the sample in order to break any association between the two.

```{r generate-null-2-sample, warning = FALSE, message = FALSE}
# generate the null distribution with randomization
null_distribution_2_sample_permute <- gss %>%
  specify(hours ~ college) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "t", order = c("degree", "no degree"))
```

We can also generate the null distribution using the theoretical $t$ distribution. The code to do that looks like this:

```{r generate-null-2-sample-theoretical, warning = FALSE, message = FALSE}
# generate the null distribution with the theoretical t
null_distribution_2_sample_theoretical <- gss %>%
  specify(hours ~ college) %>%
  hypothesize(null = "independence") %>%
  # generate() isn't used for the theoretical version!
  calculate(stat = "t", order = c("degree", "no degree"))
```

Again, note that, in the lines `specify(hours ~ college)` in the above two chunks, we could have used the syntax `specify(response = hours, explanatory = college)` instead!

To get a sense for what these distributions look like, and where our observed statistic falls, we can use `visualize()`. First, to do this with the randomization-based distribution:

```{r visualize-2-sample, warning = FALSE, message = FALSE}
# visualize the randomization-based null distribution and test statistic!
null_distribution_2_sample_permute %>%
  visualize() + 
  shade_p_value(observed_statistic,
                direction = "two-sided")
```

The syntax to do the same thing with the theory-based null distribution looks similar---we just pipe in the theory-based null distribution instead, and drop in `method = "theoretical"` to `visualize()`.

```{r visualize-2-sample-theoreticql, warning = FALSE, message = FALSE}
# visualize the theoretical null distribution and test statistic!
null_distribution_2_sample_theoretical %>%
  visualize(method = "theoretical") + 
  shade_p_value(observed_statistic,
                direction = "two-sided")
```

We could also visualize both null distributions to get a sense of how the randomization-based and theoretical versions relate. To do so, start off with the randomization-based null distribution, and then drop in `method = "both"` to `visualize()`.

```{r visualize-2-sample-both, warning = FALSE, message = FALSE}
# visualize both null distributions and test statistic!
null_distribution_2_sample_permute %>%
  visualize(method = "both") + 
  shade_p_value(observed_statistic,
                direction = "two-sided")
```

Regardless, it looks like our observed statistic of `r observed_statistic` would be really unlikely if there was truly no relationship between degree status and number of hours worked. More exactly, we can calculate the p-value; theoretical p-values are not yet supported, so we'll use the randomization-based null distribution to do calculate the p-value.

```{r p-value-2-sample, warning = FALSE, message = FALSE}
# calculate the p value from the randomization-based null 
# distribution and the observed statistic
p_value_2_sample <- null_distribution_2_sample_permute %>%
  get_p_value(obs_stat = observed_statistic,
              direction = "two-sided")

p_value_2_sample
```

Thus, if there were really no relationship between the number of hours worked a week and whether one has a college degree, the probability that we would see a statistic as or more extreme than `r observed_statistic` is approximately `r p_value_2_sample`.

Note that, equivalently to the steps shown above, the package supplies a wrapper function, `t_test`, to carry out 2-sample t-Tests on tidy data. The syntax looks like this:

```{r 2-sample-t-test-wrapper, message = FALSE, warning = FALSE}
t_test(x = gss, 
       formula = hours ~ college, 
       order = c("degree", "no degree"),
       alternative = "two-sided")
```

In the above example, we specified the relationship with the syntax `formula = hours ~ college`; we could have also written `response = hours, explanatory = college`.
