## Random Forest for multiclass classification (tidymodels inside)

We now move on from binary to multiclass classification, and put together also the use of `tidymodels`. We use the same dataset on diabetes and metabolomics that we used for the Lasso model using `tidymodels`

```{r}
library("vip")
library("ggplot2")
library("tidyverse")
library("tidymodels")
library("data.table")
library("randomForest")
```

```{r}
mtbsl1 <- fread("../data/MTBSL1.tsv")
names(mtbsl1)[c(4:ncol(mtbsl1))] <- paste("mtbl",seq(1,ncol(mtbsl1)-3), sep = "_")
```

#### Creating the multinomial variable

We combine the variables `Gender` and `Metabolic_syndrom` to create a synthetic outcome variable with four classes:

```{r}
mtbsl1$gender_status <- paste(mtbsl1$Gender,mtbsl1$Metabolic_syndrome,sep="_")
mtbsl1 %>% group_by(gender_status) %>%
    summarise(N=n())
```

#### Data splitting

We first split the data in the training and test sets (stratifying by the categorical outcome):

```{r}
diab_dt <- select(mtbsl1, -c(`Primary ID`, Gender, Metabolic_syndrome))
mtbsl1_split <- initial_split(diab_dt, strata = gender_status, prop = 0.75)
mtbsl1_train <- training(mtbsl1_split)
mtbsl1_test <- testing(mtbsl1_split)

nrow(mtbsl1_train)
nrow(mtbsl1_test)
```

#### Preprocessing

We use tidymodels to build a recipe for data preprocessing:

-   remove correlated variables
-   remove non informative variables (zero variance)
-   standardize all variables
-   impute missing data (Random Forest does not handle missing data)

```{r}
diab_recipe <- mtbsl1_train %>%
  recipe(gender_status ~ .) %>%
  step_corr(all_predictors(), threshold = 0.9) %>%
  step_zv(all_numeric(), -all_outcomes()) %>%
  step_normalize(all_numeric(), -all_outcomes()) %>%
  step_impute_knn(all_numeric(), neighbors = 5) ## there are no missing data here, but in case!
```

```{r}
prep_diab <- prep(diab_recipe)
print(prep_diab)
```

```{r}
training_set <- juice(prep_diab)
head(training_set)
```

#### Model building

We now specify the structure of our model:

-   hyperparameters to tune: `mtry` (number of features to sample for each tree) and `min_n` (minimum number of data points in a node to allow further splitting)
-   number of trees in the forest
-   the problem at hand (classification)
-   the engine (R package)

Then we put this in a workflow together with the preprocessing recipe

```{r}
tune_spec <- rand_forest(
  mtry = tune(),
  trees = 100,
  min_n = tune()
) %>%
  set_mode("classification") %>%
  set_engine("randomForest")
```

```{r}
tune_wf <- workflow() %>%
  add_formula(gender_status ~ .) %>%
  add_model(tune_spec)
```

#### Tuning the hyperparameters

We use k-fold cross-validation to tune the hyperparameters in the training set

```{r}
trees_folds <- vfold_cv(training_set, v = 5, repeats = 3)
```

```{r}
print(trees_folds)
```

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

tune_res <- tune_grid(
  tune_wf,
  resamples = trees_folds,
  grid = 30 ## n. of tuning combinations
)

```

```{r}
print(tune_res)
```

```{r}
library("repr")
options(repr.plot.width=14, repr.plot.height=8)

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

We now try to start from $\sqrt{p}$ (classification problem)

```{r}
m <- round(sqrt(ncol(training_set)-1),0)
print(m)
rf_grid <- grid_regular(
  mtry(range = c(m-4, m+16)),
  min_n(range = c(8, 12)),
  levels = c(8,4)
)
```

```{r}
print(rf_grid)
```

```{r}
regular_res <- tune_grid(
  tune_wf,
  resamples = trees_folds,
  grid = rf_grid
)
```

```{r}
print(regular_res)
```

```{r}
regular_res %>%
  collect_metrics() %>%
  filter(.metric == "roc_auc") %>%
  mutate(min_n = factor(min_n)) %>%
  ggplot(aes(mtry, mean, color = min_n)) +
  geom_line(alpha = 0.5, size = 1.5) +
  geom_point() +
  labs(y = "AUC")
```

#### Final model

We now select the best model from the hyperparameters tuning, and fit it to the training set:

1.  selecting the best model based on AUC:

```{r}
best_auc <- select_best(regular_res, "roc_auc")
# print(best_auc)
show_best(regular_res, metric = "roc_auc", n=1)
```

```{r}
best_auc <- select_best(tune_res, "roc_auc")
# print(best_auc)
show_best(tune_res, metric = "roc_auc", n=1)
```

2.  finalise the model:

```{r}
final_rf <- finalize_model(
  tune_spec,
  best_auc
)

print(final_rf)
```

3.  finalise the workflow and fit it to the initial split (training and test data):

```{r}
final_wf <- workflow() %>%
  add_recipe(diab_recipe) %>%
  add_model(final_rf)

final_res <- final_wf %>%
  last_fit(mtbsl1_split)
```

4.  evaluate the fine-tuned RF model:

```{r}
print(final_res)
final_res %>%
  collect_metrics()
```

5.  get variable importance:

```{r}
final_res %>% 
  pluck(".workflow", 1) %>%   
  pull_workflow_fit() %>% 
  #vip(num_features = 20, geom = "point")
  vip(num_features = 25)
```

#### Predictions

We collect the predictions on the test set: for each test observations we get the probabilities of belonging to each of the four classes.

```{r}
final_res %>%
  collect_predictions()
```

```{r}
cm <- final_res %>%
  collect_predictions() %>%
  conf_mat(gender_status, .pred_class)

row.names(cm$table) <-c("FCG","FDM","MCG","MDM")
colnames(cm$table) <- c("FCG","FDM","MCG","MDM")

print(cm)
```

```{r}
autoplot(cm, type = "heatmap")
```

### Cohen's k

We briefly introduce one additional metric to evaluate multiclass classification models: **Cohen's kappa** (or $\kappa$ coefficient).

We start from the results of our RF model for the combination of gender and metabolic condition (4 classes):

```{r}
cm
```


```{r}
temp <- data.matrix(cm[[1]])

true_fcg = temp[1,1]
true_fdm = temp[2,2]
true_mcg = temp[3,3]
true_mdm = temp[4,4]

n = sum(temp) ## total n. of predictions (test size)
```

We can now calculate the total accuracy (n. of correct predictions over total number of predictions):

```{r}
accuracy = (true_fcg + true_fdm + true_mcg + true_mdm)/n
print(accuracy)
```

Is this the best metric to evaluate the performance of multiclass classification models?

We already saw with binary classification that it is important also to look at TPR and TNR:

- unbalanced data
- specific classes may be more relevant

This is especially important when you have more than two classes.

Furthermore, the total ("raw") accuracy does not consider that **correct predictions may also be obtained by chance**! 
In binary classification, we know that (if data are balanced) the chance accuracy is 0.5.

But what about multiclass classification? This is where **Cohen's kappa** comes into play. The kappa coefficient tries to consider how much better the **predictive performance is over chance accuracy**.

To do so, we need to get some measure of the **expected value for chance accuracy**.

#### Let's calculate chance accuracy!

We use the frequentist definition of probabilities:

- chance predictions (relative frequencies of predictions, per class)
- chance observations (relative class frequencies from the observed values)

We sum by rows: the predictions!

```{r}
## sum of predictions per class over total number of predictions (n)
chance_preds_fcg = sum(temp[1,])/n
chance_preds_fdm = sum(temp[2,])/n
chance_preds_mcg = sum(temp[3,])/n
chance_preds_mdm = sum(temp[4,])/n
```

```{r}
print(chance_preds_fcg)
print(chance_preds_fdm)
print(chance_preds_mcg)
print(chance_preds_mdm)
```

Now, chance observations:

```{r}
chance_obs_fcg = sum(temp[,1])/n
chance_obs_fdm = sum(temp[,2])/n
chance_obs_mcg = sum(temp[,3])/n
chance_obs_mdm = sum(temp[,4])/n
```

```{r}
print(chance_obs_fcg)
print(chance_obs_fdm)
print(chance_obs_mcg)
print(chance_obs_mdm)
```

Now we have all the ingredients to calculate chance accuracy:

```{r}
chance_accuracy = chance_preds_fcg*chance_obs_fcg + chance_preds_fdm*chance_obs_fdm + chance_preds_mcg*chance_obs_mcg + chance_preds_mdm*chance_obs_mdm
print(chance_accuracy)
```


#### Cohen's kappa

$$
\kappa = \frac{\text{accuracy} - \text{chance_accuracy}}{1- \text{chance_accuracy}}
$$

```{r}
kappa = (accuracy-chance_accuracy) / (1 - chance_accuracy)
print(kappa)
```


- the numerator calculates the difference between accuracy and chance accuracy
- if accuracy = 1, we have perfect predictive performance (the confusion matrix is diagonal) $\rightarrow \kappa = 1$, regardless of chance accuracy
- if accuracy = chance accuracy, $\kappa = 0 \rightarrow$ correct predictions are by chance
- if accuracy < chance accuracy, $\kappa < 0$ (negative) $\rightarrow$ accuracy is lower than what it would be by chance.
 