homework-3.Rmd

---
title: "Homework 3: Classification - Analysis of Titanic data"
author: "Narjes Mathlouthi for PSTAT 131/231"
output:
    html_document: 
      toc: yes
      toc_float: yes
      code_folding: show
      fig_caption: yes
---

```{r setup, include=FALSE, message=FALSE, warning = FALSE}
knitr::opts_chunk$set(
	echo = TRUE,
	message = FALSE,
	warning = FALSE
)
library(tinytex)
library(GGally)
library(kableExtra)
library(DataExplorer)
library(here)
library(tidymodels)
library(ISLR)
library(ISLR2)
library(tidyverse)
library(discrim) #lda
tidymodels_prefer()
```

## Classification

For this assignment, we will be working with part of a [Kaggle data set](https://www.kaggle.com/c/titanic/overview) that was the subject of a machine learning competition and is often used for practicing ML models. The goal is classification; specifically, to predict which passengers would survive the [Titanic shipwreck](https://en.wikipedia.org/wiki/Titanic).

![Fig. 1: RMS Titanic departing Southampton on April 10, 1912.](images/RMS_Titanic.jpg){width="363"}

Load the data from `data/titanic.csv` into *R* and familiarize yourself with the variables it contains using the codebook (`data/titanic_codebook.txt`).

Notice that `survived` and `pclass` should be changed to factors. When changing `survived` to a factor, you may want to reorder the factor so that *"Yes"* is the first level.

Make sure you load the `tidyverse` and `tidymodels`!

*Remember that you'll need to set a seed at the beginning of the document to reproduce your results.*

```{r, results='hide'}

titanic_raw <- read_csv(here("data", "titanic.csv"))

titanic <- titanic_raw %>% 
 mutate(pclass = as.factor(pclass)) %>% 
  mutate(survived = as.factor(survived)) %>% 
  mutate(ticket = as.numeric(ticket)) %>% 
  select(-name) 

titanic$survived <- recode_factor(titanic$survived, No = "yes", Yes = "no") 

```

### Question 1

Split the data, stratifying on the outcome variable, `survived.` You should choose the proportions to split the data into. Verify that the training and testing data sets have the appropriate number of observations. Take a look at the training data and note any potential issues, such as missing data.

```{r, results='hide'}
#set a seed before splitting to ensure that we can reproduce the results
set.seed(3435)
#Split the data, stratifying on the outcome variable, `survived.`
titanic_split <- initial_split(titanic, strata =survived, prop = 0.7)

#training set 
titanic_train <- training(titanic_split)

#testing set
titanic_test <- testing(titanic_split)

#Verify that the training and testing data sets have the appropriate number of observations
dim(titanic_train)
dim(titanic_test)


#What does missing data look like in our training set
#rowSums(is.na(titanic_train)) <= 0


```




*Why is it a good idea to use stratified sampling for this data?*

The `strata` argument makes sure that both sides of the split have roughly the same distribution for each value of strata. If a numeric variable is passed to `strata`, then it is binned and distributions are matched within bins.

### Question 2

Using the **training** data set, explore/describe the distribution of the outcome variable `survived`.

```{r, fig.tab="Figure 1.Overview of which columns contain missing data"}
titanic_train <- training(titanic_split)

#Looking at missing data 
DataExplorer::plot_missing(
  titanic_train,
  group = list(Good = 0.05, OK = 0.4, Bad = 0.8, Remove = 1),
  missing_only = FALSE,
  geom_label_args = list(),
  title = NULL,
  ggtheme = theme_gray(),
  theme_config = list(legend.position = c("bottom"))
)

```

```{r}

ggplot(titanic_train, aes(x = survived)) +
  geom_bar(width=0.5, fill = "blue", alpha= 0.3) +
  geom_text(stat='count', aes(label=stat(count)), vjust=-0.5) +
  theme_classic()


```

- From **Figure1.** Graph of missing rows per variable, we can see that `survived` observations do not contain any **Na's**

- From  **Figure2.** Graph depicting survival rates among Titanic passengers 

### Question 3

Using the **training** data set, create a correlation matrix of all continuous variables. Create a visualization of the matrix, and describe any patterns you see. Are any predictors correlated with each other? Which ones, and in which direction?



```{r}
ggcorr(titanic_train,
       nbreaks = 6,
       label = TRUE,
       label_size = 3,
       color = "grey50")

```

**Figure3. Correlation Matrix**

- # of siblings / spouses aboard the Titanic is correlated with # of parents / children aboard the Titanic
- #Ticket vs. number of parents / children aboard the Titanic are correlated
- age vs. # of siblings / spouses aboard the Titanic are negatively correlated 


### Question 4

Using the **training** data, create a recipe predicting the outcome variable `survived`. Include the following predictors: ticket class, sex, age, number of siblings or spouses aboard, number of parents or children aboard, and passenger fare.

Recall that there were missing values for `age`. To deal with this, add an imputation step using `step_impute_linear()`. Next, use `step_dummy()` to **dummy** encode categorical predictors. Finally, include interactions between:

-   Sex and passenger fare, and
-   Age and passenger fare.

You'll need to investigate the `tidymodels` documentation to find the appropriate step functions to use.




```{r recipe}
titanic_recipe <-
  recipe(survived ~ pclass + ticket + age + sib_sp +parch +fare , data = titanic_train) %>% 
  step_impute_linear() %>% 
  step_dummy() %>% 
  step_interact(~ sex:fare + age:fare) # because a*b would generate a + b + a:b step_interact takes care of the interaction only
  
```

### Question 5

Specify a **logistic regression** model for classification using the `"glm"` engine. Then create a workflow. Add your model and the appropriate recipe. Finally, use `fit()` to apply your workflow to the **training** data.

***Hint: Make sure to store the results of `fit()`. You'll need them later on.***

```{r model glm specification}
glm_spec <- logistic_reg() %>%
  set_engine("glm")%>% 
  set_mode("classification")
```



```{r workflow,echo=FALSE}
glm_workflow <- workflow() %>%
  add_recipe(titanic_recipe) %>%
  add_model(glm_spec)

log_fit <- fit(glm_workflow, titanic_train)
# log_fit %>% 
#   tidy()
```
### Assessing Model Performance

We can use the model to generate probability predictions for the training data:

```{r}
predict(log_fit, new_data = titanic_train, type = "prob")
```

Each row represents the probability predicted by the model that a given observation belongs to a given class. Notice this is redundant, because one could be calculated directly from the other, but it's useful in multiclass situations.

However, it's more useful to summarize the predicted values. We can use `augment()` to attach the predicted values to the daata, then generate a confusion matrix:

```{r}
augment(log_fit, new_data = titanic_train) %>%
  conf_mat(truth = survived, estimate = .pred_class)
```
Or we can create a visual representation of the confusion matrix:

```{r}
augment(log_fit, new_data = titanic_train) %>%
  conf_mat(truth = survived, estimate = .pred_class) %>%
  autoplot(type = "heatmap")
```

Let's calculate the accuracy of this model, or the average number of correct predictions it made on the **training** data. This is the model's **training error rate**.


```{r}
log_reg_acc <- augment(log_fit, new_data = titanic_train) %>%
  accuracy(truth = survived, estimate = .pred_class)
log_reg_acc
```


### Question 6

**Repeat Question 5**, but this time specify a linear discriminant analysis model for classification using the `"MASS"` engine.

```{r}

lda_mod <- discrim_linear() %>% 
  set_mode("classification") %>% 
  set_engine("MASS")

lda_wkflow <- workflow() %>% 
  add_model(lda_mod) %>% 
  add_recipe(titanic_recipe) 

lda_fit<-fit(lda_wkflow, titanic_train) #remove character 'name' and make tickert as numeric
  
```


### Assessing Performance

This can be done almost exactly the same way:

```{r}
predict(lda_fit, new_data = titanic_train, type = "prob")
```

We can view a confidence matrix and calculate accuracy on the **training data**:

```{r}
augment(lda_fit, new_data = titanic_train) %>%
  conf_mat(truth = survived, estimate = .pred_class) 
```

```{r}
lda_acc <- augment(lda_fit, new_data = titanic_train) %>%
  accuracy(truth = survived, estimate = .pred_class)
lda_acc
```

### Question 7

**Repeat Question 5**, but this time specify a quadratic discriminant analysis model for classification using the `"MASS"` engine.

```{r, echo= FALSE, message=FALSE, warning=FALSE, paged.print=FALSE}
discrim_regularized(frac_common_cov = 0, frac_identity = 0) %>% 
  set_engine("klaR") %>% 
  translate()

quad_mod <- discrim_regularized(frac_common_cov = 0, frac_identity = 0) %>% 
  set_engine("klaR") %>% 
  set_mode("classification")
  
  quad_wkflow <- workflow() %>% 
  add_model(quad_mod) %>% 
  add_recipe(titanic_recipe) 

qda_fit<-fit(quad_wkflow, titanic_train)
  
```
### Assessing Performance

And again:

```{r}
predict(qda_fit, new_data = titanic_train, type = "prob")
```

We can view a confidence matrix and calculate accuracy on the **training data**:

```{r}
augment(qda_fit, new_data = titanic_train) %>%
  conf_mat(truth = survived, estimate = .pred_class) 
```

```{r}
qda_acc <- augment(qda_fit, new_data = titanic_train) %>%
  accuracy(truth = survived, estimate = .pred_class)
qda_acc
```



### Question 8

**Repeat Question 5**, but this time specify a naive Bayes model for classification using the `"klaR"` engine. Set the `usekernel` argument to `FALSE`.

## Naive Bayes

Finally, we'll fit a Naive Bayes model to the **training data**. For this, we will be using the `naive_bayes()` function to create the specification and set the `usekernel` argument to `FALSE`. This means that we are assuming that the predictors are drawn from Gaussian distributions.

```{r}

nb_mod <- naive_Bayes() %>% 
  set_mode("classification") %>% 
  set_engine("klaR") %>% 
  set_args(usekernel = FALSE) 

nb_wkflow <- workflow() %>% 
  add_model(nb_mod) %>% 
  add_recipe(titanic_recipe)

# nb_fit <- fit(nb_wkflow, titanic_train)

nb_fit <- fit(nb_wkflow , data = titanic_train %>%  na.omit())
```


### Assessing Performance

And again:

```{r}
predict(nb_fit, new_data = titanic_train, type = "prob")
```

We can view a confidence matrix and calculate accuracy on the **training data**:

```{r}
augment(nb_fit, new_data = titanic_train) %>%
  conf_mat(truth = survived, estimate = .pred_class) 
```

```{r}
nb_acc <- augment(nb_fit, new_data = titanic_train) %>%
  accuracy(truth = survived, estimate = .pred_class)
nb_acc
```


### Question 9

Now you've fit four different models to your training data.

Use `predict()` and `bind_cols()` to generate predictions using each of these 4 models and your **training** data. Then use the *accuracy* metric to assess the performance of each of the four models.

Which model achieved the highest accuracy on the training data?

## Comparing Model Performance

Now we can make a table of the accuracy rates from these four models to choose the model that produced the highest accuracy on the training data:

```{r}
accuracies <- c(log_reg_acc$.estimate, lda_acc$.estimate, 
                nb_acc$.estimate, qda_acc$.estimate)
models <- c("Logistic Regression", "LDA", "Naive Bayes", "QDA")
results <- tibble(accuracies = round(accuracies,3), models = models)
kableExtra::kable(
  results %>%
  arrange(-accuracies), caption ="The logistic regression model is the model that performed the best")
```

### Question 10

Fit the model with the highest training accuracy to the **testing** data. Report the accuracy of the model on the **testing** data.

Again using the **testing** data, create a confusion matrix and visualize it. Plot an ROC curve and calculate the area under it (AUC).

How did the model perform? Compare its training and testing accuracies. If the values differ, why do you think this is so?

## Fitting to Testing Data

Since the Naive Bayes model performed slightly better, we'll go ahead and fit it to the testing data. In future weeks, we'll cover how to use cross-validation to try out different values for models' tuning parameters, but for now, this is a general overview of the process.

```{r}
predict(nb_fit, new_data = titanic_test, type = "prob")
```

We can view the confusion matrix on the **testing** data:

```{r}
augment(nb_fit, new_data = titanic_test) %>%
  conf_mat(truth = survived, estimate = .pred_class) 
```

We can also look at the **testing** accuracy. Here, we add two other metrics, sensitivity and specificity, out of curiosity:

```{r}
multi_metric <- metric_set(accuracy, sensitivity, specificity)

augment(nb_fit, new_data = titanic_test) %>%
  multi_metric(truth = survived, estimate = .pred_class)
```






### Required for 231 Students

In a binary classification problem, let $p$ represent the probability of class label $1$, which implies that $1 - p$ represents the probability of class label $0$. The *logistic function* (also called the "inverse logit") is the cumulative distribution function of the logistic distribution, which maps a real number *z* to the open interval $(0, 1)$.

### Question 11

Given that:

$$
p(z)=\frac{e^z}{1+e^z}
$$

Prove that the inverse of a logistic function is indeed the *logit* function:

$$
z(p)=ln\left(\frac{p}{1-p}\right)
$$

$$ 
p(z)=\frac{e^z}{1+e^z} \\
p({1+e^z}) = e^z \\
p\cdot 1 + p\cdot e^z = e^z \\
e^z - p\cdot e^z = p \\
e^z \cdot (1-p) = p \\
e^z = \frac{p}{1-p} \\
z(p) = \log_e\left(\frac{p}{1-p}\right) \\
z(p) = ln\left(\frac{p}{1-p}\right) 
$$

### Question 12

Assume that $z = \beta_0 + \beta_{1}x_{1}$ and $p = logistic(z)$. How do the odds of the outcome change if you increase $x_{1}$ by two? Demonstrate this.

Increasing `x` by a unit of 1 means multiply the odds by $e^_{B}$
Therefore, if we increase if you increase $x_{1}$ by two we get :

\$\$ $z = \beta_0 + \beta_{1}(x_{1}+2)$ \$\$
\$\$ $z = \beta_0 + \beta_{1}x_{1}+ 2\beta_{1}$ \$\$

Assume now that $\beta_1$ is negative. What value does $p$ approach as $x_{1}$ approaches $\infty$? What value does $p$ approach as $x_{1}$ approaches $-\infty$?

(Lecture slide 12)

As $x_{1}$ approaches $\infty$ p will approach 0

As $x_{1}$ approaches $-\infty$ p will approach 1

References:

Katie Coburn. (2022, April 18). *Course: PSTAT 131/PSTAT 231 - STAT MACHINE LEARN - Spring 2022*. <https://gauchospace.ucsb.edu/courses/course/view.php?id=16385>