---
title: "Full infer pipeline examples using `nycflights13` `flights` data"
author: "Chester Ismay"
date: "Updated on 2018-06-14"
output: 
  rmarkdown::html_vignette:
    df_print: kable
vignette: |
  %\VignetteIndexEntry{Full infer pipeline examples}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---


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

## Data preparation

```{r message=FALSE, warning=FALSE}
library(nycflights13)
library(dplyr)
library(ggplot2)
library(stringr)
library(infer)
set.seed(2017)
fli_small <- flights %>% 
  na.omit() %>%
  sample_n(size = 500) %>% 
  mutate(season = case_when(
    month %in% c(10:12, 1:3) ~ "winter",
    month %in% c(4:9) ~ "summer"
  )) %>% 
  mutate(day_hour = case_when(
    between(hour, 1, 12) ~ "morning",
    between(hour, 13, 24) ~ "not morning"
  )) %>% 
  select(arr_delay, dep_delay, season, 
         day_hour, origin, carrier)
```

* Two numeric - `arr_delay`, `dep_delay`
* Two categories 
    - `season` (`"winter"`, `"summer"`), 
    - `day_hour` (`"morning"`, `"not morning"`)
* Three categories - `origin` (`"EWR"`, `"JFK"`, `"LGA"`)
* Sixteen categories - `carrier`

***

# Hypothesis tests

### One numerical variable (mean)

**Observed stat**
```{r}
( x_bar <- fli_small %>%
  specify(response = dep_delay) %>%
  calculate(stat = "mean") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = x_bar, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = x_bar, direction = "two_sided")
```

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

**Observed stat**
```{r}
( t_bar <- fli_small %>%
  specify(response = dep_delay) %>%
  calculate(stat = "t") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = t_bar, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = t_bar, direction = "two_sided")
```


### One numerical variable (median)

**Observed stat**

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

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

visualize(null_distn) +
  shade_p_value(obs_stat = x_tilde, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = x_tilde, direction = "two_sided")
```

### One categorical (one proportion)

**Observed stat**

```{r}
( p_hat <- fli_small %>%
  specify(response = day_hour, success = "morning") %>%
  calculate(stat = "prop") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = p_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = p_hat, direction = "two_sided")
```

Logical variables will be coerced to factors:

```{r}
null_distn <- fli_small %>%
  mutate(day_hour_logical = (day_hour == "morning")) %>%
  specify(response = day_hour_logical, success = "TRUE") %>%
  hypothesize(null = "point", p = .5) %>%
  generate(reps = 1000) %>%
  calculate(stat = "prop")
```

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

Not yet implemented.

### Two categorical (2 level) variables

**Observed stat**

```{r}
( d_hat <- fli_small %>% 
  specify(day_hour ~ season, success = "morning") %>%
  calculate(stat = "diff in props", order = c("winter", "summer")) )
```


```{r}
null_distn <- fli_small %>%
  specify(day_hour ~ season, success = "morning") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "diff in props", order = c("winter", "summer"))

visualize(null_distn) +
  shade_p_value(obs_stat = d_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = d_hat, direction = "two_sided")
```

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

**Standardized observed stat**

```{r}
( z_hat <- fli_small %>% 
  specify(day_hour ~ season, success = "morning") %>%
  calculate(stat = "z", order = c("winter", "summer")) )
```

```{r}
null_distn <- fli_small %>%
  specify(day_hour ~ season, success = "morning") %>%
  hypothesize(null = "independence") %>% 
  generate(reps = 1000) %>% 
  calculate(stat = "z", order = c("winter", "summer"))

visualize(null_distn) +
  shade_p_value(obs_stat = z_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = z_hat, direction = "two_sided")
```

Note the similarities in this plot and the previous one.

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

**Observed stat**

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

```{r}
( Chisq_hat <- fli_small %>%
  specify(response = origin) %>%
  hypothesize(null = "point", 
              p = c("EWR" = .33, "JFK" = .33, "LGA" = .34)) %>% 
  calculate(stat = "Chisq") )
```

```{r}
null_distn <- fli_small %>%
  specify(response = origin) %>%
  hypothesize(null = "point", 
              p = c("EWR" = .33, "JFK" = .33, "LGA" = .34)) %>% 
  generate(reps = 1000, type = "simulate") %>% 
  calculate(stat = "Chisq")

visualize(null_distn) +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
null_distn %>%
  get_p_value(obs_stat = Chisq_hat, direction = "greater")
```

### Two categorical (>2 level) variables

**Observed stat**

```{r}
( Chisq_hat <- fli_small %>%
  specify(formula = day_hour ~ origin) %>% 
  calculate(stat = "Chisq") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = Chisq_hat, direction = "greater")
null_distn %>%
  get_p_value(obs_stat = Chisq_hat, direction = "greater")
```

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

**Observed stat**

```{r}
( d_hat <- fli_small %>% 
  specify(dep_delay ~ season) %>% 
  calculate(stat = "diff in means", order = c("summer", "winter")) )
```

```{r}
null_distn <- fli_small %>%
  specify(dep_delay ~ season) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "diff in means", order = c("summer", "winter"))

visualize(null_distn) +
  shade_p_value(obs_stat = d_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = d_hat, direction = "two_sided")
```

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

**Standardized observed stat**

```{r}
( t_hat <- fli_small %>% 
  specify(dep_delay ~ season) %>% 
  calculate(stat = "t", order = c("summer", "winter")) )
```

```{r}
null_distn <- fli_small %>%
  specify(dep_delay ~ season) %>%
  hypothesize(null = "independence") %>%
  generate(reps = 1000, type = "permute") %>%
  calculate(stat = "t", order = c("summer", "winter"))

visualize(null_distn) +
  shade_p_value(obs_stat = t_hat, direction = "two_sided")
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)

**Observed stat**

```{r}
( d_hat <- fli_small %>% 
  specify(dep_delay ~ season) %>% 
  calculate(stat = "diff in medians", order = c("summer", "winter")) )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = d_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = d_hat, direction = "two_sided")
```

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

**Observed stat**

```{r}
( F_hat <- fli_small %>% 
  specify(arr_delay ~ origin) %>%
  calculate(stat = "F") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = F_hat, direction = "greater")
null_distn %>%
  get_p_value(obs_stat = F_hat, direction = "greater")
```

### Two numerical vars - SLR 

**Observed stat**

```{r}
( slope_hat <- fli_small %>% 
  specify(arr_delay ~ dep_delay) %>% 
  calculate(stat = "slope") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = slope_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = slope_hat, direction = "two_sided")
```

### Two numerical vars - correlation

**Observed stat**

```{r}
( correlation_hat <- fli_small %>% 
  specify(arr_delay ~ dep_delay) %>% 
  calculate(stat = "correlation") )
```

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

visualize(null_distn) +
  shade_p_value(obs_stat = correlation_hat, direction = "two_sided")
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.


```{r echo=FALSE, eval=FALSE}
# **Standardized observed stat**
( t_hat <- fli_small %>% 
  specify(arr_delay ~ dep_delay) %>% 
  calculate(stat = "t") )
```

```{r echo=FALSE, eval=FALSE}
null_distn <- fli_small %>%
   specify(arr_delay ~ dep_delay) %>% 
   hypothesize(null = "independence") %>%
   generate(reps = 1000, type = "permute") %>%
   calculate(stat = "t")

visualize(null_distn) +
  shade_p_value(obs_stat = t_hat, direction = "two_sided")
null_distn %>%
  get_p_value(obs_stat = t_hat, direction = "two_sided")
```


## Confidence intervals

### One numerical (one mean)

**Point estimate**

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

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

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( 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)

**Point estimate**

```{r}
( t_hat <- fli_small %>% 
  specify(response = arr_delay) %>%
  calculate(stat = "t") )
```

```{r}
boot <- fli_small %>%
   specify(response = arr_delay) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "t")
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = t_hat) )

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


### One categorical (one proportion)

**Point estimate**

```{r}
( p_hat <- fli_small %>% 
   specify(response = day_hour, success = "morning") %>%
   calculate(stat = "prop") )
```

```{r}
boot <- fli_small %>%
 specify(response = day_hour, success = "morning") %>%
 generate(reps = 1000, type = "bootstrap") %>%
 calculate(stat = "prop")
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, 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)

**Point estimate**

```{r}
( d_hat <- fli_small %>%
  specify(arr_delay ~ season) %>%
  calculate(stat = "diff in means", order = c("summer", "winter")) )
```

```{r}
boot <- fli_small %>%
   specify(arr_delay ~ season) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "diff in means", order = c("summer", "winter"))
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = d_hat) )

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

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

**Standardized point estimate**

```{r}
( t_hat <- fli_small %>%
  specify(arr_delay ~ season) %>%
  calculate(stat = "t", order = c("summer", "winter")) )
```

```{r}
boot <- fli_small %>%
   specify(arr_delay ~ season) %>%
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "t", order = c("summer", "winter"))
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = t_hat) )

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


### Two categorical variables (diff in proportions)

**Point estimate**

```{r}
( d_hat <- fli_small %>% 
  specify(day_hour ~ season, success = "morning") %>%
  calculate(stat = "diff in props", order = c("summer", "winter")) )
```

```{r}
boot <- fli_small %>%
  specify(day_hour ~ season, success = "morning") %>%
  generate(reps = 1000, type = "bootstrap") %>% 
  calculate(stat = "diff in props", order = c("summer", "winter"))
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = d_hat) )

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

### Two categorical variables (z)

**Standardized point estimate**

```{r}
( z_hat <- fli_small %>% 
  specify(day_hour ~ season, success = "morning") %>%
  calculate(stat = "z", order = c("summer", "winter")) )
```

```{r}
boot <- fli_small %>%
  specify(day_hour ~ season, success = "morning") %>%
  generate(reps = 1000, type = "bootstrap") %>% 
  calculate(stat = "z", order = c("summer", "winter"))
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = z_hat) )

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


### Two numerical vars - SLR

**Point estimate**

```{r}
( slope_hat <- fli_small %>% 
  specify(arr_delay ~ dep_delay) %>%
  calculate(stat = "slope") )
```

```{r}
boot <- fli_small %>%
   specify(arr_delay ~ dep_delay) %>% 
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "slope")
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = slope_hat) )

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

### Two numerical vars - correlation

**Point estimate**

```{r}
( correlation_hat <- fli_small %>% 
  specify(arr_delay ~ dep_delay) %>%
  calculate(stat = "correlation") )
```

```{r}
boot <- fli_small %>%
   specify(arr_delay ~ dep_delay) %>% 
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "correlation")
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, 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.


```{r eval=FALSE, echo=FALSE}
# **Point estimate**
( t_hat <- fli_small %>% 
  specify(arr_delay ~ dep_delay) %>%
  calculate(stat = "t") )
```

```{r eval=FALSE, echo=FALSE}
boot <- fli_small %>%
   specify(arr_delay ~ dep_delay) %>% 
   generate(reps = 1000, type = "bootstrap") %>%
   calculate(stat = "t")
( percentile_ci <- get_ci(boot) )

visualize(boot) +
  shade_confidence_interval(endpoints = percentile_ci)
( standard_error_ci <- get_ci(boot, type = "se", point_estimate = t_hat) )

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