07_high_dimensional_data.Rmd

# High-dimensional data {#machine_learning}

```{r include = FALSE}
# Caching this markdown file
#knitr::opts_chunk$set(cache = TRUE)
```

## The Big Picture

- The rise of high-dimensional data. The new data frontiers in social sciences---text ([Gentzkow et al. 2019](https://web.stanford.edu/~gentzkow/research/text-as-data.pdf); [Grimmer and Stewart 2013](https://www.jstor.org/stable/pdf/24572662.pdf?casa_token=SQdSI4R_VdwAAAAA:4QiVLhCXqr9f0qNMM9U75EL5JbDxxnXxUxyIfDf0U8ZzQx9szc0xVqaU6DXG4nHyZiNkvcwGlgD6H0Lxj3y0ULHwgkf1MZt8-9TPVtkEH9I4AHgbTg)) and and image ([Joo and Steinert-Threlkeld 2018](https://arxiv.org/pdf/1810.01544))---are all high-dimensional data. 

    - 1000 common English words for 30-word tweets: $1000^{30}$ similar to N of atoms in the universe ([Gentzkow et al. 2019](https://web.stanford.edu/~gentzkow/research/text-as-data.pdf))

    - Belloni, Alexandre, Victor Chernozhukov, and Christian Hansen. ["High-dimensional methods and inference on structural and treatment effects."](https://pubs.aeaweb.org/doi/pdfplus/10.1257/jep.28.2.29) *Journal of Economic Perspectives 28*, no. 2 (2014): 29-50.

- The rise of the new approach: statistics + computer science = machine learning 

- Statistical inference 

    - $y$ <- some probability models (e.g., linear regression, logistic regression) <- $x$
       
    - $y$ = $X\beta$ + $\epsilon$
        
    - The goal is to estimate $\beta$

- Machine learning 

    - $y$ <- unknown <- $x$ 
    
    - $y$ <-> decision trees, neutral nets <-> $x$
        
    - For the main idea behind prediction modeling, see Breiman, Leo (Berkeley stat faculty who passed away in 2005). ["Statistical modeling: The two cultures (with comments and a rejoinder by the author)."](https://projecteuclid.org/euclid.ss/1009213726) *Statistical science* 16, no. 3 (2001): 199-231.
    
    - "The problem is to find an algorithm $f(x)$ such that for future $x$ in a test set, $f(x)$ will be a good predictor of $y$."
    
    - "There are **two cultures** in the use of statistical modeling to reach conclusions from data. One assumes that the data are generated by a **given** **stochastic data model**. The other uses **algorithmic models** and treats the data mechanism as **unknown**."


- How ML differs from econometrics? 

- A review by Athey, Susan, and Guido W. Imbens. ["Machine learning methods that economists should know about."](https://www.annualreviews.org/doi/full/10.1146/annurev-economics-080217-053433) *Annual Review of Economics* 11 (2019): 685-725.
        
- Stat:
  
    - Specifying a target (i.e., an estimand)

    - Fitting a model to data using an objective function (e.g., the sum of squared errors)

    - Reporting point estimates (effect size) and standard errors (uncertainty)

    - Validation by yes-no using goodness-of-fit tests and residual examination

- ML: 

    - Developing algorithms (estimating *f(x)*)

    - Prediction power, not structural/causal parameters

    - Basically, high-dimensional data statistics (N < P)

    - The major problem is to avoid ["the curse of dimensionality"](https://en.wikipedia.org/wiki/Curse_of_dimensionality) ([too many features - > overfitting](https://towardsdatascience.com/the-curse-of-dimensionality-50dc6e49aa1e))

    - Validation: out-of-sample comparisons (cross-validation) not in-sample goodness-of-fit measures

    - So, it's curve-fitting, but the primary focus is unseen (test data), not seen data (training data)

- A quick review on ML lingos for those trained in econometrics 

    - Sample to estimate parameters = Training sample
  
    - Estimating the model = Being trained 
    
    - Regressors, covariates, or predictors = Features 
  
    - Regression parameters = weights 
  
    - Prediction problems = Supervised (some $y$ are known) + Unsupervised ($y$ unknown)
    
![How to teach machines. Based on [vas3k blog](https://vas3k.com/blog/machine_learning/). Many images in this chapter come from vas3k blog.](https://i.vas3k.ru/7w9.jpg)

![The main types of machine learning. Based on [vas3k blog](https://vas3k.com/blog/machine_learning/)](https://i.vas3k.ru/7vz.jpg)

![The map of the machine learning universe. Based on [vas3k blog](https://vas3k.com/blog/machine_learning/)](https://i.vas3k.ru/7vx.jpg)

![Classical machine learning. Based on [vas3k blog](https://vas3k.com/blog/machine_learning/)](https://i.vas3k.ru/7w1.jpg)

## Dataset  

- [Heart disease data from UCI](https://archive.ics.uci.edu/ml/datasets/heart+Disease)

- One of the popular datasets used in machine learning competitions 

```{r}
# Load packages

## CRAN packages
pacman::p_load(
  here,
  tidyverse,
  tidymodels,
  doParallel, # parallel processing
  patchwork, # arranging ggplots
  remotes,
  SuperLearner,
  vip,
  tidymodels,
  glmnet,
  xgboost,
  rpart,
  ranger,
  conflicted
)

remotes::install_github("ck37/ck37r")

conflicted::conflict_prefer("filter", "dplyr")
```


```{r}

## Jae's custom functions
source(here("functions", "ml_utils.r"))

# Import the dataset

data_original <- read_csv(here("data", "heart.csv"))

glimpse(data_original)

# Createa a copy
data <- data_original

theme_set(theme_minimal())
```

## Workflow 

```{r echo=FALSE, results='asis'}

pacman::p_load(glue)

workflow_list <- c(
  "Preprocessing",
  "Model building",
  "Model fitting",
  "Model evaluation",
  "Model tuning",
  "Prediction"
)

glue("- {seq(workflow_list)}. {workflow_list}")
```

## tidymodels 

- Like `tidyverse`, `tidymodels` is a collection of packages.

    - [`rsample`](https://rsample.tidymodels.org/): for data splitting 
    
    - [`recipes`](https://recipes.tidymodels.org/index.html): for pre-processing
    
    - [`parsnip`](https://www.tidyverse.org/blog/2018/11/parsnip-0-0-1/): for model building 
    
        - [`tune`](https://github.com/tidymodels/tune): hyperparameter tuning 
    
    - [`yardstick`](https://github.com/tidymodels/yardstick): for model evaluations 
    
    - [`workflows`](https://github.com/tidymodels/workflows): for bundling a pieplne that bundles together preprocessing, modeling, and post-processing requests 
    
- Why taking a tidyverse approach to machine learning?

- Benefits 

    - Readable code 
    
    - Reusable data structures 
    
    - Extendable code

![Tidymodels. From RStudio.](https://rviews.rstudio.com/post/2019-06-14-a-gentle-intro-to-tidymodels_files/figure-html/ds.png)

> tidymodels are an **integrated, modular, extensible** set of packages that implement a framework that facilitates creating predicative stochastic models. - Joseph Rickert@RStudio

- Currently, 238 models are [available](https://topepo.github.io/caret/available-models.html) 

- The following materials are based on [the machine learning with tidymodels workshop](https://github.com/dlab-berkeley/Machine-Learning-with-tidymodels) I developed for D-Lab. [The original workshop](https://github.com/dlab-berkeley/Machine-Learning-in-R) was designed by [Chris Kennedy](https://ck37.com/) and [Evan Muzzall](https://dlab.berkeley.edu/people/evan-muzzall.

## Pre-processing

- [`recipes`](https://recipes.tidymodels.org/index.html): for pre-processing

- [`textrecipes`](https://github.com/tidymodels/textrecipes) for text pre-processing

- Step 1: `recipe()` defines target and predictor variables (ingredients).

- Step 2: `step_*()` defines preprocessing steps to be taken (recipe).

    The preprocessing steps list draws on the vignette of the [`parsnip`](https://www.tidymodels.org/find/parsnip/) package.

    - dummy: Also called one-hot encoding

    - zero variance: Removing columns (or features) with a single unique value  

    - impute: Imputing missing values

    - decorrelate: Mitigating correlated predictors (e.g., principal component analysis)

    - normalize: Centering and/or scaling predictors (e.g., log scaling). Scaling matters because many algorithms (e.g., lasso) are scale-variant (except tree-based algorithms). Remind you that normalization (sensitive to outliers) = $\frac{X - X_{min}}{X_{max} - X_{min}}$ and standardization (not sensitive to outliers) = $\frac{X - \mu}{\sigma}$

    - transform: Making predictors symmetric 

- Step 3: `prep()` prepares a dataset to base each step on.

- Step 4: `bake()` applies the preprocessing steps to your datasets. 

In this course, we focus on two preprocessing tasks. 

- One-hot encoding (creating dummy/indicator variables)

```{r}

# Turn selected numeric variables into factor variables
data <- data %>%
  dplyr::mutate(across(c("sex", "ca", "cp", "slope", "thal"), as.factor))

glimpse(data)
```
- Imputation 

```{r}

# Check missing values

map_df(data, ~ is.na(.) %>% sum())

# Add missing values

data$oldpeak[sample(seq(data), size = 10)] <- NA

# Check missing values

# Check the number of missing values
data %>%
  map_df(~ is.na(.) %>% sum())

# Check the rate of missing values
data %>%
  map_df(~ is.na(.) %>% mean())
```
 
### Regression setup 

#### Outcome variable 

```{r}

# Continuous variable
data$age %>% class()
```

#### Data splitting using random sampling 

```{r}

# for reproducibility
set.seed(1234)

# split
split_reg <- initial_split(data, prop = 0.7)

# training set
raw_train_x_reg <- training(split_reg)

# test set
raw_test_x_reg <- testing(split_reg)
```

#### recipe 

```{r}

# Regression recipe
rec_reg <- raw_train_x_reg %>%
  # Define the outcome variable
  recipe(age ~ .) %>%
  # Median impute oldpeak column
  step_impute_median(oldpeak) %>%
  # Expand "sex", "ca", "cp", "slope", and "thal" features out into dummy variables (indicators).
  step_dummy(c("sex", "ca", "cp", "slope", "thal"))

# Prepare a dataset to base each step on
prep_reg <- rec_reg %>% prep(retain = TRUE)
```

```{r}

# x features
train_x_reg <- juice(prep_reg, all_predictors())

test_x_reg <- bake(
  object = prep_reg,
  new_data = raw_test_x_reg, all_predictors()
)

# y variables
train_y_reg <- juice(prep_reg, all_outcomes())$age %>% as.numeric()
test_y_reg <- bake(prep_reg, raw_test_x_reg, all_outcomes())$age %>% as.numeric()

# Checks
names(train_x_reg) # Make sure there's no age variable!

class(train_y_reg) # Make sure this is a continuous variable!
```
- Note that other imputation methods are also available. 

```{r}
grep("impute", ls("package:recipes"), value = TRUE)
```

- You can also create your own `step_` functions. For more information, see [tidymodels.org](https://www.tidymodels.org/learn/develop/recipes/).

### Classification setup 

#### Outcome variable 

```{r}

data$target %>% class()

data$target <- as.factor(data$target)

data$target %>% class()
```

#### Data splitting using stratified random sampling

```{r}

# split
split_class <- initial_split(data %>%
  mutate(target = as.factor(target)),
prop = 0.7,
strata = target
)

# training set
raw_train_x_class <- training(split_class)

# testing set
raw_test_x_class <- testing(split_class)
```

#### recipe 

```{r}

# Classification recipe
rec_class <- raw_train_x_class %>%
  # Define the outcome variable
  recipe(target ~ .) %>%
  # Median impute oldpeak column
  step_impute_median(oldpeak) %>%
  # Expand "sex", "ca", "cp", "slope", and "thal" features out into dummy variables (indicators).
  step_normalize(age) %>%
  step_dummy(c("sex", "ca", "cp", "slope", "thal"))

# Prepare a dataset to base each step on
prep_class <- rec_class %>% prep(retain = TRUE)
```

```{r}

# x features
train_x_class <- juice(prep_class, all_predictors())
test_x_class <- bake(prep_class, raw_test_x_class, all_predictors())

# y variables
train_y_class <- juice(prep_class, all_outcomes())$target %>% as.factor()
test_y_class <- bake(prep_class, raw_test_x_class, all_outcomes())$target %>% as.factor()

# Checks
names(train_x_class) # Make sure there's no target variable!
class(train_y_class) # Make sure this is a factor variable!
```

## Supervised learning

x -> f - > y (defined)

### OLS and Lasso

#### parsnip 

- Build models (`parsnip`)

1. Specify a model 
2. Specify an engine 
3. Specify a mode 

```{r}
# OLS spec
ols_spec <- linear_reg() %>% # Specify a model
  set_engine("lm") %>% # Specify an engine: lm, glmnet, stan, keras, spark
  set_mode("regression") # Declare a mode: regression or classification
```

![Source: http://ethen8181.github.io](http://ethen8181.github.io/machine-learning/regularization/images/lasso_ridge_coefficients.png)

Lasso is one of the regularization techniques along with ridge and elastic-net. 

```{r}
# Lasso spec
lasso_spec <- linear_reg(
  penalty = 0.1, # tuning hyperparameter
  mixture = 1
) %>% # 1 = lasso, 0 = ridge
  set_engine("glmnet") %>%
  set_mode("regression")

# If you don't understand parsnip arguments
lasso_spec %>% translate() # See the documentation
```

- Fit models 

```{r}

ols_fit <- ols_spec %>%
  fit_xy(x = train_x_reg, y = train_y_reg)
# fit(train_y_reg ~ ., train_x_reg) # When you data are not preprocessed

lasso_fit <- lasso_spec %>%
  fit_xy(x = train_x_reg, y = train_y_reg)
```

#### yardstick 

- Visualize model fits 

```{r}

map2(list(ols_fit, lasso_fit), c("OLS", "Lasso"), visualize_fit)
```

```{r}

# Define performance metrics
metrics <- yardstick::metric_set(rmse, mae, rsq)

# Evaluate many models
evals <- purrr::map(list(ols_fit, lasso_fit), evaluate_reg) %>%
  reduce(bind_rows) %>%
  mutate(type = rep(c("OLS", "Lasso"), each = 3))

# Visualize the test results
evals %>%
  ggplot(aes(x = fct_reorder(type, .estimate), y = .estimate)) +
  geom_point() +
  labs(
    x = "Model",
    y = "Estimate"
  ) +
  facet_wrap(~ glue("{toupper(.metric)}"), scales = "free_y")
```
- For more information, read [Tidy Modeling with R](https://www.tmwr.org/) by Max Kuhn and Julia Silge.

#### tune 

**Hyper**parameters are parameters that control the learning process.

##### tune ingredients 

-   Search space for hyperparameters

1. Grid search: a grid of hyperparameters 

2. Random search: random sample points from a bounded domain

![](https://www.programmersought.com/images/523/7e44435f20fe514c11ca0d930af8547b.png)

```{r}
# tune() = placeholder

tune_spec <- linear_reg(
  penalty = tune(), # tuning hyperparameter
  mixture = 1
) %>% # 1 = lasso, 0 = ridge
  set_engine("glmnet") %>%
  set_mode("regression")

tune_spec

# penalty() searches 50 possible combinations

lambda_grid <- grid_regular(penalty(), levels = 50)
```

![Source: Kaggle](https://www.googleapis.com/download/storage/v1/b/kaggle-forum-message-attachments/o/inbox%2F4788946%2F82b5a41b6693a313b246f02d79e972d5%2FK%20FOLD.png?generation=1608195745131795&alt=media)

```{r}
# 10-fold cross-validation

set.seed(1234) # for reproducibility

rec_folds <- vfold_cv(train_x_reg %>% bind_cols(tibble(age = train_y_reg)))
```

##### Add these elements to a workflow 

```{r}

# Workflow
rec_wf <- workflow() %>%
  add_model(tune_spec) %>%
  add_formula(age ~ .)
```

```{r}

# Tuning results
rec_res <- rec_wf %>%
  tune_grid(
    resamples = rec_folds,
    grid = lambda_grid
  )
```

##### Visualize 

```{r}

# Visualize

rec_res %>%
  collect_metrics() %>%
  ggplot(aes(penalty, mean, col = .metric)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err
  ),
  alpha = 0.3
  ) +
  geom_line(size = 2) +
  scale_x_log10() +
  labs(x = "log(lambda)") +
  facet_wrap(~ glue("{toupper(.metric)}"),
    scales = "free",
    nrow = 2
  ) +
  theme(legend.position = "none")
```

##### Select 

```{r}
conflict_prefer("filter", "dplyr")

top_rmse <- show_best(rec_res, metric = "rmse")

best_rmse <- select_best(rec_res, metric = "rmse")

best_rmse

glue('The RMSE of the intiail model is
     {evals %>%
  filter(type == "Lasso", .metric == "rmse") %>%
  select(.estimate) %>%
  round(2)}')

glue('The RMSE of the tuned model is {rec_res %>%
  collect_metrics() %>%
  filter(.metric == "rmse") %>%
  arrange(mean) %>%
  dplyr::slice(1) %>%
  select(mean) %>%
  round(2)}')
```

- Finalize your workflow and visualize [variable importance](https://koalaverse.github.io/vip/articles/vip.html)

```{r}

finalize_lasso <- rec_wf %>%
  finalize_workflow(best_rmse)

finalize_lasso %>%
  fit(train_x_reg %>% bind_cols(tibble(age = train_y_reg))) %>%
  pull_workflow_fit() %>%
  vip::vip()
```

##### Test fit 

- Apply the tuned model to the test dataset 

```{r}

test_fit <- finalize_lasso %>%
  fit(test_x_reg %>% bind_cols(tibble(age = test_y_reg)))

evaluate_reg(test_fit)
```

### Decision tree 

#### parsnip 

- Build a model 

1. Specify a model 
2. Specify an engine 
3. Specify a mode 

```{r}

# workflow
tree_wf <- workflow() %>% add_formula(target ~ .)

# spec
tree_spec <- decision_tree(

  # Mode
  mode = "classification",

  # Tuning hyperparameters
  cost_complexity = NULL,
  tree_depth = NULL
) %>%
  set_engine("rpart") # rpart, c5.0, spark

tree_wf <- tree_wf %>% add_model(tree_spec)
```

- Fit a model

```{r}

tree_fit <- tree_wf %>% fit(train_x_class %>% bind_cols(tibble(target = train_y_class)))
```

#### yardstick 

- Let's formally test prediction performance. 

1.  Confusion matrix

A confusion matrix is often used to describe the performance of a classification model. The below example is based on a binary classification model.

|                 | Predicted: YES      | Predicted: NO       |
|-----------------|---------------------|---------------------|
| **Actual: YES** | True positive (TP)  | False negative (FN) |
| **Actual: NO**  | False positive (FP) | True negative (TN)  |

2.  Metrics

-   `accuracy`: The proportion of the data predicted correctly ($\frac{TP + TN}{total}$). 1 - accuracy = misclassification rate.

-   `precision`: Positive predictive value. *When the model predicts yes, how correct is it?* ($\frac{TP}{TP + FP}$)

-   `recall` (sensitivity): True positive rate (e.g., healthy people healthy). *When the actual value is yes, how often does the model predict yes?* ($\frac{TP}{TP + FN}$)

-   `F-score`: A weighted average between precision and recall. 

-   `ROC Curve` (receiver operating characteristic curve): a plot that shows the relationship between true and false positive rates at different classification thresholds. y-axis indicates the true positive rate and x-axis indicates the false positive rate. What matters is the AUC (Area under the ROC Curve), which is a cumulative probability function of ranking a random "positive" - "negative" pair (for the probability of AUC, see [this blog post](https://www.alexejgossmann.com/auc/)).

![Source: Google Machine Learning Crash Course](https://developers.google.com/machine-learning/crash-course/images/ROCCurve.svg)

-   To learn more about other metrics, check out the yardstick package [references](https://yardstick.tidymodels.org/reference/index.html).

```{r}

# Define performance metrics

metrics <- yardstick::metric_set(accuracy, precision, recall)

# Visualize

tree_fit_viz_metr <- visualize_class_eval(tree_fit)

tree_fit_viz_metr

tree_fit_viz_mat <- visualize_class_conf(tree_fit)

tree_fit_viz_mat
```

#### tune 

##### tune ingredients 

Decisions trees tend to overfit. We need to consider two things to reduce this problem: how to split and when to stop a tree.

- **complexity parameter**: a high CP means a simple decision tree with few splits. 

- **tree_depth** 

```{r}

tune_spec <- decision_tree(
  cost_complexity = tune(), # how to split
  tree_depth = tune(), # when to stop
  mode = "classification"
) %>%
  set_engine("rpart")

tree_grid <- grid_regular(cost_complexity(),
  tree_depth(),
  levels = 5
) # 2 hyperparameters -> 5*5 = 25 combinations

tree_grid %>%
  count(tree_depth)

# 10-fold cross-validation

set.seed(1234) # for reproducibility

tree_folds <- vfold_cv(train_x_class %>% bind_cols(tibble(target = train_y_class)),
  strata = target
)
```

##### Add these elements to a workflow 

```{r}

# Update workflow
tree_wf <- tree_wf %>% update_model(tune_spec)

# Determine the number of cores
no_cores <- detectCores() - 1

# Initiate
cl <- makeCluster(no_cores)

registerDoParallel(cl)

# Tuning results
tree_res <- tree_wf %>%
  tune_grid(
    resamples = tree_folds,
    grid = tree_grid,
    metrics = metrics
  )
```

##### Visualize 

- The following plot draws on the [vignette](https://www.tidymodels.org/start/tuning/) of the tidymodels package. 

```{r}

tree_res %>%
  collect_metrics() %>%
  mutate(tree_depth = factor(tree_depth)) %>%
  ggplot(aes(cost_complexity, mean, col = .metric)) +
  geom_point(size = 3) +
  # Subplots
  facet_wrap(~tree_depth,
    scales = "free",
    nrow = 2
  ) +
  # Log scale x
  scale_x_log10(labels = scales::label_number()) +
  # Discrete color scale
  scale_color_viridis_d(option = "plasma", begin = .9, end = 0) +
  labs(
    x = "Cost complexity",
    col = "Tree depth",
    y = NULL
  ) +
  coord_flip()
```

##### Select 

```{r}

# Optimal hyperparameter
best_tree <- select_best(tree_res, "recall")

# Add the hyperparameter to the workflow
finalize_tree <- tree_wf %>%
  finalize_workflow(best_tree)
```

```{r}

tree_fit_tuned <- finalize_tree %>%
  fit(train_x_class %>% bind_cols(tibble(target = train_y_class)))

# Metrics
(tree_fit_viz_metr + labs(title = "Non-tuned")) / (visualize_class_eval(tree_fit_tuned) + labs(title = "Tuned"))

# Confusion matrix
(tree_fit_viz_mat + labs(title = "Non-tuned")) / (visualize_class_conf(tree_fit_tuned) + labs(title = "Tuned"))
```

- Visualize variable importance 

```{r}

tree_fit_tuned %>%
  pull_workflow_fit() %>%
  vip::vip()
```

##### Test fit

- Apply the tuned model to the test dataset 

```{r}

test_fit <- finalize_tree %>%
  fit(test_x_class %>% bind_cols(tibble(target = test_y_class)))

evaluate_class(test_fit)
```

In the next subsection, we will learn variants of ensemble models that improve decision tree models by putting models together.

### Bagging (Random forest)

Key idea applied across all ensemble models (bagging, boosting, and stacking): 
single learner -> N learners (N > 1) 

Many learners could perform better than a single learner as this approach reduces the **variance** of a single estimate and provides more stability.

Here we focus on the difference between bagging and boosting. In short, boosting may reduce bias while increasing variance. On the other hand, bagging may reduce variance but has nothing to do with bias. Please check out [What is the difference between Bagging and Boosting?](https://quantdare.com/what-is-the-difference-between-bagging-and-boosting/) by aporras.

**bagging**

- Data: Training data will be randomly sampled with replacement (bootstrapping samples + drawing random **subsets** of features for training individual trees)

- Learning: Building models in parallel (independently)

- Prediction: Simple average of the estimated responses (majority vote system)


![From Sebastian Raschka's blog](https://sebastianraschka.com/images/faq/bagging-boosting-rf/bagging.png)


**boosting** 


- Data: Weighted training data will be random sampled

- Learning: Building models sequentially (mispredicted cases would receive more weights) 

- Prediction: Weighted average of the estimated responses 


![From Sebastian Raschka's blog](https://sebastianraschka.com/images/faq/bagging-boosting-rf/boosting.png)


#### parsnip 

- Build a model 

1. Specify a model 
2. Specify an engine 
3. Specify a mode 

```{r}

# workflow
rand_wf <- workflow() %>% add_formula(target ~ .)

# spec
rand_spec <- rand_forest(

  # Mode
  mode = "classification",

  # Tuning hyperparameters
  mtry = NULL, # The number of predictors to available for splitting at each node
  min_n = NULL, # The minimum number of data points needed to keep splitting nodes
  trees = 500
) %>% # The number of trees
  set_engine("ranger",
    # We want the importance of predictors to be assessed.
    seed = 1234,
    importance = "permutation"
  )

rand_wf <- rand_wf %>% add_model(rand_spec)
```

- Fit a model

```{r}

rand_fit <- rand_wf %>% fit(train_x_class %>% bind_cols(tibble(target = train_y_class)))
```

#### yardstick 

```{r}

# Define performance metrics
metrics <- yardstick::metric_set(accuracy, precision, recall)

rand_fit_viz_metr <- visualize_class_eval(rand_fit)

rand_fit_viz_metr
```

- Visualize the confusion matrix. 
  
```{r}

rand_fit_viz_mat <- visualize_class_conf(rand_fit)

rand_fit_viz_mat
```

#### tune 

##### tune ingredients 

We focus on the following two hyperparameters:

- `mtry`: The number of predictors available for splitting at each node.

- `min_n`: The minimum number of data points needed to keep splitting nodes. 

```{r}

tune_spec <-
  rand_forest(
    mode = "classification",

    # Tuning hyperparameters
    mtry = tune(),
    min_n = tune()
  ) %>%
  set_engine("ranger",
    seed = 1234,
    importance = "permutation"
  )

rand_grid <- grid_regular(mtry(range = c(1, 10)),
  min_n(range = c(2, 10)),
  levels = 5
)

rand_grid %>%
  count(min_n)
```

```{r}
# 10-fold cross-validation

set.seed(1234) # for reproducibility

rand_folds <- vfold_cv(train_x_class %>% bind_cols(tibble(target = train_y_class)),
  strata = target
)
```

##### Add these elements to a workflow 

```{r}

# Update workflow
rand_wf <- rand_wf %>% update_model(tune_spec)

# Tuning results
rand_res <- rand_wf %>%
  tune_grid(
    resamples = rand_folds,
    grid = rand_grid,
    metrics = metrics
  )
```

##### Visualize 

```{r}
rand_res %>%
  collect_metrics() %>%
  mutate(min_n = factor(min_n)) %>%
  ggplot(aes(mtry, mean, color = min_n)) +
  # Line + Point plot
  geom_line(size = 1.5, alpha = 0.6) +
  geom_point(size = 2) +
  # Subplots
  facet_wrap(~.metric,
    scales = "free",
    nrow = 2
  ) +
  # Log scale x
  scale_x_log10(labels = scales::label_number()) +
  # Discrete color scale
  scale_color_viridis_d(option = "plasma", begin = .9, end = 0) +
  labs(
    x = "The number of predictors to be sampled",
    col = "The minimum number of data points needed for splitting",
    y = NULL
  ) +
  theme(legend.position = "bottom")
```

```{r}
# Optimal hyperparameter
best_tree <- select_best(rand_res, "accuracy")

best_tree

# Add the hyperparameter to the workflow
finalize_tree <- rand_wf %>%
  finalize_workflow(best_tree)
```

```{r}
rand_fit_tuned <- finalize_tree %>%
  fit(train_x_class %>% bind_cols(tibble(target = train_y_class)))

# Metrics
(rand_fit_viz_metr + labs(title = "Non-tuned")) / (visualize_class_eval(rand_fit_tuned) + labs(title = "Tuned"))

# Confusion matrix
(rand_fit_viz_mat + labs(title = "Non-tuned")) / (visualize_class_conf(rand_fit_tuned) + labs(title = "Tuned"))
```

- Visualize variable importance 

```{r}
rand_fit_tuned %>%
  pull_workflow_fit() %>%
  vip::vip()
```

##### Test fit

- Apply the tuned model to the test dataset 

```{r}
test_fit <- finalize_tree %>%
  fit(test_x_class %>%
    bind_cols(tibble(target = test_y_class)))

evaluate_class(test_fit)
```

### Boosting (XGboost)

#### parsnip 

- Build a model 

1. Specify a model 
2. Specify an engine 
3. Specify a mode 

```{r}
# workflow
xg_wf <- workflow() %>% add_formula(target ~ .)

# spec
xg_spec <- boost_tree(

  # Mode
  mode = "classification",

  # Tuning hyperparameters

  # The number of trees to fit, aka boosting iterations
  trees = c(100, 300, 500, 700, 900),
  # The depth of the decision tree (how many levels of splits).
  tree_depth = c(1, 6),
  # Learning rate: lower means the ensemble will adapt more slowly.
  learn_rate = c(0.0001, 0.01, 0.2),
  # Stop splitting a tree if we only have this many obs in a tree node.
  min_n = 10L
) %>%
  set_engine("xgboost")

xg_wf <- xg_wf %>% add_model(xg_spec)
```

- Fit a model

```{r}
xg_fit <- xg_wf %>% fit(train_x_class %>% bind_cols(tibble(target = train_y_class)))
```

#### yardstick 

```{r}
metrics <- metric_set(
  yardstick::accuracy,
  yardstick::precision,
  yardstick::recall
)

evaluate_class(xg_fit)
```

```{r}
xg_fit_viz_metr <-
  visualize_class_eval(xg_fit)

xg_fit_viz_metr
```

- Visualize the confusion matrix. 
  
```{r}
xg_fit_viz_mat <-
  visualize_class_conf(xg_fit)

xg_fit_viz_mat
```

#### tune 

##### tune ingredients 

- We focus on the following hyperparameters: `trees,` `tree_depth,` `learn_rate,` `min_n,` `mtry,` `loss_reduction,` and `sample_size`

```{r}
tune_spec <-
  xg_spec <- boost_tree(

    # Mode
    mode = "classification",

    # Tuning hyperparameters

    # The number of trees to fit, aka boosting iterations
    trees = tune(),
    # The depth of the decision tree (how many levels of splits).
    tree_depth = tune(),
    # Learning rate: lower means the ensemble will adapt more slowly.
    learn_rate = tune(),
    # Stop splitting a tree if we only have this many obs in a tree node.
    min_n = tune(),
    loss_reduction = tune(),
    # The number of randomly selected hyperparameters
    mtry = tune(),
    # The size of the data set used for modeling within an iteration
    sample_size = tune()
  ) %>%
  set_engine("xgboost")

# Space-filling hyperparameter grids
xg_grid <- grid_latin_hypercube(
  trees(),
  tree_depth(),
  learn_rate(),
  min_n(),
  loss_reduction(),
  sample_size = sample_prop(),
  finalize(mtry(), train_x_class),
  size = 30
)

# 10-fold cross-validation

set.seed(1234) # for reproducibility

xg_folds <- vfold_cv(train_x_class %>% bind_cols(tibble(target = train_y_class)),
  strata = target
)
```

##### Add these elements to a workflow 

```{r}
# Update workflow
xg_wf <- xg_wf %>% update_model(tune_spec)

# Tuning results
xg_res <- xg_wf %>%
  tune_grid(
    resamples = xg_folds,
    grid = xg_grid,
    control = control_grid(save_pred = TRUE)
  )
```

##### Visualize 

```{r}

conflict_prefer("filter", "dplyr")

xg_res %>%
  collect_metrics() %>%
  filter(.metric == "roc_auc") %>%
  pivot_longer(mtry:sample_size,
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(x = value, y = mean, color = parameter)) +
  geom_point(alpha = 0.8, show.legend = FALSE) +
  facet_wrap(~parameter, scales = "free_x") +
  labs(
    y = "AUC",
    x = NULL
  )
```

```{r}
# Optimal hyperparameter
best_xg <- select_best(xg_res, "roc_auc")

best_xg

# Add the hyperparameter to the workflow
finalize_xg <- xg_wf %>%
  finalize_workflow(best_xg)
```

```{r}
xg_fit_tuned <- finalize_xg %>%
  fit(train_x_class %>% bind_cols(tibble(target = train_y_class)))

# Metrics
(xg_fit_viz_metr + labs(title = "Non-tuned")) / (visualize_class_eval(xg_fit_tuned) + labs(title = "Tuned"))

# Confusion matrix
(xg_fit_viz_mat + labs(title = "Non-tuned")) / (visualize_class_conf(xg_fit_tuned) + labs(title = "Tuned"))
```

- Visualize variable importance 

```{r}
xg_fit_tuned %>%
  pull_workflow_fit() %>%
  vip::vip()
```

##### Test fit

- Apply the tuned model to the test dataset 

```{r}
test_fit <- finalize_xg %>%
  fit(test_x_class %>% bind_cols(tibble(target = test_y_class)))

evaluate_class(test_fit)
```

### Stacking (SuperLearner)

This stacking part of the book heavily relies on [Chris Kennedy's notebook](https://github.com/dlab-berkeley/Machine-Learning-in-R/blob/master/07-ensembles.Rmd).

#### Overview

##### Stacking

Wolpert, D.H., 1992. [Stacked generalization](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.56.1533). *Neural networks*, 5(2), pp.241-259.

Breiman, L., 1996. [Stacked regressions]((https://statistics.berkeley.edu/sites/default/files/tech-reports/367.pdf). *Machine learning*, 24(1), pp.49-64.

##### SuperLearner 

The ["SuperLearner" R package](https://cran.r-project.org/web/packages/SuperLearner/index.html) is a method that simplifies ensemble learning by allowing you to simultaneously evaluate the cross-validated performance of multiple algorithms and/or a single algorithm with differently tuned hyperparameters. This is a generally advisable approach to machine learning instead of fitting single algorithms.

Let's see how the four classification algorithms you learned in this workshop (1-lasso, 2-decision tree, 3-random forest, and 4-gradient boosted trees) compare to each other and also to 5-binary logistic regression (`glm`) and the 6-mean of Y as a benchmark algorithm, in terms of their cross-validated error!

A "wrapper" is a short function that adapts an algorithm for the SuperLearner package. Check out the different algorithm wrappers offered by SuperLearner:

#### Choose algorithms

```{r}
# Review available models
SuperLearner::listWrappers()
```

```{r}
# Compile the algorithm wrappers to be used.
sl_lib <- c(
  "SL.mean", # Marginal mean of the outcome ()
  "SL.glmnet", # GLM with lasso/elasticnet regularization
  "SL.rpart", # Decision tree
  "SL.ranger", # Random forest
  "SL.xgboost"
) # Xgbboost
```

#### Fit model

Fit the ensemble!

```{r}
# This is a seed that is compatible with multicore parallel processing.
# See ?set.seed for more information.
set.seed(1, "L'Ecuyer-CMRG")

# This will take a few minutes to execute - take a look at the .html file to see the output!
cv_sl <- SuperLearner::CV.SuperLearner(
  Y = as.numeric(as.character(train_y_class)),
  X = train_x_class,
  family = binomial(),
  # For a real analysis we would use V = 10.
  cvControl = list(V = 5L, stratifyCV = TRUE),
  SL.library = sl_lib,
  verbose = FALSE
)
```

#### Risk

Risk is the average loss, and loss is how far off the prediction was for an individual observation. The lower the risk, the fewer errors the model makes in its prediction. SuperLearner's default loss metric is squared error $(y_{actual} - y_{predicted})^2$, so the risk is the mean-squared error (just like in ordinary least *squares* regression). View the summary, plot results, and compute the Area Under the ROC Curve (AUC)!

##### Summary 

* `Discrete SL` chooses the best single learner (in this case, `SL.glmnet` or `lasso`).
* `SuperLearner` takes a weighted average of the **models** using the coefficients (importance of each learner in the overall ensemble). Coefficient 0 means that learner is not used at all.
* `SL.mean_All` (the weighted mean of $Y$) is a benchmark algorithm (ignoring features). 

```{r}
summary(cv_sl)
```

##### Plot

```{r cvsl_review}
# Plot the cross-validated risk estimate with 95% CIs.

plot(cv_sl)
```

#### Compute AUC for all estimators

**ROC**

ROC: a ROC (receiver operating characteristic curve) plots the relationship between True Positive Rate (Y-axis) and FALSE Positive Rate (X-axis). 

![Area Under the ROC Curve](https://developers.google.com/machine-learning/crash-course/images/AUC.svg)

**AUC** 

AUC: Area Under the ROC Curve 

1 = perfect 

0.5 = no better than chance 

```{r}

ck37r::auc_table(cv_sl)
```

##### Plot the ROC curve for the best estimator (DiscretSL)

```{r}
ck37r::plot_roc(cv_sl)
```

##### Review weight distribution for the SuperLearner

```{r}
print(ck37r::cvsl_weights(cv_sl), row.names = FALSE)
```

The general stacking approach is available in the tidymodels framework through [`stacks`](https://github.com/tidymodels/stacks) package (developmental stage). 

However, SuperLearner is currently not available in the tidymodels framework. You can easily build and add a parsnip model if you'd like to. If you are interested in knowing more about it, please look at [this vignette](https://www.tidymodels.org/learn/develop/models/) of the tidymodels.

### Applications 

#### Bandit algorithm (optimizing an experiment)

#### Causal forest (estimating heterogeneous treatment effect)

## Unsupervised learning

x -> f - > y (not defined)

### Dimension reduction

![Projecting 2D-data to a line (PCA). From vas3k.com](https://i.stack.imgur.com/Q7HIP.gif)

#### Correlation analysis 

This dataset is a good problem for PCA as some features are highly correlated. 

Again, think about what the dataset is about. The following data dictionary comes from [this site](http://rstudio-pubs-static.s3.amazonaws.com/24341_184a58191486470cab97acdbbfe78ed5.html).

* age - age in years
* sex - sex (1 = male; 0 = female)
* cp - chest pain type (1 = typical angina; 2 = atypical angina; 3 = non-anginal pain; 4 = asymptomatic)
* trestbps - resting blood pressure (in mm Hg on admission to the hospital)
* chol - serum cholestoral in mg/dl
* fbs - fasting blood sugar > 120 mg/dl (1 = true; 0 = false)
* restecg - resting electrocardiographic results (0 = normal; 1 = having ST-T; 2 = hypertrophy)
* thalach - maximum heart rate achieved
* exang - exercise induced angina (1 = yes; 0 = no)
* oldpeak - ST depression induced by exercise relative to rest
slope - the slope of the peak exercise ST segment (1 = upsloping; 2 = flat; 3 = downsloping)
* ca - number of major vessels (0-3) colored by flourosopy
* thal - 3 = normal; 6 = fixed defect; 7 = reversable defect
* num - the predicted attribute - diagnosis of heart disease (angiographic disease status) (Value 0 = < 50% diameter narrowing; Value 1 = > 50% diameter narrowing)

```{r}
data_original %>%
  select(-target) %>%
  corrr::correlate()
```

#### Descriptive statistics 

Notice the scaling issues? PCA is not scale-invariant. So, we need to fix this problem.

```{r}
min_max <- list(
  min = ~ min(.x, na.rm = TRUE),
  max = ~ max(.x, na.rm = TRUE)
)

data_original %>%
  select(-target) %>%
  summarise(across(where(is.numeric), min_max))
```

#### Preprocessing 

`recipe` is essential for preprocessing multiple features at once.

```{r}
pca_recipe <- recipe(~., data = data_original) %>%
  # Imputing NAs using mean
  step_impute_mean(all_predictors()) %>%
  # Normalize some numeric variabless
  step_normalize(c("age", "trestbps", "chol", "thalach", "oldpeak"))
```

#### PCA analysis 

```{r}
pca_res <- pca_recipe %>%
  step_pca(all_predictors(),
    id = "pca"
  ) %>% # id argument identifies each PCA step
  prep()

pca_res %>%
  tidy(id = "pca")
```

##### Screeplot

```{r}
# To avoid conflicts
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")

pca_recipe %>%
  step_pca(all_predictors(),
    id = "pca"
  ) %>% # id argument identifies each PCA step
  prep() %>%
  tidy(id = "pca", type = "variance") %>%
  filter(terms == "percent variance") %>%
  ggplot(aes(x = component, y = value)) +
  geom_col() +
  labs(
    x = "PCAs of heart disease",
    y = "% of variance",
    title = "Scree plot"
  )
```

##### View factor loadings 

Loadings are the covariances between the features and the principal components (=eigenvectors).

```{r}
pca_recipe %>%
  step_pca(all_predictors(),
    id = "pca"
  ) %>% # id argument identifies each PCA step
  prep() %>%
  tidy(id = "pca") %>%
  filter(component %in% c("PC1", "PC2")) %>%
  ggplot(aes(
    x = fct_reorder(terms, value), y = value,
    fill = component
  )) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(
    x = "Terms",
    y = "Contribtutions",
    fill = "PCAs"
  )
```

**The key lesson**

You can use these low-dimensional data to solve the curse of dimensionality problem. Compressing feature space via dimension reduction techniques is called feature extraction. PCA is one way of doing this.

### Topic modeling 

#### Setup 

```{r}
pacman::p_load(
  tidytext, # tidy text analysis
  glue, # paste string and objects
  stm, # structural topic modeling
  gutenbergr
) # toy datasets
```

#### Dataset 

The data munging process draws on [Julia Silge's blog post](https://juliasilge.com/blog/sherlock-holmes-stm/).

```{r}
sherlock_raw <- gutenberg_download(1661)

sherlock <- sherlock_raw %>%
  # Mutate story using a conditional statement
  mutate(
    story = ifelse(str_detect(text, "ADVENTURE"), text, NA)
  ) %>%
  # Fill in missing values with next value
  tidyr::fill(story, .direction = "down") %>%
  # Filter
  dplyr::filter(story != "THE ADVENTURES OF SHERLOCK HOLMES") %>%
  # Factor
  mutate(story = factor(story, levels = unique(story)))

sherlock <- sherlock[, 2:3] # no id
```

#### Key ideas 

![Source: paperswithcode.com](https://paperswithcode.com/media/thumbnails/task/task-0000000179-fd3a1d11_fGQkZCJ.jpg)

-   Main papers: See [Latent Dirichlet Allocation](https://proceedings.neurips.cc/paper/2001/file/296472c9542ad4d4788d543508116cbc-Paper.pdf) by David M. Blei, Andrew Y. Ng and Michael I. Jordan (then all Berkeley) and this [follow-up paper](http://www.cse.cuhk.edu.hk/irwin.king/_media/presentations/latent_dirichlet_allocation.pdf) with the same title.

-   Topics as **distributions** of words ($\beta$ distribution)

-   Documents as **distributions** of topics ($\alpha$ distribution)

-   What distributions?

    -   Probability

    -   Multinominal

-   Words lie on a lower-dimensional space (dimension reduction akin to PCA)

-   Co-occurrence of words (clustering)

-   Bag of words (feature engineering)

    -   Upside: easy and fast (also working quite well)
    -   Downside: ignored grammatical structures and rich interactions among words (Alternative: word embeddings. Please check out [text2vec](http://text2vec.org/))

-   Documents are exchangeable (sequencing won't matter).

-   Topics are independent (uncorrelated). If you don't think this assumption holds, use Correlated Topics Models by [Blei and Lafferty (2007)](https://arxiv.org/pdf/0708.3601.pdf#:~:text=The%20correlated%20topic%20model%20(CTM)%20is%20a%20hierarchical%20model%20of,are%20document%2D%20specific%20random%20variables.).

#### Exploratory data analysis 

```{r}
sherlock_n <- sherlock %>%
  unnest_tokens(
    output = word,
    input = text
  ) %>%
  count(story, word, sort = TRUE)

sherlock_total_n <- sherlock_n %>%
  group_by(story) %>%
  summarise(total = sum(n))

sherlock_words <- sherlock_n %>%
  left_join(sherlock_total_n)

sherlock_words %>%
  mutate(freq = n / total) %>%
  group_by(story) %>%
  top_n(10) %>%
  ggplot(aes(
    x = fct_reorder(word, freq),
    y = freq,
    fill = story
  )) +
  geom_col() +
  coord_flip() +
  facet_wrap(~story,
    ncol = 2,
    scales = "free_y"
  ) +
  scale_fill_viridis_d() +
  labs(
    x = "",
    fill = "Story"
  ) +
  theme(legend.position = "bottom")
```

#### STM

[Structural Topic Modeling](https://www.structuraltopicmodel.com/) by Roberts, Stewart, and Tingley helps estimate how topics' proportions vary by covariates. If you don't use covariates, this approach is close to CTM. The other useful (and very recent) topic modeling package is Keyword Assisted Topic Models ([keyATM](https://keyatm.github.io/keyATM/)) by Shusei, Imai, and Sasaki.

Also, note that we didn't cover other important techniques in topic modeling, such as dynamic and hierarchical topic modeling.

![](https://warin.ca/shiny/stm/images/fig02.png)

##### Turn text into document-term matrix

`stm` package has its preprocessing function.

```{r}
dtm <- textProcessor(
  documents = sherlock$text,
  metadata = sherlock,
  removestopwords = TRUE,
  verbose = FALSE
)
```

##### Tuning K

- K is the number of topics. 
- Let's try K = 5, 10, 15.

```{r}
test_res <- searchK(
  dtm$documents,
  dtm$vocab,
  K = c(5, 10, 15),
  prevalence = ~story,
  data = dtm$meta
)
```

##### Evaludating models 

Several metrics assess topic models' performance: the held-out likelihood, residuals, semantic coherence, and exclusivity. Here we examine the relationship between semantic coherence and exclusivity to understand the trade-off involved in selecting K.

-   Semantic coherence: high probability words for a topic co-occur in documents

-   Exclusivity: keywords of one topic are not likely to appear as keywords in other topics.

> In Roberts et al. 2014 we proposed using the Mimno et al. 2011 semantic coherence metric for helping with topic model selection. However, we found that semantic coherence alone is relatively easy to achieve by having only a couple of topics that dominate the most common words. Thus we also proposed an exclusivity measure.

> Our exclusivity measure includes some information on word frequency as well. It is based on the FREX labeling metric (calcfrex) with the weight set to .7 in favor of exclusivity by default.

```{r}
test_res$results %>%
  unnest(c(K, exclus, semcoh)) %>%
  select(K, exclus, semcoh) %>%
  mutate(K = as.factor(K)) %>%
  ggplot(aes(x = exclus, y = semcoh)) +
  geom_point() +
  geom_text(
    label = glue("K = {test_res$results$K}"),
    size = 5,
    color = "red",
    position = position_jitter(width = 0.05, height = 0.05)
  ) +
  labs(
    x = "Exclusivity",
    y = "Semantic coherence",
    title = "Exclusivity and semantic coherence"
  )
```

##### Finalize 

```{r}
final_stm <- stm(dtm$documents,
  dtm$vocab,
  K = 10, prevalence = ~story,
  max.em.its = 75,
  data = dtm$meta,
  init.type = "Spectral",
  seed = 1234567,
  verbose = FALSE
)
```

##### Explore the results 

- Using the `stm` package. 

```{r}
plot(final_stm)
```

- Using ggplot2 

In LDA distribution, $\alpha$ represents document-topic density and $\beta$ represents topic-word density. 

```{r}
# tidy
tidy_stm <- tidy(final_stm)

# top terms
tidy_stm %>%
  group_by(topic) %>%
  top_n(10, beta) %>%
  ungroup() %>%
  ggplot(aes(fct_reorder(term, beta), beta, fill = as.factor(topic))) +
  geom_col(alpha = 0.8, show.legend = FALSE) +
  facet_wrap(~topic, scales = "free_y") +
  coord_flip() +
  scale_y_continuous(labels = scales::percent) +
  scale_fill_viridis_d()
```

## References

### Books 

- *An Introduction to Statistical Learning - with Applications in R (2013)* by Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani. Springer: New York. [Amazon](https://www.amazon.com/Introduction-Statistical-Learning-Applications-Statistics/dp/1461471370) or [free PDF](http://www-bcf.usc.edu/~gareth/ISL/). 

- *Hands-On Machine Learning with R (2020)* by Bradley Boehmke & Brandon Greenwell. [CRC Press](https://www.routledge.com/Hands-On-Machine-Learning-with-R/Boehmke-Greenwell/p/book/9781138495685) or [Amazon](https://www.amazon.com/gp/product/1138495689?pf_rd_p=ab873d20-a0ca-439b-ac45-cd78f07a84d8&pf_rd_r=JBRX0ZJ1WFSR9T3JPTQE)

- *Applied Predictive Modeling (2013)* by Max Kuhn and Kjell Johnson. Springer: New York. [Amazon](https://www.amazon.com/Applied-Predictive-Modeling-Max-Kuhn/dp/1461468485?SubscriptionId=0ENGV10E9K9QDNSJ5C82&tag=apm0a-20&linkCode=xm2&camp=2025&creative=165953&creativeASIN=1461468485) 

- *Feature Engineering and Selection: A Practical Approach for Predictive Models (2019)* by Kjell Johnson and Max Kuhn. Taylor & Francis. [Amazon](http://www.feat.engineering/) or [free HTML](http://www.feat.engineering/). 
- *[Tidy Modeling with R](https://www.tmwr.org/) (2020)* by Max Kuhn and Julia Silge (work-in-progress)

### Lecture slides 

- [An introduction to supervised and unsupervised learning (2015)](https://www.nber.org/econometrics_minicourse_2015/nber_slides11.pdf) by Susan Athey and Guido Imbens 

- [Introduction Machine Learning with the Tidyverse](https://education.rstudio.com/blog/2020/02/conf20-intro-ml/) by Alison Hill

### Blog posts 

- ["Using the recipes package for easy pre-processing"](http://www.rebeccabarter.com/blog/2019-06-06_pre_processing/) by Rebecca Barter