---
title: "Full infer Pipeline Examples"
description: "A near-exhaustive demonstration of the functionality in infer."
output: 
  rmarkdown::html_vignette:
    df_print: kable
    toc: true
vignette: |
  %\VignetteIndexEntry{Full infer pipeline examples}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

#### Introduction

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

This vignette is intended to provide a set of examples that nearly exhaustively demonstrate the functionalities provided by `infer`. Commentary on these examples is limited---for more discussion of the intuition behind the package, see the "Getting to Know infer" vignette, accessible by calling `vignette("infer")`.

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 load-packages, echo = FALSE}
library(dplyr)
devtools::load_all()
```


```{r load-gss}
# load in the dataset
data(gss)

# take a look at its structure
dplyr::glimpse(gss)
```

## Hypothesis tests

### One numerical variable (mean)

Calculating the observed statistic,

```{r}
x_bar <- gss %>%
  specify(response = hours) %>%
  calculate(stat = "mean")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  generate(reps = 1000) %>%
  calculate(stat = "mean")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = x_bar, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = x_bar, direction = "two-sided")
```

### One numerical variable (standardized mean $t$)

Calculating the observed statistic,

```{r}
t_bar <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  calculate(stat = "t")
```

Alternatively, using the wrapper to calculate the test statistic, 

```{r}
t_bar <- gss %>%
  t_stat(response = hours, mu = 40)
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  generate(reps = 1000) %>%
  calculate(stat = "t")
```

Alternatively, finding the null distribution using theoretical methods by skipping the `generate()` step,

```{r}
null_distn_theoretical <- gss %>%
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  calculate(stat = "t")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = t_bar, direction = "two-sided")
```

Alternatively, visualizing the observed statistic using the theory-based null distribution,

```{r}
visualize(null_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = t_bar, direction = "two-sided")
```

Alternatively, visualizing the observed statistic using both of the null distributions,

```{r}
visualize(null_distn, method = "both") +
  shade_p_value(obs_stat = t_bar, direction = "two-sided")
```

Note that the above code makes use of the randomization-based null distribution.

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = t_bar, direction = "two-sided")
```

Alternatively, using the `t_test` wrapper:

```{r}
gss %>%
  t_test(response = hours, mu = 40)
```


### One numerical variable (median)

Calculating the observed statistic,

```{r}
x_tilde <- gss %>%
  specify(response = age) %>%
  calculate(stat = "median")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(response = age) %>%
  hypothesize(null = "point", med = 40) %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "median")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = x_tilde, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = x_tilde, direction = "two-sided")
```

### One categorical (one proportion)

Calculating the observed statistic,

```{r}
p_hat <- gss %>%
  specify(response = sex, success = "female") %>%
  calculate(stat = "prop")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(response = sex, success = "female") %>%
  hypothesize(null = "point", p = .5) %>%
  generate(reps = 1000) %>%
  calculate(stat = "prop")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = p_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = p_hat, direction = "two-sided")
```

Note that logical variables will be coerced to factors:

```{r}
null_distn <- gss %>%
  dplyr::mutate(is_female = (sex == "female")) %>%
  specify(response = is_female, success = "TRUE") %>%
  hypothesize(null = "point", p = .5) %>%
  generate(reps = 1000) %>%
  calculate(stat = "prop")
```

### One categorical variable (standardized proportion $z$)

Calculating the observed statistic,

```{r}
p_hat <- gss %>%
  specify(response = sex, success = "female") %>%
  hypothesize(null = "point", p = .5) %>%
  calculate(stat = "z")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(response = sex, success = "female") %>%
  hypothesize(null = "point", p = .5) %>%
  generate(reps = 1000, type = "simulate") %>%
  calculate(stat = "z")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = p_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = p_hat, direction = "two-sided")
```

The package also supplies a wrapper around `prop.test` for tests of a single proportion on tidy data.

```{r prop_test_1_grp}
prop_test(gss,
          college ~ NULL,
          p = .2)
```

### Two categorical (2 level) variables

The `infer` package provides several statistics to work with data of this type. One of them is the statistic for difference in proportions. 

Calculating the observed statistic,

```{r}
d_hat <- gss %>% 
  specify(college ~ sex, success = "no degree") %>%
  calculate(stat = "diff in props", order = c("female", "male"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(college ~ sex, success = "no degree") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "diff in props", order = c("female", "male"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = d_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = d_hat, direction = "two-sided")
```

`infer` also provides functionality to calculate ratios of proportions. The workflow looks similar to that for `diff in props`.

Calculating the observed statistic,

```{r}
r_hat <- gss %>% 
  specify(college ~ sex, success = "no degree") %>%
  calculate(stat = "ratio of props", order = c("female", "male"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(college ~ sex, success = "no degree") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "ratio of props", order = c("female", "male"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = r_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = r_hat, direction = "two-sided")
```

In addition, the package provides functionality to calculate odds ratios. The workflow also looks similar to that for `diff in props`.

Calculating the observed statistic,

```{r}
or_hat <- gss %>% 
  specify(college ~ sex, success = "no degree") %>%
  calculate(stat = "odds ratio", order = c("female", "male"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(college ~ sex, success = "no degree") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "odds ratio", order = c("female", "male"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = or_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = or_hat, direction = "two-sided")
```

### Two categorical (2 level) variables (z)

Finding the standardized observed statistic,

```{r}
z_hat <- gss %>% 
  specify(college ~ sex, success = "no degree") %>%
  calculate(stat = "z", order = c("female", "male"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(college ~ sex, success = "no degree") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "z", order = c("female", "male"))
```

Alternatively, finding the null distribution using theoretical methods by skipping the `generate()` step,

```{r}
null_distn_theoretical <- gss %>%
  specify(college ~ sex, success = "no degree") %>%
  hypothesize(null = "independence") %>%  
  calculate(stat = "z", order = c("female", "male"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = z_hat, direction = "two-sided")
```

Alternatively, visualizing the observed statistic using the theory-based null distribution,

```{r}
visualize(null_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = z_hat, direction = "two-sided")
```

Alternatively, visualizing the observed statistic using both of the null distributions,

```{r}
visualize(null_distn, method = "both") +
  shade_p_value(obs_stat = z_hat, direction = "two-sided")
```

Note that the above code makes use of the randomization-based null distribution.

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = z_hat, direction = "two-sided")
```

Note the similarities in this plot and the previous one.

The package also supplies a wrapper around `prop.test` to allow for tests of equality of proportions on tidy data.

```{r prop_test_2_grp}
prop_test(gss, 
          college ~ sex,  
          order = c("female", "male"))
```

### One categorical (>2 level) - GoF

Calculating the observed statistic,

Note the need to add in the hypothesized values here to compute the observed statistic.

```{r}
Chisq_hat <- gss %>%
  specify(response = finrela) %>%
  hypothesize(null = "point",
              p = c("far below average" = 1/6,
                    "below average" = 1/6,
                    "average" = 1/6,
                    "above average" = 1/6,
                    "far above average" = 1/6,
                    "DK" = 1/6)) %>%
  calculate(stat = "Chisq")
```

Alternatively, using the `chisq_stat` wrapper to calculate the test statistic,

```{r}
Chisq_hat <- gss %>%
  chisq_stat(response = finrela,
             p = c("far below average" = 1/6,
                   "below average" = 1/6,
                   "average" = 1/6,
                   "above average" = 1/6,
                   "far above average" = 1/6,
                   "DK" = 1/6))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(response = finrela) %>%
  hypothesize(null = "point",
              p = c("far below average" = 1/6,
                    "below average" = 1/6,
                    "average" = 1/6,
                    "above average" = 1/6,
                    "far above average" = 1/6,
                    "DK" = 1/6)) %>%
  generate(reps = 1000, type = "simulate") %>%
  calculate(stat = "Chisq")
```

Alternatively, finding the null distribution using theoretical methods by skipping the `generate()` step,

```{r}
null_distn_theoretical <- gss %>%
  specify(response = finrela) %>%
  hypothesize(null = "point",
              p = c("far below average" = 1/6,
                    "below average" = 1/6,
                    "average" = 1/6,
                    "above average" = 1/6,
                    "far above average" = 1/6,
                    "DK" = 1/6)) %>%
  calculate(stat = "Chisq")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Alternatively, visualizing the observed statistic using the theory-based null distribution,

```{r}
visualize(null_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Alternatively, visualizing the observed statistic using both of the null distributions,

```{r}
visualize(null_distn_theoretical, method = "both") +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Note that the above code makes use of the randomization-based null distribution.

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Alternatively, using the `chisq_test` wrapper:

```{r}
chisq_test(gss, 
           response = finrela,
           p = c("far below average" = 1/6,
                 "below average" = 1/6,
                 "average" = 1/6,
                 "above average" = 1/6,
                 "far above average" = 1/6,
                 "DK" = 1/6))
```

### Two categorical (>2 level): Chi-squared test of independence

Calculating the observed statistic,

```{r}
Chisq_hat <- gss %>%
  specify(formula = finrela ~ sex) %>% 
  calculate(stat = "Chisq")
```

Alternatively, using the wrapper to calculate the test statistic,

```{r}
Chisq_hat <- gss %>%
  chisq_stat(formula = finrela ~ sex)
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(finrela ~ sex) %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000, type = "permute") %>% 
  calculate(stat = "Chisq")
```

Alternatively, finding the null distribution using theoretical methods by skipping the `generate()` step,

```{r}
null_distn_theoretical <- gss %>%
  specify(finrela ~ sex) %>%
  hypothesize(null = "independence") %>% 
  calculate(stat = "Chisq")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Alternatively, visualizing the observed statistic using the theory-based null distribution,

```{r}
visualize(null_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Alternatively, visualizing the observed statistic using both of the null distributions,

```{r}
visualize(null_distn, method = "both") +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Note that the above code makes use of the randomization-based null distribution.

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = Chisq_hat, direction = "greater")
```

Alternatively, using the wrapper to carry out the test,

```{r}
gss %>%
  chisq_test(formula = finrela ~ sex)
```

### One numerical variable, one categorical (2 levels) (diff in means)

Calculating the observed statistic,

```{r}
d_hat <- gss %>% 
  specify(age ~ college) %>% 
  calculate(stat = "diff in means", order = c("degree", "no degree"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(age ~ college) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("degree", "no degree"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = d_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = d_hat, direction = "two-sided")
```

### One numerical variable, one categorical (2 levels) (t)

Finding the standardized observed statistic,

```{r}
t_hat <- gss %>% 
  specify(age ~ college) %>% 
  calculate(stat = "t", order = c("degree", "no degree"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(age ~ college) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "t", order = c("degree", "no degree"))
```

Alternatively, finding the null distribution using theoretical methods by skipping the `generate()` step,

```{r}
null_distn_theoretical <- gss %>%
  specify(age ~ college) %>%
  hypothesize(null = "independence") %>%
  calculate(stat = "t", order = c("degree", "no degree"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = t_hat, direction = "two-sided")
```

Alternatively, visualizing the observed statistic using the theory-based null distribution,

```{r}
visualize(null_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = t_hat, direction = "two-sided")
```

Alternatively, visualizing the observed statistic using both of the null distributions,

```{r}
visualize(null_distn, method = "both") +
  shade_p_value(obs_stat = t_hat, direction = "two-sided")
```

Note that the above code makes use of the randomization-based null distribution.

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = t_hat, direction = "two-sided")
```

Note the similarities in this plot and the previous one.

### One numerical variable, one categorical (2 levels) (diff in medians)

Calculating the observed statistic,

```{r}
d_hat <- gss %>% 
  specify(age ~ college) %>% 
  calculate(stat = "diff in medians", order = c("degree", "no degree"))
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
  specify(age ~ college) %>% # alt: response = age, explanatory = season
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in medians", order = c("degree", "no degree"))
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = d_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = d_hat, direction = "two-sided")
```

### One numerical, one categorical (>2 levels) -  ANOVA

Calculating the observed statistic,

```{r}
F_hat <- gss %>% 
  specify(age ~ partyid) %>%
  calculate(stat = "F")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
   specify(age ~ partyid) %>%
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "F")
```

Alternatively, finding the null distribution using theoretical methods by skipping the `generate()` step,

```{r}
null_distn_theoretical <- gss %>%
   specify(age ~ partyid) %>%
   hypothesize(null = "independence") %>%
   calculate(stat = "F")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = F_hat, direction = "greater")
```

Alternatively, visualizing the observed statistic using the theory-based null distribution,

```{r}
visualize(null_distn_theoretical, method = "theoretical") +
  shade_p_value(obs_stat = F_hat, direction = "greater")
```

Alternatively, visualizing the observed statistic using both of the null distributions,

```{r}
visualize(null_distn, mdthod = "both") +
  shade_p_value(obs_stat = F_hat, direction = "greater")
```

Note that the above code makes use of the randomization-based null distribution.

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = F_hat, direction = "greater")
```

### Two numerical vars - SLR 

Calculating the observed statistic,

```{r}
slope_hat <- gss %>% 
  specify(hours ~ age) %>% 
  calculate(stat = "slope")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
   specify(hours ~ age) %>% 
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "slope")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = slope_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = slope_hat, direction = "two-sided")
```

### Two numerical vars - correlation

Calculating the observed statistic,

```{r}
correlation_hat <- gss %>% 
  specify(hours ~ age) %>% 
  calculate(stat = "correlation")
```

Then, generating the null distribution,

```{r}
null_distn <- gss %>%
   specify(hours ~ age) %>% 
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "correlation")
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(null_distn) +
  shade_p_value(obs_stat = correlation_hat, direction = "two-sided")
```

Calculating the p-value from the null distribution and observed statistic,

```{r}
null_distn %>%
  get_p_value(obs_stat = correlation_hat, direction = "two-sided")
```


### Two numerical vars - SLR (t)

Not currently implemented since $t$ could refer to standardized slope or standardized correlation.

## Confidence intervals

### One numerical (one mean)

Finding the observed statistic,

```{r}
x_bar <- gss %>% 
  specify(response = hours) %>%
  calculate(stat = "mean")
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
   specify(response = hours) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "mean")
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- get_ci(boot, type = "se", point_estimate = x_bar)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```

### One numerical (one mean - standardized)

Finding the observed statistic,

```{r}
t_hat <- gss %>% 
  specify(response = hours) %>%
  hypothesize(null = "point", mu = 40) %>%
  calculate(stat = "t")
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
   specify(response = hours) %>%
   hypothesize(null = "point", mu = 40) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "t")
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = t_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```


### One categorical (one proportion)

Finding the observed statistic,

```{r}
p_hat <- gss %>% 
   specify(response = sex, success = "female") %>%
   calculate(stat = "prop")
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
 specify(response = sex, success = "female") %>%
 generate(reps = 1000, type = "bootstrap") %>%
 calculate(stat = "prop")
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = p_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```

### One categorical variable (standardized proportion $z$)

Not yet implemented.

### One numerical variable, one categorical (2 levels) (diff in means)

Finding the observed statistic,

```{r}
d_hat <- gss %>%
  specify(hours ~ college) %>%
  calculate(stat = "diff in means", order = c("degree", "no degree"))
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
   specify(hours ~ college) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "diff in means", order = c("degree", "no degree"))
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = d_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```

### One numerical variable, one categorical (2 levels) (t)

Finding the standardized point estimate,

```{r}
t_hat <- gss %>%
  specify(hours ~ college) %>%
  calculate(stat = "t", order = c("degree", "no degree"))
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
   specify(hours ~ college) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "t", order = c("degree", "no degree"))
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = t_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```


### Two categorical variables (diff in proportions)

Finding the observed statistic,

```{r}
d_hat <- gss %>% 
  specify(college ~ sex, success = "degree") %>%
  calculate(stat = "diff in props", order = c("female", "male"))
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
  specify(college ~ sex, success = "degree") %>%
  generate(reps = 1000, type = "bootstrap") %>% 
  calculate(stat = "diff in props", order = c("female", "male"))
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = d_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```

### Two categorical variables (z)

Finding the standardized point estimate,

```{r}
z_hat <- gss %>% 
  specify(college ~ sex, success = "degree") %>%
  calculate(stat = "z", order = c("female", "male"))
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
  specify(college ~ sex, success = "degree") %>%
  generate(reps = 1000, type = "bootstrap") %>% 
  calculate(stat = "z", order = c("female", "male"))
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = z_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```


### Two numerical vars - SLR

Finding the observed statistic,

```{r}
slope_hat <- gss %>% 
  specify(hours ~ age) %>%
  calculate(stat = "slope")
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
   specify(hours ~ age) %>% 
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "slope")
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = slope_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```

### Two numerical vars - correlation

Finding the observed statistic,

```{r}
correlation_hat <- gss %>% 
  specify(hours ~ age) %>%
  calculate(stat = "correlation")
```

Then, generating the null distribution,

```{r}
boot <- gss %>%
   specify(hours ~ age) %>% 
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "correlation")
```

Use the null distribution to find a confidence interval,

```{r}
percentile_ci <- get_ci(boot)
```

Visualizing the observed statistic alongside the null distribution,

```{r}
visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
```

Alternatively, use the null distribution to find a confidence interval using the standard error,

```{r}
standard_error_ci <- boot %>%
  get_ci(type = "se", point_estimate = correlation_hat)

visualize(boot) +
  shade_confidence_interval(endpoints = standard_error_ci)
```


### Two numerical vars - t

Not currently implemented since $t$ could refer to standardized slope or standardized correlation.
