---
title: "Tidymodels"
subtitle: "Lecture 23"
author: "Dr. Colin Rundel"
footer: "Sta 323 - Spring 2024"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
    self-contained: true
execute: 
  echo: true
  message: true
---

```{r setup}
#| message: False
#| warning: False
#| include: False


knitr::opts_chunk$set(
  fig.align = "center", fig.retina = 2, dpi = 150,
  out.width = "100%"
)

library(tidyverse)
library(patchwork)

ggplot2::theme_set(ggplot2::theme_bw())

options(
  width=90
)

library(rstanarm)
```

#

![](imgs/hex_tidymodels.png){fig-align="center" width="50%"}

## Tidymodels

::: {.small}
```{r message=TRUE}
library(tidymodels)
```
:::

## Book data

```{r books}
#| output-location: column
(books = DAAG::allbacks |>
  as_tibble() |>
  select(-area) |>
  mutate(
    cover = forcats::fct_recode(
      cover, 
      "hardback" = "hb", 
      "paperback" = "pb"
    )
  )
)
```

##

```{r}
ggplot(books, aes(x=volume, y=weight, color = cover)) +
  geom_point(size=2)
```


## Building a tidymodel

```{r}
linear_reg()
```

. . .

```{r}
linear_reg() |>
  set_engine("lm")
```


## Building a tidymodel

```{r}
#| include: false
options(width=60)
```

:::: {.columns .small}
::: {.column width='50%'}
```{r}
linear_reg() |>
  set_engine("lm") |>
  fit(weight ~ volume * cover, data = books)
```
:::

::: {.column width='50%' .fragment}
```{r}
lm(weight ~ volume * cover, data = books)
```
:::
::::

## Tidy model objects

:::: {.columns .small}
::: {.column width='50%'}
```{r}
lm_tm = linear_reg() |>
  set_engine("lm") |>
  fit(weight ~ volume * cover, data = books)
```
:::

::: {.column width='50%'}
```{r}
lm_b = lm(weight ~ volume * cover, data = books)
```
:::
::::

:::: {.columns .small}
::: {.column width='50%'}
```{r}
summary(lm_tm)
```
:::

::: {.column width='50%'}
```{r}
summary(lm_b)
```
:::
::::

#

![](imgs/hex_broom.png){fig-align="center" width="75%"}




## Tidy coefficients

```{r}
#| include: false
options(width=50)
```

:::: {.columns .small}
::: {.column width='50%'}
```{r}
broom::tidy(lm_tm)
```
:::

::: {.column width='50%'}
```{r}
broom::tidy(lm_b)
```
:::
::::

```{r}
#| include: false
options(width=85)
```



## Tidy statistics

::: {.small}
```{r}
broom::glance(lm_tm)
```
:::

. . .

::: {.small}
```{r}
broom::glance(lm_b)
```
:::


## Tidy prediction

```{r}
broom::augment(lm_tm, new_data = books)
```


## Putting it together

::: {.small}
```{r fig.height=3}
#| output-location: fragment
lm_tm |>
  augment(
    new_data = tidyr::expand_grid(
      volume = seq(0, 1500, by=5),
      cover = c("hardback", "paperback") |> as.factor()
    )
  ) |>
  rename(weight = .pred) |>
  ggplot(aes(x = volume, y = weight, color = cover, group = cover)) +
    geom_line() +
    geom_point(data = books)
```
:::


#

![](imgs/hex_parsnip.png){fig-align="center" width="75%"}


## Why do we care?

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{r}
show_engines("linear_reg")
```
:::

::: {.column width='50%' .fragment}
```{r}
(bayes_tm = linear_reg() |> 
  set_engine(
    "stan", 
    prior_intercept = rstanarm::student_t(df = 1), 
    prior = rstanarm::student_t(df = 1),
    seed = 1234
  ) 
)
```
:::
::::


## Fitting with `rstanarm`

::: {.small}
```{r}
(bayes_tm = bayes_tm |>
  fit(weight ~ volume * cover, data = books)
)
```
:::


::: {.aside}
See `?details_linear_reg_stan` for details within `parsnip`
:::

## What was actually run?

::: {.small}
```{r}
linear_reg() |> 
  set_engine(
    "stan", 
    prior_intercept = rstanarm::student_t(df = 1), 
    prior = rstanarm::student_t(df = 1),
    seed = 1234
  ) |>
  translate()
```
:::


## Back to broom

::: {.small}
```{r}
#| error: True
broom::tidy(bayes_tm)
```
:::

. . .


::: {.small}
```{r}
broom.mixed::tidy(bayes_tm)
```

```{r}
broom.mixed::glance(bayes_tm)
```
:::

## Augment

::: {.medium}
```{r}
augment(bayes_tm, new_data=books)
```
:::


## Predictions

::: {.small}
```{r fig.height=3}
bayes_tm |>
  augment(
    new_data = tidyr::expand_grid(
      volume = seq(0, 1500, by=5),
      cover = c("hardback", "paperback") |> as.factor()
    )
  ) |>
  rename(weight = .pred) |>
  ggplot(aes(x = volume, y = weight, color = cover, group = cover)) +
    geom_line() +
    geom_point(data = books)
```
:::


#

![](imgs/hex_yardstick.png){fig-align="center" width="75%"}

## Performance

:::: {.columns .small}
::: {.column width='50%'}
```{r}
lm_tm |>
  augment(new_data = books) |>
  yardstick::rmse(weight, .pred)
```
:::

::: {.column width='50%'}
```{r}
bayes_tm |>
  augment(new_data = books) |>
  yardstick::rmse(weight, .pred)
```
:::
::::

::: {.aside}
More on combining these type of results later, see [`workflows`](https://workflows.tidymodels.org/)
:::


# Cross validation and<br/>Feature engineering

## The Office & IMDB

::: {.medium}
```{r}
#| message: false
(office_ratings = read_csv("data/office_ratings.csv"))
```
:::

::: {.aside}
These data are from [data.world](https://data.world/anujjain7/the-office-imdb-ratings-dataset), by way of [TidyTuesday](https://github.com/rfordatascience/tidytuesday/blob/master/data/2020/2020-03-17/readme.md).
:::


## Rating vs Air Date

```{r}
#| echo: false
ggplot(office_ratings, aes(x = air_date, y = imdb_rating, color = as_factor(season), size = total_votes)) +
  geom_point() +
  scale_size(trans = "log10", range=c(0.5,4)) +
  scale_color_discrete() +
  labs(color = "season")
```

#

![](imgs/hex_rsample.png){fig-align="center" width="75%"}


## Test-train split

```{r}
#| include: false
options(width=50)
options(tibble.width = 50)
```

::: {.small}
```{r}
set.seed(123)
(office_split = initial_split(office_ratings, prop = 0.8))
```
:::

. . .


:::: {.columns .small}
::: {.column width='50%'}
```{r}
(office_train = training(office_split))
```
:::

::: {.column width='50%'}
```{r}
(office_test = testing(office_split))
```
:::
::::


## Feature engineering with dplyr

```{r}
#| echo: false
options(width=100)
options(tibble.width=100)
```

::: {.small}
```{r}
office_train |>
  mutate(
    season = as_factor(season),
    month = lubridate::month(air_date),
    wday = lubridate::wday(air_date),
    top10_votes = as.integer(total_votes > quantile(total_votes, 0.9))
  )
```
:::

. . .

::: {.center}
Anyone see a potential problem with the code above?
:::


#

![](imgs/hex_recipes.png){fig-align="center" width="75%"}

## Better living through recipes

::: {.small}
```{r}
r = recipe(imdb_rating ~ ., data = office_train)
```  
:::

. . .

::: {.small}
```{r}
summary(r)
```
:::



## Recipe roles

::: {.small}
```{r}
r = recipe(
  imdb_rating ~ ., data = office_train
) |> 
  update_role(title, new_role = "ID")
```

```{r}
summary(r)
```
:::


## Adding features (month & day of week)

:::: {.small}
```{r}
r = recipe(
  imdb_rating ~ ., data = office_train
) |> 
  update_role(title, new_role = "ID") |>
  step_date(air_date, features = c("dow", "month"))
```
```{r}
summary(r)
```
::::


## Adding Holidays

::: {.small}
```{r}
r = recipe(
  imdb_rating ~ ., data = office_train
) |> 
  update_role(title, new_role = "ID") |>
  step_date(air_date, features = c("dow", "month")) |>
  step_holiday(
    air_date, 
    holidays = c("USThanksgivingDay", "USChristmasDay", "USNewYearsDay", "USIndependenceDay"), 
    keep_original_cols = FALSE
  )
```
:::

::: {.small}
```{r}
summary(r)
```
:::


## Seasons as factors

::: {.small}
```{r}
r = recipe(
  imdb_rating ~ ., data = office_train
) |> 
  update_role(title, new_role = "ID") |>
  step_date(air_date, features = c("dow", "month")) |>
  step_holiday(
    air_date, 
    holidays = c("USThanksgivingDay", "USChristmasDay", "USNewYearsDay", "USIndependenceDay"), 
    keep_original_cols = FALSE
  ) |>
  step_num2factor(season, levels = as.character(1:9))
```
:::


::: {.small}
```{r}
summary(r)
```
:::


## Dummy coding

::: {.small}
```{r}
r = recipe(
  imdb_rating ~ ., data = office_train
) |> 
  update_role(title, new_role = "ID") |>
  step_date(air_date, features = c("dow", "month")) |>
  step_holiday(
    air_date, 
    holidays = c("USThanksgivingDay", "USChristmasDay", "USNewYearsDay", "USIndependenceDay"), 
    keep_original_cols = FALSE
  ) |>
  step_num2factor(season, levels = as.character(1:9)) |>
  step_dummy(all_nominal_predictors())
```
:::


::: {.small}
```{r}
summary(r)
```
:::

## `top10_votes`

::: {.small}
```{r}
r = recipe(
  imdb_rating ~ ., data = office_train
) |> 
  update_role(title, new_role = "ID") |>
  step_date(air_date, features = c("dow", "month")) |>
  step_holiday(
    air_date, 
    holidays = c("USThanksgivingDay", "USChristmasDay", "USNewYearsDay", "USIndependenceDay"), 
    keep_original_cols = FALSE
  ) |>
  step_num2factor(season, levels = as.character(1:9)) |>
  step_dummy(all_nominal_predictors()) |>
  step_percentile(total_votes) |>
  step_mutate(top10 = as.integer(total_votes >= 0.9)) |>
  step_rm(total_votes)
```
:::


::: {.small}
```{r}
summary(r)
```
:::


## Preparing a recipe

::: {.small}
```{r}
#| message: true
prep(r)
```
:::


## Baking a recipe

::: {.small}
```{r}
prep(r) |>
  bake(new_data = office_train)
```
:::


## Informative features?

::: {.small}
```{r}
prep(r) |>
  bake(new_data = office_train) |>
  map_int(~ length(unique(.x)))
```
:::


## Removing zero variance predictors

::: {.small}
```{r}
r = recipe(
    imdb_rating ~ ., data = office_train
  ) |> 
  update_role(title, new_role = "ID") |>
  step_date(air_date, features = c("dow", "month")) |>
  step_holiday(
    air_date, 
    holidays = c("USThanksgivingDay", "USChristmasDay", "USNewYearsDay", "USIndependenceDay"), 
    keep_original_cols = FALSE
  ) |>
  step_num2factor(season, levels = as.character(1:9)) |>
  step_dummy(all_nominal_predictors()) |>
  step_percentile(total_votes) |>
  step_mutate(top10 = as.integer(total_votes >= 0.9)) |>
  step_rm(total_votes) |>
  step_zv(all_predictors())
```
:::


##

::: {.small}
```{r}
prep(r) |>
  bake(new_data = office_train)
```
:::


#

![](imgs/hex_workflows.png){fig-align="center" width="75%"}

## Really putting it all together

::: {.small}
```{r}
(office_work = workflow() |>
  add_recipe(r) |>
  add_model(
    linear_reg() |>
    set_engine("lm")
  )
)
```
:::

## Workflow fit

::: {.small}
```{r}
(office_fit = office_work |>
  fit(data = office_train))
```
:::


## Performance

:::: {.columns .small}
::: {.column width='50%'}
```{r}
office_fit |>
  augment(office_train) |>
  rmse(imdb_rating, .pred)
```
:::

::: {.column width='50%'}
```{r}
office_fit |>
  augment(office_test) |>
  rmse(imdb_rating, .pred)
```
:::
::::


## k-fold cross validation

```{r echo=FALSE, out.width="75%", fig.align="center"}
knitr::include_graphics("imgs/kfold-cv.png")
```

## Creating folds

```{r}
#| include: false
options(width=55)
```

:::: {.columns .small}
::: {.column width='50%'}
```{r}
set.seed(123)
(folds = vfold_cv(office_train, v=5))
```
:::

::: {.column width='50%' .fragment}
```{r}
(office_fit_folds = office_work |>
  fit_resamples(folds)
)
```
:::
::::


## Fold performance

```{r}
#| include: false
options(width=85)
```

::: {.small}
```{r}
tune::collect_metrics(office_fit_folds)
```
:::

. . .

::: {.small}
```{r}
tune::collect_metrics(office_fit_folds, summarize = FALSE) |>
  filter(.metric == "rmse")
```
:::

::: {.aside}
More on the `tune` package next time
:::