---
title: "💻 LSE DS202W 2023: Week 08 - Lab"
author: "MY NAME (MY CANDIDATE NUMBER)"
date: 04 March 2024
format: html
self-contained: true
---

# ⚙️ Setup

**Import required libraries:**

```{r}
# Tidyverse packages we will use
library(ggplot2)
library(dplyr)     
library(tidyr)     
library(readr)     

# Tidymodel packages we will use
library(rsample)
library(yardstick) 
library(parsnip)   
library(recipes)   
library(workflows) 
library(rpart)
library(rpart.plot)
library(tune)
library(randomForest)
library(vip)

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

#New package for boosting (xgboost)
library(xgboost)
```

**Read the data set:**

It is the dataset you've downloaded last week.

```{r}
# Modify the filepath if needed
filepath <- "data/corruption_data_2019_nomissing.csv"
corruption_data_2019 <- read_csv(filepath)
```

🧑‍🏫 **TEACHING MOMENT:**

(*Your class teacher will guide you through this section. Just run all the code chunks below together with your class teacher.*)

🎯 **ACTION POINTS:**

1.  Create the `corruption_poor` column by running the following code.

```{r}
corruption_data_2019 <- 
  corruption_data_2019 %>% 
  mutate(corruption_poor = if_else(cpi_score < 50, "poor", "not poor"),
         corruption_poor = as.factor(corruption_poor))
```

2.  Select the columns we will use to analyse the data.

```{r}

corruption_data_2019 <- 
  corruption_data_2019 %>% 
  select(corruption_poor, property_rights:education_index)

```

3.  Let's again randomly split our dataset into a training set (containing 70% of the rows in our data) and a test set (including 30% of the rows in our data), and retrieve our resulting training and testing sets.

```{r}
set.seed(123)
# Randomly split the initial data frame into training and testing sets (70% and 30% of rows, respectively)
split <- initial_split(corruption_data_2019, prop = 0.7, strata = corruption_poor)
training_data <- training(split)
testing_data <- testing(split)
```

Just like in last week's lab, we use stratified cross-validation (i.e a cross-validation procedure where we ensure that, for each fold, both training and test sets have the feature of interest - in our example feature `corruption_poor`- sampled in the same proportion as in the original dataset). Stratified cross-validation is particularly useful in case of class imbalance to ensure that all classes in the data are represented appropriately during the training procedure and to make sure there is no fold where one of the classes we're trying to predict is not sampled at all!

4.  Create a recipe, a decision tree model specification, and wrap them up into a workflow.

In this decision tree specification, we're tuning the `min_n` hyperparameter (i.e a parameter that refers to the minimum number of observations required in a terminal (leaf) node).

```{r}
# Create a recipe
cpi_rec <-
  recipe(corruption_poor ~ .,
         data = training_data)

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

wflow <- 
  workflow() %>% 
  add_recipe(cpi_rec) %>% 
  add_model(dt_spec)
```

5.  The bias-variance tradeoff is an essential concept to consider when choosing a model. Bias is a concept describing the difference between the model’s average prediction and the expected value. A model with high bias is said to be underfitting to the data. It makes simplistic assumptions about the training data, which makes it difficult to learn the underlying pattern. Variance captures the generalisability of the model. High variance typically means that we are overfitting to our training data, finding patterns and complexity that are a product of randomness as opposed to some real trend. Ideally, we are looking for a model with low bias, and low variance.

In order to understand the bias-variance tradeoff, we will be varying a hyperparameter of the decision tree to produce more complex decision tree models for comparison. We will do this using the tune_grid() function from the tune package, which uses grid search to train different models based on the parameter values you have chosen. The hyperparameter we are varying here is the `min_n` - a parameter that refers to the minimum number of observations required in a terminal (leaf) node.

```{r}
set.seed(234)
folds <- vfold_cv(training_data, v = 5)

ctrl <- control_grid(verbose = FALSE, save_pred = TRUE)

# Create a grid specifying the min number we want to try
grid_search <- expand_grid(
  min_n = c(2:15)
)

# This will take a little while
dt_res <- tune_grid(
  wflow,
  # This computes k-fold CV during tuning
  resamples = folds,
  grid = grid_search,
  # Making sure we keep the out-of-sample predictions for each resample during tuning
  control = control_grid(verbose = FALSE, save_pred = TRUE)
)
```

If we want to find out which parameter created the best model, we can run the following command.

```{r}
show_best(dt_res, metric = "roc_auc")
```

6.  Plot the values for the roc_auc metric for each value of `min_n` we have tried with the grid search.

```{r}
g <- dt_res %>% autoplot(metric="roc_auc")

# Customise the look of the plot
g + 
  theme_bw(base_size=20) + 
  labs(y = "ROC-AUC" , title="Decision Tree performance (training set)") +
  ylim(0.8, 1) +  
  scale_x_continuous(breaks = seq(2, 15, 1))
  
```

Note that we set the y-scale of our AUC curve plot as `[0.8-1]` to better visualise the fluctuations of AUC as the `min_n` hyperparameter varies: AUC values don't below 0.8 for decision tree models on the training dataset (as specified above) regardless of `min_n` value. To make the rest of the AUC plots comparable to this one, we'll be using the same y-scale for all AUC plots from now on.

7.  Now visualise the out-of-sample predictions for each value of `min_n` we have tried - e.g. the average validation set predictions.

```{r}
g <- 
  collect_predictions(dt_res)  %>%
  group_by(min_n) %>%
  roc_auc(truth=corruption_poor, .pred_poor, event_level="second") %>% # Remember to set event_level="second" (predicted class = poor)
  ggplot() +
  geom_point(aes(x=min_n, y=.estimate)) +
  geom_line(aes(x=min_n, y=.estimate))

# Customise the look of the plot
g + 
  theme_bw(base_size=20) + 
  labs(y = "ROC-AUC", title="Decision Tree performance (test set)") +
  ylim(0.8, 1) + 
  scale_x_continuous(breaks = seq(2, 15, 1))
```

🗣️ **DISCUSSION:**

We trained different models using a variety of values for our `min_n` parameter. Looking at the training vs testing performance at each value, how do you think the parameter changes the overfitting or underfitting of the model?

## Part II - Comparing models (30 min)

🧑‍🏫 **TEACHING MOMENT:** Your class teacher will briefly explain the concept of decision boundaries and how different models produce different decision boundaries.

So far we have looked at decision trees. But how well do other models compare to predicting a country's level of corruption?

Create a new recipe, model specification, and workflow for the support vector machine model. Note that this is using a linear kernel, SVM models can also be used with polynomial, radial and sigmoid kernels.

```{r}
# SVMs do not accept NAs or categorical variables
# Create the specification of a support vector machine (SVM) model
cpi_rec <-
  recipe(corruption_poor ~ ., data = training_data) %>%
  # What if there are missing values? We input the missing values using the k-nearest neighbors algorithm
  step_impute_knn(all_predictors(), neighbors = 3)

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

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

Have a look at the recipe here: it includes a rather unfamiliar step `step_impute_knn`. What does that correspond to? Your data could contain missing values (and the test data actually does in the `education_index` column!). Instead of potentially distorting the distribution of the data by simply removing the missing values, we perform missing value imputation, i.e we're trying to estimate the values of the blanks.

One of the methods for doing this is the k-nearest neighbours method, which estimates missing values based on the characteristics of similar neighboring data points. In practice, this means that:

-   the subset of data points that have non-missing values for the target variable are identified (i.e if the target variable is `education_index` , the subset of data points which have a value for `education_index` are identified)

-   calculate the Euclidean distance between the data point with the missing value and each data point in the subset that was identified in the previous step. The Euclidean distance calculation ignores the feature that has a missing value and only uses the features that have complete values (e.g if `education_index` is the feature with missing value, `education_index` is not taken into account in the Euclidean distance calculations)

-   select the K closest neighbours (K is a parameter passed to the imputation algorithm and gives the number of neighbours to take into account in the final calculation) based on the Euclidean distance previously calculated.

-   The final imputed value is obtained as follows:

    $$
    x_j=\frac{\sum_{i=1}^K w_i v_i}{\sum_{i=1}^K w_i}
    $$

    where

    $$
    \begin{split}
    &K \quad \text{number of neighbours} \\
    &w_i \quad \text{Euclidian distance between i-th number and data point with missing value}\\
    &v_i \quad \text{value of feature containing missing value in i-th neighbour}
    \end{split}
    $$

The `step_impute_knn` step performs (as you can expect from the name), k-nearest neighbour imputation (also abbreviated as KNN imputation). And in our case, we choose K to be 3 and choose to impute all predictors, i.e to try and estimate missing values in all predictor columns.

Now that you have a workflow to fit an SVM model:

1.  Fit the model to a testing and training set, calculating an appropriate metric (e.g. confusion matrix, ROC/AUC curve). How does it compare to the decision tree?

```{r}
model <- 
  wflow_svm %>% 
  fit(training_data)

model %>% 
  augment(new_data =testing_data) %>% 
  f_meas(truth = corruption_poor, estimate = .pred_class)
```

2.  Train an SVM model, using a radial kernel.\

```{r}
svm_spec <- 
  svm_rbf() %>%
  set_mode("classification") %>% 
  set_engine("kernlab")

wflow_svm <- 
  workflow() %>% 
  add_recipe(cpi_rec) %>% 
  add_model(svm_spec)

model <- 
  wflow_svm %>% 
  fit(training_data)

model %>% 
  augment(new_data = testing_data) %>% 
  f_meas(truth = corruption_poor, estimate = .pred_class)

```

## Part III - Ensemble methods (40 min)

🧑‍🏫 **TEACHING MOMENT:** Your class teacher will briefly explain the concept of ensemble methods.

While we can tweak hyperparameters to reduce overfitting and underfitting to try and improve the bias-variance tradeoff in decision trees, we also have techniques such as 'ensemble methods' which can also help to improve modelling results. Ensemble learning helps to improve predictions by combining several models, which can lead to better predictive performance than compared to a single model. The basic idea is to learn a set of classifiers (experts) and to allow them to vote.

Ensemble algorithms, such as bootstrapping aggregation (bagging) and boosting, which aim to reduce variance at the small cost of bias in decision trees.

We suggest you have these tabs open in your browser:

1.  The [`tidymodels` documentation](https://www.tidymodels.org) page (you can open tabs with documentation pages for each package if you need to)
2.  The [`tidyverse` documentation](https://www.tidyverse.org/) page (you can open tabs with documentation pages for each package if you need to)

This is a model specification for boosting decision trees (this one is for XGBoost).

```{r}
boost_spec <- boost_tree(trees = 200, tree_depth = 4) %>%
  set_engine("xgboost") %>%
  set_mode("classification")
```

1.  Apply boosting to our dataset using a workflow.\

```{r}
cpi_rec <-
  recipe(corruption_poor ~ .,
         data = training_data) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) 


wflow_boost <- 
  workflow() %>% 
  add_recipe(cpi_rec) %>% 
  add_model(boost_spec)

boost_model <- 
  wflow_boost %>% 
  fit(data = training_data)

boost_model %>% 
  augment(training_data) %>%
  conf_mat(truth = corruption_poor, .pred_class) %>%
  autoplot(type = "heatmap")

boost_model %>% 
  augment(testing_data) %>%
  conf_mat(truth = corruption_poor, .pred_class) %>%
  autoplot(type = "heatmap")
```

2.  How does this compare to decision trees?

3.  Tune the hyperparameters using a grid search to improve your model.

```{r}
 rec_boost <-
  recipe(corruption_poor ~ .,
         data = training_data) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE)

# Add the tune() function to the parameter you want to tune
boost_spec <- 
  boost_tree(trees = tune(), tree_depth = 4) %>%
  set_engine("xgboost") %>%
  set_mode("classification")

wflow_boost <- 
  workflow() %>% 
  add_recipe(rec_boost) %>% 
  add_model(boost_spec)

set.seed(234)
folds <- vfold_cv(training_data, v = 5,strata=corruption_poor)

ctrl <- control_grid(verbose = FALSE, save_pred = TRUE)

# Create a grid specifying the min number we want to try
grid_search = expand_grid(
  trees = floor(seq(1, 30, length.out = 5))
)

# This will take a little while
boost_res <- tune_grid(
  wflow_boost,
  resamples = folds,
  grid = grid_search,
  control = ctrl
)
```

```{r}
g <- boost_res %>% autoplot(metric="roc_auc")

# Customise the look of the plot
g + 
  theme_bw(base_size=20) + 
  labs(y = "ROC-AUC") +
  ylim(0.8, 1) + 
  scale_x_continuous(breaks = floor(seq(1, 30, length.out = 5)))
```

```{r}
# Show best model
show_best(boost_res, metric="roc_auc")
```

```{r}
# The model with the best performing value of n can be selected using select_best()
best_complexity <- select_best(boost_res, metric = "roc_auc")
boost_final <- finalize_workflow(wflow_boost, best_complexity)
boost_final_fit <- fit(boost_final, data = training_data)
boost_final_fit
```

As expected, boosted trees consistently outperform decision trees: a plateau of AUC performance (AUC of 1 based on cross-validation results!) is reached for a number of trees set between 8 and 22, with the value of AUC slightly dipping when the number of trees is either set to 1 (0.957) or 30 (0.998) but remaining above the best AUC value achieved by the decision tree model (AUC of 0.949 for tree numbers set 8, 9 or 10).

## (Bonus)

4.  Apply bagging to decision trees to try and reduce overfitting. Can you also do it while tuning hyperparameters?

```{r}

cpi_rec <-
  recipe(corruption_poor ~ .,
         data = training_data) %>%
  step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
  step_impute_knn(all_predictors(),skip=FALSE)

  # Create the specification of a model
bagging_spec <- 
   rand_forest(mtry = .cols()) %>%
  set_engine("randomForest", importance = TRUE) %>%
  set_mode("classification")

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

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

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

conf_mat_training %>%
  summary()

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

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

conf_mat_testing %>% summary()

# Uncomment the following line to see the confusion matrix
# conf_mat_testing %>% autoplot(type = "heatmap") 
```

What if we tune the `trees` and `min_n` hyperparameters?

```{r}

tune_spec <- rand_forest(
  mtry = 13,
  trees = tune(),
  min_n = tune()
) %>%
  set_mode("classification") %>%
  set_engine("randomForest", importance=TRUE)

tune_wf <- workflow() %>%
  add_recipe(cpi_rec) %>%
  add_model(tune_spec)

set.seed(456)
folds <- vfold_cv(training_data, v = 5, strata=corruption_poor)

ctrl <- control_grid(verbose = FALSE, save_pred = TRUE)

# Create a grid specifying the min number we want to try
grid_search = expand_grid(
  trees = seq(1,25, length.out = 5),
  min_n = seq(1,20)
)

# This will take a little while
rand_res <- tune_grid(
  tune_wf,
  resamples = folds,
  grid = grid_search,
  control = ctrl
)
```

```{r}
g <- rand_res %>% autoplot(metric="roc_auc")

# Customise the look of the plot
g + 
  theme_bw(base_size=20) + 
  theme(legend.position = "right") + 
  labs(y = "ROC-AUC") +
  ylim(0.8, 1) + 
  scale_x_continuous(breaks = seq(0,20, 2), limits = c(1, 20))
```

```{r}
# Show best model
show_best(rand_res)
```

Based on these metrics, random forest comes second best after boosted trees and before decision trees for this particular data set...

```{r}
# The model with the best performing value of n can be selected using select_best()
best_complexity <- select_best(rand_res, metric="roc_auc")
rand_final <- finalize_workflow(tune_wf, best_complexity)
rand_final_fit <- fit(rand_final, data = training_data)
rand_final_fit
```