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

# 📋 Lab Tasks

## Part 0

- [ ] Export your chat logs

# ⚙️ Setup

**Import required libraries:**

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

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

**Read the 2019 data set:**

It is the first brand-new dataset you've downloaded. 

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

::: {.callout-important}
We've been doing a lot of training vs testing splits in this course, primarily by hand and focused on a particular point in time (e.g., test samples are defined to start in 2019).

But not all splits need to be done this way. Some data sets don't have a time-series component and can be split randomly for testing purposes. To make things more robust, we don't just simply use one single train versus test split, but we use a technique called cross-validation and split the data into multiple train vs test splits.
:::

In this lab, we'll learn how to use cross-validation to build a more robust model for changes in the data. Let's kick things off by cleaning up our data and randomly splitting it into a training set (70% of data points/rows) and a test set (30% of data points/rows).

🧑‍🏫 **TEACHING MOMENT:**

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

Our goal in this lab is to predict the level of corruption in a given country. The Corruption Perception Index (CPI), already present in our data set (in the `cpi_score` column), sorts countries according to their perceived corruption levels. According to [@lima_predicting_2020],"[t]he index captures the assessments of domain experts on corrupt behavioral information, originating a scale from 0 to 100 where economies close to 0 are perceived as highly corrupt while economies close to 100 are perceived as less corrupt". In other words, the CPI already gives us a way to classify countries in scales of corruption. To go from a scale of 0 to 100 (which is the CPI scale), to a categorical variable (which is what we need for our classification), we simply need the following definition:

1. **Corruption class:** We define the following levels of corruption based on the CPI scale:
    - if the CPI is lower than 50, then the corruption level is `poor`
    - if the CPI is between 50 (included) and 70 (excluded), then the corruption level is `average`
    - if the CPI is higher than 70, then the corruption level is `good`


    We store the result in the `corruption_class` column

2. **Convert the `corruption_class` column from character to factor** 

🎯 **ACTION POINTS:** 

1. Create the `corruption_class` column and convert it to `factor`. 

2. Now, let's 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)

```{.r style="margin-left:2em;"}
#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)
```

3. **What is in the training and testing set?**

    To get the actual data assigned to either set, use the `rsample::training()` and `rsample::testing()` functions:

```{.r style="margin-left:2em;"}
training_data <- training(split)
testing_data <- testing(split)
```

For curiosity, let's confirm which unique countries are represented in our training and test sets and by how many data points. We can uncover this mystery by counting the unique values in the `economy` column (you could do the same per `region`).

```{.r style="margin-left:2em;"}
# tallying the number of rows per country in the training set
training_data %>% 
    group_by(economy) %>% 
    tally()

# tallying the number of rows per country in the test set
testing_data %>% 
    group_by(economy) %>% 
    tally()
```

Just how many non-empty records are there in our datasets per country?

```{.r style="margin-left:2em;"}
# tallying the number of non-empty rows per country in the training set
training_data %>%
    drop_na() %>%
    group_by(economy) %>%
    tally() %>%
    knitr::kable()

# tallying the number of non-empty rows per country in the test set
testing_data %>%
    drop_na() %>%
    group_by(economy) %>%
    tally() %>%
    knitr::kable()
```

To browse the original Ease of Doing Business dataset, including a description of the variables (though some were recoded for ease of comprehension) - included in the `Metadata` sheet of the file-, you can download the Excel file below:

<a href="/2023/winter-term/download/ease_doing_business_data.xlsx" style="margin-left:2em;" download>
    <button class="button-61">Download Ease of Doing Business original data + documentation </button>
</a>

For more information about other variables of the dataset (i.e variables related to the Economic Freedom Index, Education Index or CPI), check out the explanations in [@lima_predicting_2020].

🗣️ **DISCUSSION:**

We defined corruption (and corruption levels) in a specific way, and we will be building a model to predict corruption in a way that fits this definition. Do you think that the definition we gave of corruption is satisfactory? Is it a good modelling objective? If not, what would you do differently? (You can look at the dataset documentation for further ideas). 



## Part II - Introduction to decision trees (40 min)

🧑🏻‍🏫 **TEACHING MOMENT:** In this part, you'll be learning about a new classification model called **decision trees**, which is more suited to handle large numbers of features of varying types that potentially contain missing data than logistic regression. 


Our dataset is a real-world dataset and ticks all these boxes:

- We mentioned earlier that our goal is to build a model to predict `corruption_class`. Aside from a few variables that don't look like likely predictors (e.g. `db_year`, `economy`,`region`,`income_group`(though this one is debatable),`cpi_score`(it is correlated with the outcome variable so can't be used to predict it!)), we have a large potential number of features/predictors to choose from for our model (100! for a dataset that only has 140 rows/data points).
- Many variables include missing values
- Though we smoothed things out, the various variables had varying types (e.g., **categorical** values or **strings** or **numerical** values)

Your class teacher will explain the basic principle of a decision tree.

### As usual, we start with a recipe


For computational reasons (too long to run!), we'll be building our model with a random subset of features from the original dataset. We'll randomly choose 25 variables from the original columns of the dataset, excluding the columns `db_year`, `economy`,`region`,`income_group` and `cpi_score`.

We construct a vector `predictors`, which contains the names of our chosen predictors:

```{r}
  remove <- c("country_code","economy","region","income_group","db_year","cpi_score") # list of the columns we exclude
  predictors <- corruption_data_2019 %>% 
                colnames() %>%
                str_remove_all(., paste(remove, collapse = "|"))%>%
                .[!(.=="")] %>%
                sample(.,25)
```

Then, we write our recipe:
```{r}

formula_string <- paste("corruption_class ~", paste(predictors, collapse = " + "))
formula <- as.formula(formula_string)
impute_rec <- recipe(formula, data = training_data) %>% 
  step_impute_median(all_numeric(), -all_outcomes()) %>%
  prep()
```

Here, the only pre-processing step we perform before fitting our model is a simple median imputation in order to fill in missing values (in numeric columns) per column with the median value of their respective column.

### Now, we can fit a decision tree model on our data

You can create the model specification for a decision tree using this scaffolding code (you need to install the `rpart` library using the `install.packages("rpart")` command to be able to use this code and of course, load the `rpart` library):

```{r}
# Create the specification of a model but don't fit it yet
dt_spec <- decision_tree(mode = "classification", tree_depth = 5) %>% 
  set_engine("rpart")
```

Now that you have the model specification:

1. Fit the model to the training set using a workflow and evaluate it its performance with an appropriate metric (e.g AUC/ROC curve)

```{r}
#Delete this line and write your code here
```

2. Fit the model to the test set using a workflow and evaluate it its performance with an appropriate metric (e.g AUC/ROC curve)

```{r}
#Delete this line and write your code here
```

🤔 What happens if you choose a different set of features to train your model?

```{r}
#Delete this line and write your code here
```


### Part III - Cross-validation (30 min)

🧑🏻‍🏫 **TEACHING MOMENT:** Your class teacher will briefly explain the concept of cross-validation.


1. Question: Can you retrain your decision tree model using cross-validation?

-   Use the `initial_split()` and `training()` functions to split your data and extract your training set. Documentation: [rsample](https://rsample.tidymodels.org/reference/initial_split.html).

```{r}
#Delete this line and write your code here
```

-   Employ a `v_fold_cv()` function to the data with 10-fold cross-validation (10-fold by default).

```{r}
#Delete this line and write your code here
```

-   Fit the model using `fit_resamples()` and collect your metrics (as in the W05 lecture notebook). How does the model performance compare to previously?

```{r}
#Delete this line and write your code here
```

-   What happens when you tweak tree parameters and/or cross-validation parameters?

```{r}
#Delete this line and write your code here
```

2. Question: How does your model perform on (a subset of the) 2020 dataset?

```{r}
#Delete this line and write your code here
```