---
title: "🗓️ Week 05: Non-linear algorithms"
subtitle: "Theme: Supervised Learning"
author: "Dr. Ghita Berrada"
format: 
  html:
    from: markdown+emoji
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 05! We wrap up our exploration of supervised learning this week.

But, let's start with some basic setup as usual...

# :gear: Setup

Before you can load the libraries below, you might have to do the following package installs (feel free to skip if you already have done the installs before)

``` r
install.packages("remotes")
remotes::install_github("grantmcdermott/parttree", force = TRUE)
install.packages("vip")
install.packages("kernlab")
install.packages("LiblineaR")
install.packages("rsample")
install.packages("rpart")
install.packages("rpart.plot")
install.packages("randomForest")
```

You can now proceed to loading your libraries and defining your functions as usual.

```{r, echo=FALSE, message=FALSE, warning=FALSE}
# Tidyverse packages
library(tidyr)
library(dplyr)
library(lubridate)
library(ggplot2)
library(readr)
library(readxl)
library(broom)
library(ggsci)

# Tidymodels packages, spelled out individually for clarity
library(rsample)
library(tune)
library(parsnip)
library(yardstick)
library(recipes)
library(workflows)
library(probably)

#new libraries for decision trees
library(rpart)
library(rpart.plot)
library(parttree)

# New packages for SVM
library(LiblineaR)
library(kernlab)

# new library for feature importance
library(vip)
#new library for random forest
library(randomForest)
library(doParallel)
```

```{r}
# ggplot2 plot themes

theme_histogram <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank())
  
}

theme_boxplot <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.position = "none")
  
}

theme_line <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        legend.position = "none")
  
}

theme_bar <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        panel.grid.major.y = element_blank())
  
}

theme_scatter <- function() {
  
  theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        legend.position = "bottom")
  
}
```

# Part 1: The case for non-linear models

:question: Why do we use linear models?

-   Linear and logistic regression are a good first shot for building ML model

    -   Easy-to-interpret coefficients
    -   Intuitive (ish) ways to assess variable importance
    -   Often good out-of-the-box prediction

-   However...

    -   Assumption that the predictors are linearly related to the outcome is restrictive
    -   Sometimes, accounting for higher order polynomial relationships or other non-linear relationships can produce better model fit

Let's take a look at a quick example: the `Wage` dataset (in the `ISLR2` library) which provides income survey data for men in the central Atlantic region of the USA. Let's load the dataset and have a very quick look at it

```{r,echo=FALSE,message=FALSE,warning=FALSE}
library(ISLR2)
glimpse(Wage)
      
```

In this dataset, we are particularly interested in the relationship between `age` and `wage` .

```{r}
Wage %>%
  ggplot(aes(age,wage)) +
  geom_point(alpha = 0.2)
```

Looking at this plot, we see that the relationship is likely not linear... And in this case, we fit a polynomial of degree 4 to this data (for more details, on the code, this example and polynomial interpolations, see [here](https://emilhvitfeldt.github.io/ISLR-tidymodels-labs/07-moving-beyond-linearity.html)).

```{r}
rec_raw_poly <- recipe(wage ~ age, data = Wage) %>%
                step_poly(age, degree = 4, options = list(raw = TRUE))

lm_spec <- linear_reg() %>%
  set_mode("regression") %>%
  set_engine("lm")

raw_poly_wf <- workflow() %>%
  add_model(lm_spec) %>%
  add_recipe(rec_raw_poly)

raw_poly_fit <- fit(raw_poly_wf, data = Wage)

tidy(raw_poly_fit)

age_range <- tibble(age = seq(min(Wage$age), max(Wage$age)))

regression_lines <- bind_cols(
  augment(raw_poly_fit, new_data = age_range),
  predict(raw_poly_fit, new_data = age_range, type = "conf_int")
)
regression_lines


Wage %>%
  ggplot(aes(age, wage)) +
  geom_point(alpha = 0.2) +
  geom_line(aes(y = .pred), color = "darkgreen",
            data = regression_lines) +
  geom_line(aes(y = .pred_lower), data = regression_lines, 
            linetype = "dashed", color = "blue") +
  geom_line(aes(y = .pred_upper), data = regression_lines, 
            linetype = "dashed", color = "blue")
```

# Part 2: Time for another case study: can we predict diabetes?

As in the week 4, we look at the diabetes prediction problem. However, we'll look at a different dataset this time and a twist on the problem.

Our dataset, this time, is the CDC Diabetes Health Indicators dataset, which comes from the [UCI Machine Learning Repository](https://archive.ics.uci.edu/dataset/891/cdc+diabetes+health+indicators) where it was made public. This dataset was created for the purpose of better understanding relationship between lifestyle and diabetes in the US. It consists of demographics features, lab test results, and answers to survey questions for each patient. The target variable (`Diabetes_012`) in this dataset encodes whether a patient has diabetes (value 1), is pre-diabetic (value 2), or is healthy (value 0).

So, we're not simply predicting whether somebody has diabetes or not, we are predicting whether someone is healthy, has pre-diabetes or has diabetes: we're going from our usual binary classification problem to a multi-class problem.

## We start by loading the data

```{r}

# set seed for reproducibility 
set.seed(123)
#loading the data
diabetes_data <- read_csv("data/cdc_diabetes.csv") %>%
  group_by(Diabetes_012) %>% 
  slice_sample(prop = 0.003) %>%  # Sample a proportion of the data set for each class
  # We have used 0.3% but you can experiment depending on
  # the capability of your machine.
  ungroup()


continuous_vars <- c("Age", "BMI")
discrete_vars <- c("GenHlth", "Education", "Income","MentHlth","PhysHlth")
categorical_vars <- c("Diabetes_012", "HighBP", "HighChol", "CholCheck",
                      "Smoker", "Stroke", "HeartDiseaseorAttack",
                      "PhysActivity", "Fruits",  "Veggies",
                      "HvyAlcoholConsump", "AnyHealthcare", "NoDocbcCost",
                      "DiffWalk", "Sex")

diabetes_data <- 
  diabetes_data %>% 
  mutate(across(all_of(continuous_vars), ~ as.numeric(.x)),
         across(all_of(discrete_vars), ~ as.integer(.x)),
         across(all_of(categorical_vars), ~ as.factor(.x)))

```

🩻 **Exploration**

1.  What variables does our data contain?

```{r}
glimpse(diabetes_data)
```

The data contains `r nrow(diabetes_data)` rows and `r ncol(diabetes_data)` columns.

2.  Can we get summary statistics for the variables in our dataset? What do we learn?

```{r}
summary(diabetes_data)
```

3.  Let's, as usual, plot some descriptive distributions in our data

```{r}
diabetes_data %>% 
  ggplot(aes(Diabetes_012, Income, fill = Diabetes_012)) +
  geom_boxplot() +
  theme_boxplot() +
  scale_fill_npg() +
  labs(x = "Diabetes", y = "Income")
```

```{r}
diabetes_data %>% 
  ggplot(aes(Diabetes_012, Education, fill = Diabetes_012)) +
  geom_boxplot() +
  theme_boxplot() +
  scale_fill_npg() +
  labs(x = "Diabetes", y = "Education")
```

```{r}
diabetes_data %>% 
  ggplot(aes(Diabetes_012, BMI, fill = Diabetes_012)) +
  geom_boxplot() +
  theme_boxplot() +
  scale_fill_npg() +
  labs(x = "Diabetes", y = "BMI")
```

## Predicting diabetes with Support Vector Machines

Now that we've looked at our variables and made sense of them, let's try our hand at some prediction.

The very first model we'll try out is called Support Vector Machines (SVM). SVM are mainly used for classification, though they are sometimes also used for regression tasks.

### What is the principle of Support Vector Machines?

Consider this example from the famous Iris dataset with two easily separated classes of of Iris flowers (according to their sepal width and sepal length characteristics)...Colours in the plots represent distinct species...

```{r, dpi=300}
theme_set(ggthemes::theme_clean())

p = iris |> 
   filter(Species != "versicolor") |> 
   ggplot(aes(x = Sepal.Length, y = Sepal.Width, colour = Species)) +
   geom_point()
p
```

The goal here is to find the line (hyperplane in dimensions higher than 2) that best separates our classes (not necessarily simply two classes). There is obviously an infinity of hyperplanes separating the classes.

```{r, dpi=300}

p = p +
   geom_segment(aes(x = 4.4, y = 2, xend = 6.5, yend = 4.5), 
                colour = 5, lty = 2, alpha = 0.7) +
   geom_segment(aes(x = 5, y = 2, xend = 6.5, yend = 4), 
                colour = 6, lty = 2, alpha = 0.7) +
geom_segment(aes(x = 4, y = 2, xend = 7, yend = 4.5), 
                colour = 4, lty = 2, alpha = 0.7)
p

```

SVM chooses the one that maximises the margin (i.e the distance) between the classes.

![](svm_iris_margin.png)

Note that the separation does not have to be linear. If the classes can't be separated linearly, we apply the so-called kernel trick (i.e we transform the data with a function called a kernel, which is often equivalent to projecting the data onto a higher dimension), so as to make it separable.

### Back to our case study

First let's split our data into training and test set as usual (70% training data and 30% test set):

```{r}
set.seed(123)
split <- initial_split(diabetes_data, prop = 0.7,strata=Diabetes_012)
training_data <- training(split)
testing_data <- testing(split)
```

SVMs require a bit of pre-processing to be applied to the data before the model (they can't deal with categorical variables out-of-the box, require normalization and have issues with missing values). So let's write a recipe that addresses this first.

```{r}
svm_rec <-
  recipe( Diabetes_012 ~ .,
         data = training_data) %>%
  step_normalize(all_numeric_predictors(),-all_outcomes())%>%
  step_dummy(all_nominal_predictors(),-all_outcomes())%>%
  step_naomit(all_predictors(),-all_outcomes(), skip = FALSE)
  
```

Then we move to our usual workflow sequence (we'll use a linear SVM for now):

```{r}
# Create the specification of a support vector machine (SVM) model
svm_spec <- 
    svm_linear(mode = "classification") %>%
    set_engine("LiblineaR")

wflow_svm <- 
  workflow() %>% 
  add_recipe(svm_rec) %>% 
  add_model(svm_spec)
```

**:dart: Action Points**

1.  Can we fit the model to the training and testing set, calculating an appropriate metric (e.g. confusion matrix, ROC/AUC curve).

    ```{r}
    doParallel::registerDoParallel()

    svm_model <- 
     wflow_svm %>% 
     fit(data = training_data)


    ```

    ```{r}
    conf_mat_svm_training <- svm_model %>% 
      augment(new_data=training_data, type.predict = "response") %>%
      conf_mat(truth = Diabetes_012, .pred_class)

    conf_mat_svm_training %>%
      summary(event_level="second")

    conf_mat_svm_training %>%
      autoplot(type = "heatmap")
    ```

    The accuracy metric is to be taken with a pinch of salt here since we have a heavily imbalanced dataset (84.2% of the data is from class 0 i.e the healthy class, 1.8% to the class 1 i.e diabetes and about 14% of the data is from class 2 i.e pre-diabetes).

    The metrics we should particularly take into account are the negative predictive value (npv), which would decrease the more false negatives (i.e people who either have pre-diabetes or diabetes but are misclassified as healthy) increase and the F1-score (which provides a balanced measure of the number of sick people/healthy people correctly categorized).

    On this score, the npv metric indicates good performance while the F1-score seems to indicate a middling performance of our classifier and our confusion matrix reveals that this classifier performs extremely poorly when it comes to predicting diabetes (all diabetes instances are either misclassified as pre-diabetes or healthy), does much better at correctly identifying the healthy class and misclassfies pre-diabetes as healthy more often than not: this is all very problematic for our application obviously but could be explained by the class imbalance.

    For more on the interpretation of the confusion matrix in the multiclass case, see [here](https://www.analyticsvidhya.com/blog/2021/06/confusion-matrix-for-multi-class-classification/#:~:text=And%20TP%3A%20the%20True%2Dpositive,predicted%20value%20are%20the%20same.&text=Let%20us%20calculate%20the%20TP,1%20is%20the%20TP%20value.).

    ```{r}
    conf_mat_svm_testing <- svm_model %>% 
      augment(new_data=testing_data, type.predict = "response") %>%
      conf_mat(truth = Diabetes_012, .pred_class)

    conf_mat_svm_testing %>%
      summary(event_level="second")


    conf_mat_svm_testing %>%
      autoplot(type = "heatmap")
    ```

    The same conclusions are mirrored in the testing data.

2.  Change a parameter of the SVM model (e.g kernel) to find out if it improved the predictions (e.g. confusion matrix, AUC/ROC curve).

    We try an tweak the `cost` hyperparameter i.e the SVM parameter is that controls the trade-off between maximizing the margin and minimizing the classification error ( a smaller cost encourages a larger margin and allows more errors in the training data, while a larger cost penalizes errors more heavily and results in a smaller margin). This parameter helps to control overfitting and underfitting in SVM models.

    ```{r}
    svm_spec_cost <- 
        svm_linear(mode = "classification", cost=2) %>%
        set_engine("LiblineaR")

    wflow_svm_cost <- 
      workflow() %>% 
      add_recipe(svm_rec) %>% 
      add_model(svm_spec_cost)

    svm_model_cost <- 
      wflow_svm_cost %>% 
      fit(data = training_data)


    ```

    ```{r}
    # Evaluate the fitted model using conf matrices and summary stats
    cm_svm_cost_training <- svm_model_cost %>% 
      augment(new_data=training_data, type.predict = "response") %>%
      conf_mat(truth = Diabetes_012, .pred_class) 

    cm_svm_cost_training %>%
      summary(event_level="second")

    cm_svm_cost_training %>%
      autoplot(type = "heatmap")
    ```

    ```{r}
    cm_svm_cost_testing <- svm_model_cost %>% 
      augment(new_data=testing_data,type.predict = "response") %>%
      conf_mat(truth = Diabetes_012, .pred_class) 

    cm_svm_cost_testing %>%
      summary(event_level="second")

    cm_svm_cost_testing %>%
      autoplot(type = "heatmap")
    ```

Tweaking the cost parameter doesn't seem to change the performance of our model much in this instance. It could be due to the rather "clean" nature of our data...

```{r}
svm_spec_rbf <- 
    svm_rbf(mode = "classification",rbf_sigma = 0.2) %>%
    set_engine("kernlab")

wflow_svm_rbf <- 
  workflow() %>% 
  add_recipe(svm_rec) %>% 
  add_model(svm_spec_rbf)

svm_model_rbf <- 
  wflow_svm_rbf %>% 
  fit(data = training_data)

```

```{r}
# Evaluate the fitted model using conf matrices and summary stats
cm_svm_rbf_training <- svm_model_rbf %>% 
  augment(new_data=training_data) %>%
  conf_mat(truth = Diabetes_012, .pred_class) 

cm_svm_rbf_training %>%
  summary(event_level="second")

cm_svm_rbf_training %>%
  autoplot(type = "heatmap")
```

```{r}
cm_svm_rbf_testing <- svm_model_rbf %>% 
  augment(new_data=testing_data) %>%
  conf_mat(truth = Diabetes_012, .pred_class) 

cm_svm_rbf_testing %>%
  summary(event_level="second")

cm_svm_rbf_testing %>%
  autoplot(type = "heatmap")
```

On the other hand, changing to a non-linear kernel seems to make a massive difference in terms of performance: the F1-score and npv are both much higher! The confusion matrix for the training set shows that all instances of the healthy class are correctly identified, that the diabetes class is much better predicted here (though a couple of instances are still predicted as healthy) though it's still the class it struggles most with (3/8 entities misclassified i.e 37.5% compared to 23/75 i.e 30.7% for pre-diabetes).

The performance degrades significantly in the test set: the F1-score is middling mainly due to a low recall. If we look at the confusion matrix, we see that the model consistently misclassifies the diabetes class and does poorly too at predicting pre-diabetes.

```{r}
  svm_model_rbf %>% 
  augment(new_data= testing_data) %>%
  roc_auc(truth = Diabetes_012, .pred_0, .pred_1, .pred_2, event_level="second") 

  svm_model_rbf %>% 
  augment(testing_data) %>%
  roc_curve(truth=Diabetes_012, .pred_0, .pred_1, .pred_2, event_level="second") %>%
  autoplot() + 
  geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
  scale_color_gradient(name="Threshold", low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
  labs(title="Meet the ROC curve",
       subtitle="This model performs better than a random guess!", 
       x="(1 - specificity) = 1 - TN/N",
       y="(sensitivity) = TP/P")
```

The ROC curve for this model seems to show that it performs better than a random guess for any given threshold...

We've highlighted earlier that SVMs had a few shortcomings:

-   they can't handle missing values (so imputation or other methods of missing value handling are needed before proceeding)

-   numerical values need to be normalized if the model is to converge

-   categorical values can't be handled out-of-the-box and a dummy step is required

Another shortcoming is that they're not explainable methods, i.e it's hard to know what features exactly contributed to the final class assignments. They would require post-hoc methods like LIME (Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016. "Why Should I Trust You?": Explaining the Predictions of Any Classifier. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD '16). Association for Computing Machinery, New York, NY, USA, 1135–1144. <https://doi.org/10.1145/2939672.2939778>) or Shapley values (Scott M. Lundberg and Su-In Lee. 2017. A unified approach to interpreting model predictions. In Proceedings of the 31st International Conference on Neural Information Processing Systems (NIPS'17). Curran Associates Inc., Red Hook, NY, USA, 4768–4777.<https://proceedings.neurips.cc/paper_files/paper/2017/file/8a20a8621978632d76c43dfd28b67767-Paper.pdf>) to become explainable.

## Predicting diabetes with decision trees

Trees, on the other hand, are inherently explainable methods. But, we'll see why after fitting a tree model to our data.

### Let's fit decision trees to our data (just as during Monday's lab)

```{r}
tree_rec <-
  recipe(Diabetes_012 ~ .,
         data = training_data)

# Create the specification of a model
dt_spec <- 
    decision_tree(mode = "classification", 
    # You must specify the parameter you want to tune
    tree_depth=5) %>% 
    set_engine("rpart")

tree_wflow <- 
  workflow() %>% 
  add_recipe(tree_rec) %>% 
  add_model(dt_spec)

tree_model <- 
  tree_wflow %>% 
  fit(data = training_data)
```

Unlike with SVMs, pre-processing was minimal. Trees can handle mixed types of data (categorical data included) fairly well and are tolerant to missing values. Like many other methods, they are affected by the curse of dimensionality, i.e data of large dimensions (too many data points and/or features). We'll see, in W09, how to overcome that.

But, for now, let's evaluate our new tree model...

```{r}
# Evaluate the fitted model using conf matrices and summary stats
cm_tree_training <- tree_model %>% 
  augment(new_data= training_data) %>%
  conf_mat(truth = Diabetes_012, .pred_class) 

cm_tree_training %>%
  summary(event_level="second")

cm_tree_training %>%
  autoplot(type = "heatmap")
```

```{r}
cm_tree_testing <- tree_model %>% 
  augment(new_data=testing_data) %>%
  conf_mat(truth = Diabetes_012, .pred_class) 

cm_tree_testing %>%
  summary(event_level="second")

cm_tree_testing %>%
  autoplot(type = "heatmap")
```

```{r}
  tree_model %>% 
  augment(new_data= testing_data) %>%
  roc_auc(truth = Diabetes_012, .pred_0, .pred_1, .pred_2, event_level="second") 

  tree_model %>% 
  augment(testing_data) %>%
  roc_curve(truth=Diabetes_012, .pred_0, .pred_1, .pred_2, event_level="second") %>%
  autoplot() + 
  geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
  scale_color_gradient(name="Threshold", low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
  labs(title="Meet the ROC curve",
       subtitle="This model also performs better than a random guess and on par with SVM!", 
       x="(1 - specificity) = 1 - TN/N",
       y="(sensitivity) = TP/P")
```

On the face of it, the performance of our tree is on par with the non-linear SVM (with the non-linear SVM having a slight advantage)...

Next comes the big advantage of decision trees: the explainability part!

We can plot a tree that gives us easily interpreted rules that would allow us to make predictions on new data (and understand the predictions made on the current data!).

To understand how the tree you'll plot below works, have a look [here](https://www.spsanderson.com/steveondata/posts/2023-09-29/index.html).

```{r}
tree_fit <- tree_model %>% extract_fit_engine()
rpart.plot(tree_fit,roundint=FALSE)
```

You can also determine which features were most important in the model (see [here](https://cran.r-project.org/web/packages/vip/vip.pdf) for the full documentation) :

```{r}
tree_fit %>%
  vip(geom = "col", aesthetics = list(fill = "midnightblue", alpha = 0.8)) +
  scale_y_continuous(expand = c(0, 0))
```

From the looks of it, the three most important features in this model are: difficulty walking and BMI. Gender also comes close. Something to think about if reducing features...

## What's going on behind the scenes?

How decision trees work:

<div>

-   Divide the predictor space into $\mathbf{J}$ distinct regions $R_1$, $R_2$,...,$R_j$.
-   Take the mean of the response values in each region

</div>

Let's take a closer look through the Iris example once again:

```{r}
set.seed(456)
iris_split <- initial_split(iris)
iris_train <- training(iris_split)
iris_test <- testing(iris_split)
# extract the training data
# We will have to forego a recipe now. 
# The visualisation package we're using, parttree, is only configured to work with the rpart engine
# not with tidymodels.
simple_model <- rpart(Species ~ Petal.Length + Petal.Width, data = iris_train)

ggplot(data = iris_train, mapping=aes(x = Petal.Length, y = Petal.Width)) +
  geom_parttree(data = simple_model, aes(fill = Species), alpha=0.25) +
  geom_point(aes(color = Species), alpha = 0.7) +
  theme_minimal() +
  labs(title="An OK model",
       subtitle="The model finds three different 'boxes' that more or less match the classes")
```

How are regions created?

-   Start from the top of the tree
-   Then perform splits at a current level of depth
-   Splits are "local" not global
-   Only cares about data in the current branch

Trees apply to both regression and classification.

**Regression case**

The tree selects a predictor $X_j$ and a cutpoint $s$ that minimises the residual sum of squares. :

-   We define two half planes $R_1(j,s) = \left\{X|X_j < s\right\}$ and $R_2(j,s) = \left\{X|X_j \geq s\right\}$ and find $j$ and $s$ by minimising:

$$
\sum_{i: x_i \in R_1(j,s)} (y_i - \hat{y}_{R_1})^2 + \sum_{i: x_i \in R_2(j,s)} (y_i - \hat{y}_{R_2})^2
$$

**Classification case**

-   The tree selects a predictor $X_j$ and a cutpoint $s$ that maximises **node purity** (i.e that maximizes the consistency of groups within a node)**.**
-   Gini index (or "measure of inequality/impurity" - the lower in a node the better): $G = \sum_{k = 1}^{K} \hat{p}_{mk}(1 - \hat{p}_{mk})$
-   Entropy (or "measure of disorder/disorganization" - the lower in a node the better to achieve consistency): $D = - \sum_{k = 1}^{K} \hat{p}_{mk}\log \hat{p}_{mk}$

### A word on overfitting

What happens if I modify the `cost_complexity` parameter (i.e the parameter that sets the penalty for tree complexity and helps us prune nodes in the tree by being set higher) and set it extremely low (to 0 for example - meaning I am not pruning any node)?

```{r}
set.seed(42)

diabetes_data_original <- read_csv("data/cdc_diabetes.csv")


continuous_vars <- c("Age", "BMI")
discrete_vars <- c("GenHlth", "Education", "Income","MentHlth","PhysHlth")
categorical_vars <- c("Diabetes_012", "HighBP", "HighChol", "CholCheck",
                      "Smoker", "Stroke", "HeartDiseaseorAttack",
                      "PhysActivity", "Fruits",  "Veggies",
                      "HvyAlcoholConsump", "AnyHealthcare", "NoDocbcCost",
                      "DiffWalk", "Sex")

diabetes_data_original <- 
  diabetes_data_original %>% 
  mutate(across(all_of(continuous_vars), ~ as.numeric(.x)),
         across(all_of(discrete_vars), ~ as.integer(.x)),
         across(all_of(categorical_vars), ~ as.factor(.x)))
# Split the data into training (70%) and testing (30%) sets
data_split_original <- initial_split(diabetes_data_original, prop = 0.7)
train_data_original <- training(data_split_original)
test_data_original <- testing(data_split_original)

# Define a decision tree model 

tree_rec_original <-
  recipe(Diabetes_012 ~ .,
         data = train_data_original)

tree_spec_original <- decision_tree(cost_complexity = 0, mode="classification") %>%
  set_engine("rpart") 

# Define the workflow
tree_workflow <- workflow() %>%
  add_recipe(tree_rec_original) %>% 
  add_model(tree_spec_original) 

# Fit the model on the training data
tree_model_original <- 
  tree_workflow %>% 
  fit(data = train_data_original)

# Visualize the tree
tree_fit_original <- tree_model_original %>% extract_fit_engine()
rpart.plot(tree_fit_original,roundint=FALSE)
```

My model overfits my training data to an extreme degree and won't generalize to unseen observations at all! Decision trees are known for been prone to overfitting if their parameters are not set correctly...

Enter cross-validation...

```{r}
cv <- vfold_cv(training_data,strata=Diabetes_012)
rec <- recipe(Diabetes_012 ~ ., data = train_data_original)
w_flow <- 
  workflow() %>% 
  add_recipe(rec) %>% 
  add_model(dt_spec)
res <- fit_resamples(w_flow, cv) %>% collect_metrics()
```

## Predicting diabetes with random forest

### What is random forest?

What is the rough idea behind it?

A single tree can overfit...In the real world, when we take important decisions, we take advice from multiple sources hoping to build up consensus from the variety and number of sources.

The idea of random forest is similar. There is strength in numbers! Instead of taking one single tree, we take multiple trees and instead of having all trees build a model from the same set of features, we have each tree randomly sample a subset of features from our inputs and build a model from it. The final model is a consensus model based on all the trees built.

1.  Step 1: In a Random Forest model, a subset of data points and a subset of features is selected to construct each decision tree. In other words, $n$ data points and $m$ features are randomly selected out of the dataset to construct each of the $B$ trees that will compose the forest

2.  Step 2: Individual decision trees are constructed for each sample.

3.  Step 3: Each decision tree will generate an output.

4.  Step 4: The final output of the forest is obtained through majority voting (for classification) or averaging of tree predictions (for regression).

Check Breiman, L. Random Forests. Machine Learning 45, 5–32 (2001). <https://doi.org/10.1023/A:1010933404324> for more details

### Back to our case study

As usual, we go through the same routine of writing a model specification, workflow before fitting the model to the training data first. Note that, in this case, aside from `mtry` (i.e the number of randomly selected predictors/features per tree), we didn't specify any hyperparameters for our random forest model and stuck with default values (the number of trees is by default 500 and the `min_n` , i.e the minimum number of data points in a node that are needed for a tree node to be split further, is set to 10). However, we ideally should tune our hyperparameters to find their optimal values: you'll explore how to do that more during the W08 lab and you can also take a look at how to do that [here](https://juliasilge.com/blog/sf-trees-random-tuning/).

```{r}
bagging_spec <- 
   rand_forest(mtry = .cols()) %>%
  set_engine("randomForest", importance = TRUE) %>%
  set_mode("classification")

wflow_forest <- 
  workflow() %>% 
  add_recipe(tree_rec) %>% 
  add_model(bagging_spec)

random_forest_model <- 
  wflow_forest %>% 
  fit(data = training_data)

```

```{r}

# Evaluate the fitted random forest model using confusion matrices & summary stats
conf_mat_training <- random_forest_model %>% 
  augment(training_data) %>%
  conf_mat(truth = Diabetes_012, .pred_class) 

conf_mat_training %>%
  summary()

conf_mat_training %>%
  autoplot(type = "heatmap")

```

The performance on training data looks stunning! Does it replicate on test data?

```{r}
conf_mat_testing <- random_forest_model %>% 
  augment(testing_data) %>%
  conf_mat(truth = Diabetes_012, .pred_class) 

conf_mat_testing %>%
  summary()

conf_mat_testing %>%
  autoplot(type = "heatmap")
```

No, it doesn't, of course, but it's more than decent and it's on par with the performances of the other models on test data.

```{r}
random_forest_fit <- random_forest_model %>% extract_fit_engine()
random_forest_fit %>%
  vip(geom = "col", aesthetics = list(fill = "midnightblue", alpha = 0.8)) +
  scale_y_continuous(expand = c(0, 0))
```

Again, just like the decision tree model, the most important feature in the model is the difficulty to walk but then the next 3 features are HighBP, Stroke and BMI instead (which makes clinical sense actually!).

```{r}
  random_forest_model %>% 
  augment(testing_data) %>%
  roc_auc(truth = Diabetes_012, .pred_0,.pred_1,.pred_2, event_level="second") 

  random_forest_model %>% 
  augment(testing_data) %>%
  roc_curve(truth=Diabetes_012, .pred_0,.pred_1,.pred_2, event_level="second") %>%
  autoplot() + 
  geom_point(aes(x=1-specificity, y=sensitivity, color=.threshold)) + 
  scale_color_gradient(name="Threshold", low = "#c6733c", high="#3cc6b8", limits=c(0, 1)) + 
  labs(title="Meet the ROC curve",
       subtitle="This model performs better than a random guess, and better than both non-linear SVM and decision tree!", 
       x="(1 - specificity) = 1 - TN/N",
       y="(sensitivity) = TP/P")
```

:warning:

-   Random Forest too can overfit! They are less prone to that, though, than decision trees (since they rely average or majority ranking).

-   Slower than decision trees

-   Less interpretable!

-   Trained on a subset of features so immune to curse of dimensionality unlike decision trees

    Random Forest are a type of tree ensemble. Other types of tree ensembles include boosted trees like XGBoost (Tianqi Chen and Carlos Guestrin. 2016. XGBoost: A Scalable Tree Boosting System. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD '16). Association for Computing Machinery, New York, NY, USA, 785–794. <https://doi.org/10.1145/2939672.2939785>) (which you'll see more of during the lab on Monday) or CatBoost (Liudmila Prokhorenkova, Gleb Gusev, Aleksandr Vorobev, Anna Veronika Dorogush, and Andrey Gulin. 2018. CatBoost: unbiased boosting with categorical features. In Proceedings of the 32nd International Conference on Neural Information Processing Systems (NIPS'18). Curran Associates Inc., Red Hook, NY, USA, 6639–6649.<https://proceedings.neurips.cc/paper_files/paper/2018/file/14491b756b3a51daac41c24863285549-Paper.pdf>).
