---
title: "Stat. 652 Midterm"
author: "Prof. Eric A. Suess"
date: today
toc: true
format: 
  html:
    embed-resources: true
---

### Midterm

For the [titanic](https://www.kaggle.com/c/titanic/data) data set try the following machine learning classification algorithms.

Use the training and test datasets from the **titanic** R package.

You should note that the *titanic_train* has the *Survived* variable and the *titanic_test* does not. So to select your best model you need to use the *titanic_train* dataset to train and test your models. So that means you will need to select a training dataset from *titanic_train* and select a testing dataset (this would be a validation dataset) from *titanic_train* to evaluate the models you try.

I have not demonstrated the use of cross-validation, once you are comfortable running all of the models see if you can figure out how to use cross-validation to pick the best model.

Once you have picked the best model you should do the following:

1. Re-run your chosen model on the full *titanic_train* dataset.
2. Then produce predictions for the *titanic_test* dataset. This is what you would submit in a .csv to kaggle in a competition.

Build **classification models** for the *Survived* variable. Pick a model scoring function and determine which model is the best. I would suggest making a confusion matrix and computing the accuracy or kappa. 

0. Null Model
1. kNN (the sample code given did not scale or normalize, if you use this model you need to do that.)
2. Boosted C5.0 
3. Random Forest
4. Logistic Regression using regularization
5. Naive Bayes

Extra Credit:

Make one plot containing all of the ROC curves for the algorithms trained.

```{r}
#| message: false
library(tidyverse)
library(tidymodels)
```

# Data

## Loading

```{r}
library(titanic)
data(titanic_train)
data(titanic_test)
```

## Exploration

```{r}
library(naniar)
naniar::gg_miss_var(titanic_train)
naniar::gg_miss_var(titanic_test)
```

## Prep

```{r}
titanic_train <- titanic_train |>
    mutate(
        Survived = as.factor(Survived),
        Sex = as.factor(Sex),
        Embarked = as.factor(Embarked),
        Pclass = as.factor(Pclass),
        SibSp = as.double(SibSp),
        Parch = as.double(Parch)
    )
titanic_test <- titanic_test |>
    mutate(
        Sex = as.factor(Sex),
        Embarked = as.factor(Embarked),
        Pclass = as.factor(Pclass),
        SibSp = as.double(SibSp),
        Parch = as.double(Parch)
    )
```


```{r}
# Training/validation split.
set.seed(108)
parts <- titanic_train |> initial_split(prop = 0.8, strata = Survived)
train <- parts |> training()
val <- parts |> testing()
```

```{r}
rec <- recipe(Survived ~ ., data = train) |>
    update_role(PassengerId, Name, Ticket, Cabin, new_role = "id") |>
    step_impute_mean(all_numeric_predictors()) |>
    step_nzv(all_predictors()) |>
    step_normalize(all_numeric_predictors())

preprocess <- rec 
```

# Models

## Null model

The null model just predicts not survived, and has a validation accuracy of 61.5%.

```{r}
model_null <- null_model() |>
    set_engine("parsnip") |>
    set_mode("classification")

fit_null <- workflow() |>
    add_formula(Survived ~ .) |>
    add_model(model_null) |>
    last_fit(parts)

fit_null |> collect_metrics()

fit_null |> collect_predictions() |> conf_mat(truth = Survived, estimate = .pred_class)
```

## k-NN

For k nearest neighbors, after doing a grid search of k = 1 to 10 along with 10-fold cross validation, the optimal k value was found to be 7, yielding an accuracy of 79.3%.

The k nearest neighbors model found k = 8 to be the optimal v

```{r}
model_knn <- nearest_neighbor(neighbors = tune()) |>
    set_engine("kknn") |>
    set_mode("classification")

workflow_knn <- workflow() |>
    add_recipe(preprocess) |>
    add_model(model_knn)

fit_knn <- workflow_knn |>
    tune_grid(
        resamples = vfold_cv(train, v = 10, strata = Survived),
        grid = expand_grid(neighbors = 1:10),
        control = control_grid(save_pred = T)
    )

best_knn <- fit_knn |> select_best(metric = "accuracy")
best_knn

workflow_knn_final <- workflow_knn |> finalize_workflow(best_knn)

fit_knn_final <- workflow_knn_final |> last_fit(parts)

fit_knn_final |> collect_metrics()

fit_knn_final |>
    collect_predictions() |>
    conf_mat(truth = Survived, estimate = .pred_class)
```

## Boosted C5.0

Wasn't able to get C5.0 to work with tidymodels tuning.

```{r}
model_c50 <- boost_tree(trees = 10) |>
    set_engine("C5.0") |>
    set_mode("classification")

workflow_c50 <- workflow() |>
    add_recipe(preprocess) |>
    add_model(model_c50)

fit_c50 <- workflow_c50 |> fit(train)

fit_c50
```

## Random forest

The random forest model has an accuracy of 81.6%.

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

workflow_rf <- workflow() |>
    add_recipe(rec) |>
    add_model(model_rf)

fit_rf <- workflow_rf |>
    tune_grid(
        resamples = vfold_cv(train, v = 10, strata = Survived),
        grid = grid_latin_hypercube(
            mtry(range = c(1, 7)),
            trees(), min_n(),
            size = 10
        ),
        control = control_grid(save_pred = T)
    )

best_rf <- fit_rf |> select_best(metric = "accuracy")
best_rf

workflow_rf_final <- workflow_rf |> finalize_workflow(best_rf)

fit_rf_final <- workflow_rf_final |> last_fit(parts)

fit_rf_final |> collect_metrics()

fit_rf_final |>
    collect_predictions() |>
    conf_mat(truth = Survived, estimate = .pred_class)
```

## Logistic regression using regularization

The logistic regression model obtained an accuracy of 78.8%.

```{r}
model_lr <- logistic_reg(penalty = tune(), mixture = tune()) |>
    set_engine("glmnet") |>
    set_mode("classification")

# Extra preprocessing since logistic regression only works with numeric or dummy
# variables.
preprocess_lr <- rec |>
    step_dummy(Pclass, Sex, Embarked) 

workflow_lr <- workflow() |>
    add_recipe(preprocess_lr) |>
    add_model(model_lr)

fit_lr <- workflow_lr |>
    tune_grid(
        resamples = vfold_cv(train, v = 10, strata = Survived),
        grid = grid_latin_hypercube(penalty(), mixture(), size = 10),
        control = control_grid(save_pred = T)
    )

best_lr <- fit_lr |> select_best(metric = "accuracy")
best_lr

workflow_lr_final <- workflow_lr |> finalize_workflow(best_lr)

fit_lr_final <- workflow_lr_final |> last_fit(parts)

fit_lr_final |> collect_metrics()

fit_lr_final |>
    collect_predictions() |>
    conf_mat(truth = Survived, estimate = .pred_class)
```

## Naive Bayes

The Naive Bayes model obtained an accuracy of 78.2%.

```{r}
# Required for Naive Bayes implementation.
library(agua)
library(discrim)

model_nb <- naive_Bayes(smoothness = tune(), Laplace = tune()) |>
    set_engine("klaR") |>
    set_mode("classification")

workflow_nb <- workflow() |>
    add_recipe(preprocess) |>
    add_model(model_nb)

fit_nb <- workflow_nb |>
    tune_grid(
        resamples = vfold_cv(train, v = 10, strata = Survived),
        grid = grid_latin_hypercube(smoothness(), Laplace(), size = 10),
        control = control_grid(save_pred = T)
    )

best_nb <- fit_nb |> select_best(metric = "accuracy")
best_nb

workflow_nb_final <- workflow_nb |> finalize_workflow(best_nb)

fit_nb_final <- workflow_nb_final |> last_fit(parts)

fit_nb_final |> collect_metrics()

fit_nb_final |>
    collect_predictions() |>
    conf_mat(truth = Survived, estimate = .pred_class)
```

# Final

The model with the highest accuracy was the random forest model, so that will be used to fit on the entire training set.

**There is an error in the next code chunk.  Needs to be updated.**

```{r}
#| eval: false
model_final <- model_rf |> finalize_model(best_rf)

workflow_final <- workflow() |>
    add_recipe(preprocess) |>
    add_model(model_final)

fit_final <- workflow_final |> fit(titanic_train)

pred_final <- fit_final |> predict(new_data = preprocess |> bake(titanic_test) )

titanic_test |>
    dplyr::select(PassengerId) |>
    mutate(Survived = pred_final$.pred_class)
```

```{r}
library(yardstick)
library(ggplot2)

roc_data <- bind_rows(
  titanic_train2_null |>
    predict(titanic_train2_testing, type = "prob") |>
    bind_cols(titanic_train2_testing) |>
    roc_curve(survived, .pred_0) |>
    mutate(model = "Null"),
  titanic_train2_knn |>
    predict(titanic_train2_testing, type = "prob") |>
    bind_cols(titanic_train2_testing) |>
    roc_curve(survived, .pred_0) |>
    mutate(model = "kNN"),
  titanic_train2_boost |>
    predict(titanic_train2_testing, type = "prob") |>
    bind_cols(titanic_train2_testing) |>
    roc_curve(survived, .pred_0) |>
    mutate(model = "Boosted C5.0"),
  titanic_train2_ranger |>
    predict(titanic_train2_testing, type = "prob") |>
    bind_cols(titanic_train2_testing) |>
    roc_curve(survived, .pred_0) |>
    mutate(model = "Random Forest"),
  titanic_train2_glm |>
    predict(titanic_train2_testing, type = "prob") |>
    bind_cols(titanic_train2_testing) |>
    roc_curve(survived, .pred_0) |>
    mutate(model = "Logistic Regression"),
  titanic_train2_nb |>
    predict(titanic_train2_testing, type = "prob") |>
    bind_cols(titanic_train2_testing) |>
    roc_curve(survived, .pred_0) |>
    mutate(model = "Naive Bayes"))
```

