OLDBank-Account-Fraud.Rmd

---
title: "Bank Account Fraud"
author: "Erik De Luca"
date: "`r Sys.Date()`"
output:   
  html_document:
    df_print: "paged"
    # code_folding: hide
    toc: true
    toc_color: "orchid"
    theme: united
    keep_md: true
    toc_float: true
    number_sections: true
---

# Librerie

```{r, include}
library(gbm)
library(randomForest)
library(ROSE)
library(dplyr)
library(tidyverse)
library(tree)
library(smotefamily)
library(ModelMetrics)
library(ggplot2)
library(viridis)
library(hrbrthemes)
library(kableExtra)
library(knitr)
library(tune)
library(discrim)
library(klaR) 
library(themis)
library(tidymodels)
library(plotly)
library(DescTools)
tidymodels_prefer()
conflicted::conflict_prefer("select", "dplyr")
```


# Importazione dei dati

I dati sono stati importati dal  sito Kaggle al seguente link https://www.kaggle.com/datasets/sgpjesus/bank-account-fraud-dataset-neurips-2022 e si riferiscono alle frodi bancarie.

L'obbiettivo del progetto è sviluppare un modello di previsione delle frodi bancarie.

Il dataset a disposizione offre dati realistici basati sul mondo reale ma protetti da privacy tramite diverse tecniche come *noise addittion* e addesstramento di un modello generativo *CTGAN*.

```{r}
df <- read.csv("data/Base.csv") %>% 
  tibble()

df %>% head(50)
```

## Trasformazione variabili

Trasformo la variabile `fraud_bool` da *integer* a *factor*.

```{r}
df = df %>% 
  mutate(fraud = factor(fraud_bool, labels = c("No fraud", "Fraud"))) %>%  
  select(-fraud_bool)
table(df$fraud)
```

# Analisi preliminare dei dati

Controllo della presenza di NA nel dataset.

```{r}
anyNA.data.frame(df)
```
Breve visione del dataframe, è composto da un milione di osservazioni e 32 variabili. 
Sarà necessario poner attenzione nell'ottimizzazione di alcuni comandi per la riduzione del tempo macchina nell'esecuzione degli stessi.

```{r}
df %>% 
  skimr::skim()
```





## Grafici delle variabili

Eseguo un'analisi esplorativa dei dati grafica per visualizzare con facilità eventuli problematiche o possibili miglioramenti che possono essere applicati al dataset.

### Istogramma delle variabili numeriche

Vengono prodotti gli istogrammi delle variabili numeriche per osservare le loro distribuzioni.

```{r, warning=FALSE}
df %>%
  select(where(is.numeric)) %>% 
  # mutate_all(scale) %>%  
  pivot_longer(cols = 1:9,
               names_to = "Variabili",
               values_to = "Valori") %>%  
  ggplot(aes(x = Valori, color = Variabili, fill = Variabili)) +
  geom_histogram(bins = 20, alpha = 0.6) + 
  scale_fill_viridis(discrete=TRUE) +
    scale_color_viridis(discrete=TRUE) +
    theme_ipsum() +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      strip.text.x = element_text(size = 8)
    ) +
    xlab("") +
    ylab("") +
    facet_wrap(~Variabili, ncol = 3, scales = "free")
```

```{r, warning=FALSE}
df %>%
  select(where(is.numeric)) %>% 
  # mutate_all(scale) %>%  
  pivot_longer(cols = 10:18,
               names_to = "Variabili",
               values_to = "Valori") %>%  
  ggplot(aes(x = Valori, color = Variabili, fill = Variabili)) +
  geom_histogram(bins = 20, alpha = 0.6) + 
  scale_fill_viridis(discrete=TRUE) +
    scale_color_viridis(discrete=TRUE) +
    theme_ipsum() +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      strip.text.x = element_text(size = 8)
    ) +
    xlab("") +
    ylab("") +
    facet_wrap(~Variabili, ncol = 3, scales = "free")
```

```{r, warning=FALSE}
df %>%
  select(where(is.numeric)) %>% 
  # mutate_all(scale) %>%  
  pivot_longer(cols = 19:26,
               names_to = "Variabili",
               values_to = "Valori") %>%  
  ggplot(aes(x = Valori, color = Variabili, fill = Variabili)) +
  geom_histogram(bins = 20, alpha = 0.6) + 
  scale_fill_viridis(discrete=TRUE) +
    scale_color_viridis(discrete=TRUE) +
    theme_ipsum() +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      strip.text.x = element_text(size = 8)
    ) +
    xlab("") +
    ylab("") +
    facet_wrap(~Variabili, ncol = 3, scales = "free")
```

### Grafici a barre delle variabili qualitative

#### Trasformazione Variabili

Alcune variabili devono essere trasformate da `character` a `factor`, allo stesso tempo anche alcune variabili quantitative discrete possono essere considerate come variabili qualitative ordinali. 

Prima della rappresentazione grafiche si sottopone il dataset alle trasformazioni sopra descritte.

Viene rimossa anche la variabile *device_fraud_count* in quanto possiede una varianza nulla.

Alcune variabili quantitative discrete vengono trasformate in variabili qualitative poiché presentano un numero limitato di distinte osservazioni e mostrano una correlazione lineare con la variabile dipendente. 
Tale approccio può migliorare la precisione delle stime, sebbene comporti un aumento nel numero dei parametri da stimare.

```{r}
df = df %>% 
  select(- device_fraud_count) %>% 
  mutate_at(vars(c(customer_age, income, email_is_free, proposed_credit_limit,
                   phone_home_valid, phone_mobile_valid, device_distinct_emails_8w,
                   foreign_request, has_other_cards, keep_alive_session, month)),
            ~ factor(.)) %>% # ordered = T
  mutate_if(is.character, factor)
df
```


```{r, warning=FALSE}
dfChar = df %>%
  select(where(is.factor)) 
# apply(matrix(1:(ncol(dfChar) - 1)),1, function(indexCol)
apply(matrix(1:2),1, function(indexCol)
{
  ggplot(dfChar, aes(x = get(colnames(dfChar)[indexCol]), group = fraud)) +
  geom_bar(aes(y = ..prop.., fill =  fraud), stat = "count") +
  scale_y_continuous(labels=scales::percent) +
  xlab(colnames(dfChar)[indexCol]) +
  ylab("Osservazioni") +
  theme_ipsum() +
  theme(legend.position="none")
}
)
```

## Variabile fraud con le varibili indipendenti

```{r, warning=FALSE}
df %>% 
  group_by(fraud) %>%
  summarise_if(is.numeric,mean) %>% 
  bind_cols(df %>% 
    group_by(fraud) %>%
    summarise_if(is.character,DescTools::Mode)) %>% 
  as.matrix() %>% 
  t() %>% 
  janitor::row_to_names(row_number = 1) %>% 
  as.data.frame() %>% 
  rownames_to_column("Variabili") %>% 
  tibble()
```

```{r, warning=FALSE}
grafici = lapply(1:4,function(x)
{
  df %>% 
  select(where(is.numeric)|fraud) %>%
  # mutate_if(is.numeric, scale) %>% 
  pivot_longer(cols = 1:ifelse(x==4,2,4)+seq(0,16,by = 4)[x],
               names_to = "Variabili",
               values_to = "Valori") %>%
  group_by(fraud) %>% 
  ggplot(aes(y = Valori, x = fraud, color = fraud, fill = fraud)) +
  geom_boxplot(alpha = 0.6) + 
  scale_fill_viridis(discrete=TRUE) +
    scale_color_viridis(discrete=TRUE) +
    theme_ipsum() +
    theme(
      legend.position="none",
      panel.spacing = unit(0.1, "lines"),
      strip.text.x = element_text(size = 8)
    ) +
    xlab("") +
    ylab("") +
    facet_wrap(~Variabili, ncol = 2, scales = "free_y")
}
)

for(i in 1:4) print(grafici[[i]])
rm(grafici)
```



## Correlazione


```{r}
# df[c(1:1000, 50000:51000),] %>% 
df %>% 
  mutate(fraud = as.numeric(fraud)) %>% 
  select(where(is.numeric)) %>% 
  # select(where(~ length(levels(factor(.))) > 2)) %>% 
  rename_all(~ str_trunc(., width = 15)) %>% 
  scale() %>% 
  cor() %>% 
  corrplot::corrplot(cl.cex = 0.5,
                     type = "lower",
                    tl.pos = "l")
```
Il seguente grafico mostra in modo diretto la correlazione della variabile fraud con le altre variabili del modello. 
Ovviamente, non è rappresentata la multicollinearità tra le altre variabili, come mostrato nel grafico soprastante.

```{r}
df %>% 
  mutate(fraud = as.numeric(fraud)) %>% 
  select(where(is.numeric)) %>% 
  # select(where(~ length(levels(factor(.))) > 2)) %>% 
  rename_all(~ str_trunc(., width = 15)) %>% 
  scale() %>% 
  cor() %>% 
  data.frame() %>%
  rownames_to_column("Variabili") %>% 
  tibble() %>% 
  select(Variabili,fraud) %>% 
  mutate(Variabili = as.factor(Variabili)) %>% 
  filter(fraud<1) %>% 
  ggplot(aes(x = Variabili, y = fraud, fill = fraud)) +
  geom_col(color = "gray45") +
  coord_flip() +
  scale_fill_gradient2(low = "mediumspringgreen",
                       high = "mediumvioletred",
                       midpoint = 0) + 
  theme_bw() +
  theme(legend.position = "none")

```


# Modelli

## Divisione dei dati

```{r}
data_split = initial_split(df, prop = 3/4, strata = "fraud")

data_train = training(data_split)
data_test = testing(data_split)
```


## Preprocessamento

`step_novel(all_nominal_predictors())`: Identifica variabili categoriche che potrebbero contenere nuove categorie non presenti nei dati di addestramento. Questo passo è utile per gestire categorie non viste durante l'addestramento del modello che invece potrebbero essere presenti nel dataset di testing del modello.

`step_normalize(all_numeric_predictors())`: Normalizza tutte le variabili numeriche in modo che abbiano una media di 0 e una deviazione standard di 1. Questo può essere utile per garantire che tutte le variabili abbiano lo stesso peso nella costruzione del modello. 
Inoltre, sarà comodo per capire il peso dei coefficienti nel modello.

`step_dummy(all_nominal_predictors())`: Crea variabili fittizie (dummy variables) per le variabili categoriche, consentendo al modello di trattare correttamente le categorie senza un ordine intrinseco.

`step_rose(fraud)`: Utilizza il metodo di sovracampionamento chiamato ROSE (Random Over-Sampling Examples) per bilanciare le classi nel dataset. Infatti i casi di frode sono nettamente inferiori (1,1%) ai casi in cui non si verifica la frode.

`step_zv(all_nominal_predictors())`: Rimuove le variabili categoriche che hanno varianza zero, il che significa che sono costanti e non portano alcuna informazione predittiva.

```{r, eval=FALSE}
df_rec = 
  recipe(fraud ~ ., data = data_train) %>% 
  step_novel(all_nominal_predictors()) %>% 
  step_normalize(all_numeric_predictors()) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  step_rose(fraud) %>% 
  step_zv(all_nominal_predictors())
```

### Campionamento e sottocampionamento

```{r, eval=FALSE}
df %>% 
  select(fraud) %>% 
  table() %>%
  as_tibble() %>% 
  mutate(perc = round(n/nrow(df)*100, digits = 1))
```

## Regressione logistica

### Specificazione modello

```{r, eval=FALSE}
log_mod = logistic_reg() %>% 
  set_engine("glm") %>% 
  set_mode("classification")

log_wflw = 
  workflow() %>% 
  add_model(log_mod) %>% 
  add_recipe(df_rec)

log_fit = fit(log_wflw, data_train)
```


```{r, echo=FALSE}
log_fit = readRDS("var/log_fit")

lev.p.value = c("very low confidence", "low confidence", "confidence", "high confidence", "very high confidence")
log_coef = log_fit %>% 
  tidy() %>%
  mutate(OR = exp(estimate),
         OR.conf.low = exp(estimate - 1.96*std.error),
         OR.conf.high = exp(estimate + 1.96*std.error),
         p.value.level = case_when(p.value > .1 ~ lev.p.value[1],
                                   p.value > .05 ~ lev.p.value[2],
                                   p.value > .01 ~ lev.p.value[3],
                                   p.value > .001 ~ lev.p.value[4],
                                   .default = lev.p.value[5])) %>% 
  mutate(p.value.level = factor(p.value.level, levels = lev.p.value, ordered = T))

logistic_output = data_test %>%
  bind_cols(predict(log_fit, new_data = data_test, type = 'prob')) %>%
  mutate(.pred_class = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = .5))

rm(log_fit)
gc()
```


```{r}
log_coef
```



```{r}
creaGruppoVariabili = function(data, varOrig, colVar)
{
  macroGruppi = matrix(nrow = nrow(data))
  for(i in 1:nrow(data))
  {
    if(data[i,colVar] == "(Intercept)")
    {
      macroGruppi[i,1] = "Intercetta"
    }else
    {
      if(sum(str_equal(data[i,colVar], varOrig))==1)
      {
        macroGruppi[i,1] = data[i,colVar] %>% as.character()
      }else
      {
        data[i,colVar] = data[i,colVar] %>% 
          str_sub(end = data[i,colVar] %>% 
                    str_locate_all(pattern = "_") %>% 
                    unlist() %>% 
                    tail(1) - 1)
        macroGruppi[i,1] = data[i,colVar] %>% as.character()
      }
    }
  }
  return(macroGruppi)
}

grafico = log_coef %>%
  mutate(Variabili = creaGruppoVariabili(log_coef, names(df), "term")) %>% 
  ggplot(aes(x = Variabili, y = OR, ymin = OR.conf.low,
             ymax = OR.conf.high, color = p.value, label = term)) +
  geom_pointrange(position = "jitter") +
  coord_flip() +
  scale_x_discrete(labels = str_sub(names(df),end = 10)) +
  scale_color_gradient2(low = "mediumspringgreen",
                        mid = "mediumvioletred",
                       high = "mediumvioletred",
                       midpoint = 0.4) +
  # scale_color_brewer(palette = palette(palette.colors(palette = "RdYlGn"))) +
  # scale_color_gradient2(low = "springgreen",
  #                       high = "tomato",
  #                       mid = "orange",
  #                       midpoint = .05, 
  #                       na.value = "gray") +
  xlab("Coefficienti") +
  geom_hline(aes(yintercept = 1), color = "tomato") +
  ylim(c(min(log_coef$OR.conf.low),2*median(log_coef$OR.conf.high, na.rm = T))) +
  theme(legend.position = "none", 
        legend.box = "horizontal",
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 7),
        axis.text.x = element_text(vjust = 1, size = 8)) 

ggplotly(grafico)
```

### Valutazione del modello

Dalle probabilità creo la predizione della classe in base al threshold impostato, per ora pongo 0.5. successivamente lo andrò ad adattare per trovare la soluzione migliore al problema.



```{r}
logistic_output %>%
  ggplot(aes(x = fraud, y = `.pred_No fraud`)) +
  geom_violin(fill = "mediumspringgreen", alpha = .6) + 
  geom_hline(yintercept = 0.5, color='mediumvioletred') +  # try changing threshold
  geom_hline(yintercept = 0.73, color='mediumvioletred') +  # try changing threshold
  labs(y = 'Predicted Probability of Outcome', x = 'Observed Outcome') +
  scale_y_percent() +
  theme_classic()
```

Osservo le metriche delle predizioni del modello con il threshold impostato al 0.5.
Le metriche scelte in questo caso sono la sensitività, specificità, accuratezza e precione. 


```{r}
log_metrics = metric_set(sens, yardstick::spec, accuracy,yardstick::precision) 

logistic_output %>% 
  log_metrics(estimate = .pred_class, truth = fraud, event_level = "second") 
```
#### Matrice di confusione

<!-- Attraverso la curva ROC vado a cercare la migliore soluzione tra FPR e TPR (recall) in modo tale da trovare il threshold che massimizza la specificità con la sensitività al più del 95%. -->

Il threshold per trovare un recall del 5% è stato impostato a 0,73.

```{r}
log_roc = logistic_output %>% 
    roc_curve(fraud, .pred_Fraud, event_level = "second")  # con second dico che la seconda classe (fraud) è la classe che è interesse di studio
  

# log_bestThreshold = log_roc %>% 
#   filter(sensitivity > 0.5) %>% 
#   filter(specificity == max(specificity))

logistic_output %>%
  # mutate(pred_best = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = log_bestThreshold$.threshold)) %>% 
  mutate(pred_best = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = 0.73)) %>% 
  select(fraud,pred_best) %>% 
  table() %>% 
  caret::confusionMatrix(positive = "Fraud")
```


#### Curva ROC

```{r}
log_roc %>% 
  mutate(FPR = 1 - specificity,
         TPR = sensitivity) %>% 
  ggplot(aes(x = FPR, y = TPR)) +
  geom_line() + 
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_vline(aes(xintercept = log_roc %>% 
                   filter(.threshold > .73) %>%
                   filter(.threshold == min(.threshold)) %>% 
                   select(specificity) %>% 
                   mutate(specificity = 1 - specificity) %>% 
                   as.numeric()),
             color = "tomato",
             linetype = "dashed") +
  geom_hline(aes(yintercept = log_roc %>% 
                   filter(.threshold > .73) %>%
                   filter(.threshold == min(.threshold)) %>% 
                   select(sensitivity) %>% 
                   as.numeric()),
             color = "tomato",
             linetype = "dashed") +
  labs(title = "Curva ROC") +
  theme_bw() +
  theme(aspect.ratio = 1)
```


# Dataframe ridotto

Per problemi computazionali, dovuti alla grande quantità di dati e di features, si è deciso di ridurre la dimensione del campione riducendo del 95% le osservazioni.
L'obiettivo del lavoro è quello di creare un modello che identifichi le frodi bancarie, perciò nella riduzione del dataframe è stato scelto di conservare tutte le osservazioni nelle quali è stata commessa una frode.

```{r}
# dfRid = ovun.sample(fraud ~.,
#                  data = df,
#                  method = "under",
#                  p = 0.5,
#                  seed = 1)$data %>% tibble
dfRid = df %>% 
  filter(fraud == "Fraud") %>% 
  add_row(df %>% 
                filter(fraud != "Fraud") %>% 
                slice_sample(n = nrow(df)*0.05 - nrow(filter(df, fraud == "Fraud"))))

dfRid
```

L'attuale dataframe contiene ancora le classi sbilanciate, si procederà quindi al ribilanciamento attraverso el diverse tecniche di ricampionamento

```{r}
table(dfRid$fraud)
```

## Variabile fraud

```{r, warning=FALSE}
apply(matrix(1:(ncol(dfChar) - 1)),1, function(indexCol)
# apply(matrix(1:2),1, function(indexCol)
{
dfRid %>% 
  select(where(is.factor)) %>%  
  ggplot(aes(x = get(colnames(dfChar)[indexCol]), group = fraud)) +
  geom_bar(aes(y = ..prop.., fill =  fraud), stat = "count") +
  scale_y_continuous(labels=scales::percent) +
  xlab(colnames(dfChar)[indexCol]) +
  ylab("Osservazioni") +
  theme_ipsum() +
  theme(legend.position="none")
}
)
```


## Correlazione


```{r}
dfRid %>% 
  mutate(fraud = as.numeric(fraud)) %>% 
  select(where(is.numeric)) %>% 
  rename_all(~ str_trunc(., width = 15)) %>% 
  scale() %>% 
  cor() %>% 
  corrplot::corrplot(cl.cex = 0.5,
                     type = "lower",
                    tl.pos = "l")
```



```{r}
dfRid %>% 
  mutate(fraud = as.numeric(fraud)) %>% 
  select(where(is.numeric)) %>% 
  # select(where(~ length(levels(factor(.))) > 2)) %>% 
  rename_all(~ str_trunc(., width = 15)) %>% 
  scale() %>% 
  cor() %>% 
  data.frame() %>%
  rownames_to_column("Variabili") %>% 
  tibble() %>% 
  select(Variabili,fraud) %>% 
  mutate(Variabili = as.factor(Variabili)) %>% 
  filter(fraud<1) %>% 
  ggplot(aes(x = Variabili, y = fraud, fill = fraud)) +
  geom_col(color = "gray45") +
  coord_flip() +
  scale_fill_gradient2(low = "mediumspringgreen",
                       high = "mediumvioletred",
                       midpoint = 0) + 
  theme_bw() +
  theme(legend.position = "none")

```

## Preprocessamento

```{r}
data_split = initial_split(dfRid, prop = 3/4, strata = "fraud")

data_train = training(data_split)
data_test = testing(data_split)
```


```{r}
df_rec = 
  recipe(fraud ~ ., data = data_train) %>% 
  step_novel(all_nominal_predictors()) %>% 
  step_normalize(all_numeric_predictors()) %>% 
  step_dummy(all_nominal_predictors()) %>% 
  step_zv(all_nominal_predictors())
```

## Regressione logistica

### Specificazione modello

```{r}
log_mod = logistic_reg() %>% 
  set_engine("glm") %>% 
  set_mode("classification")
```

Imposto il workflow, aggiungendo il modello e il preprocessamento dei dati

```{r}
log_wflw = 
  workflow() %>% 
  add_model(log_mod) %>% 
  add_recipe(df_rec)
```

### Performance modello

```{r}
log_fit_rid = fit(log_wflw, data_train)

lev.p.value = c("very low confidence", "low confidence", "confidence", "high confidence", "very high confidence")
log_coef_rid = log_fit_rid %>% 
  tidy() %>%
  mutate(OR = exp(estimate),
         OR.conf.low = exp(estimate - 1.96*std.error),
         OR.conf.high = exp(estimate + 1.96*std.error),
         p.value.level = case_when(p.value > .1 ~ lev.p.value[1],
                                   p.value > .05 ~ lev.p.value[2],
                                   p.value > .01 ~ lev.p.value[3],
                                   p.value > .001 ~ lev.p.value[4],
                                   .default = lev.p.value[5])) %>% 
  mutate(p.value.level = factor(p.value.level, levels = lev.p.value, ordered = T))
log_coef_rid
```


```{r}
creaGruppoVariabili = function(data, varOrig, colVar)
{
  macroGruppi = matrix(nrow = nrow(data))
  for(i in 1:nrow(data))
  {
    if(data[i,colVar] == "(Intercept)")
    {
      macroGruppi[i,1] = "Intercetta"
    }else
    {
      if(sum(str_equal(data[i,colVar], varOrig))==1)
      {
        macroGruppi[i,1] = data[i,colVar] %>% as.character()
      }else
      {
        data[i,colVar] = data[i,colVar] %>% 
          str_sub(end = data[i,colVar] %>% 
                    str_locate_all(pattern = "_") %>% 
                    unlist() %>% 
                    tail(1) - 1)
        macroGruppi[i,1] = data[i,colVar] %>% as.character()
      }
    }
  }
  return(macroGruppi)
}

grafico = log_coef_rid %>%
  mutate(Variabili = creaGruppoVariabili(log_coef_rid, names(dfRid), "term")) %>% 
  ggplot(aes(x = Variabili, y = OR, ymin = OR.conf.low,
             ymax = OR.conf.high, color = p.value.level, label = term)) +
  geom_pointrange(position = "jitter") +
  coord_flip() +
  scale_x_discrete(labels = str_sub(names(dfRid),end = 10)) +
  scale_color_brewer(palette = "RdYlGn") +
  # scale_color_gradient2(low = "springgreen",
  #                       high = "tomato",
  #                       mid = "orange",
  #                       midpoint = .05, 
  #                       na.value = "gray") +
  xlab("Coefficienti") +
  geom_hline(aes(yintercept = 1), color = "tomato") +
  ylim(c(min(log_coef_rid$OR.conf.low),2*median(log_coef_rid$OR.conf.high, na.rm = T))) +
  theme(legend.position = "top", 
        legend.box = "horizontal",
        legend.title = element_text(size = 9),
        legend.text = element_text(size = 7),
        axis.text.x = element_text(vjust = 1, size = 8)) 

# grafico

ggplotly(grafico)
```

### Valutazione del modello

Dalle probabilità creo la predizione della classe in base al threshold impostato, per ora pongo 0.5. successivamente lo andrò ad adattare per trovare la soluzione migliore al problema.

```{r}
log_output_rid <- data_test %>%
  bind_cols(predict(log_fit_rid, new_data = data_test, type = 'prob')) %>%
  mutate(.pred_class = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = .5))

```



```{r}
log_output_rid %>%
  ggplot(aes(x = fraud, y = `.pred_No fraud`)) +
  geom_violin(fill = "mediumspringgreen", alpha = .6) + 
  geom_hline(yintercept = 0.5, color='mediumvioletred') +  # try changing threshold
  geom_hline(yintercept = 0.6, color='mediumvioletred') +  # try changing threshold
  labs(y = 'Predicted Probability of Outcome', x = 'Observed Outcome') +
  scale_y_percent() +
  theme_classic()
```

Osservo le metriche delle predizioni del modello con il threshold impostato al 0.5.
Le metriche scelte in questo caso sono la sensitività, specificità, accuratezza e precione. 


```{r}
log_metrics = metric_set(sens, yardstick::spec, accuracy,yardstick::precision) 

logistic_output %>% 
  log_metrics(estimate = .pred_class, truth = fraud, event_level = "second") 
```

#### Matrice di confusione

<!-- Attraverso la curva ROC vado a cercare la migliore soluzione tra FPR e TPR (recall) in modo tale da trovare il threshold che massimizza la specificità con la sensitività al più del 95%. -->

Il threshold per trovare un recall del 5% è stato impostato a 0,73.

```{r}
log_roc = logistic_output %>% 
    roc_curve(fraud, .pred_Fraud, event_level = "second")  # con second dico che la seconda classe (fraud) è la classe che è interesse di studio
  

# log_bestThreshold = log_roc %>% 
#   filter(sensitivity > 0.5) %>% 
#   filter(specificity == max(specificity))

logistic_output %>%
  # mutate(pred_best = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = log_bestThreshold$.threshold)) %>% 
  mutate(pred_best = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = 0.6)) %>% 
  select(fraud,pred_best) %>% 
  table() %>% 
  caret::confusionMatrix(positive = "Fraud")
```


#### Curva ROC

```{r}
log_roc %>% 
  mutate(FPR = 1 - specificity,
         TPR = sensitivity) %>% 
  ggplot(aes(x = FPR, y = TPR)) +
  geom_line() + 
  geom_abline(aes(intercept = 0, slope = 1)) +
  geom_vline(aes(xintercept = log_roc %>% 
                   filter(.threshold > .73) %>%
                   filter(.threshold == min(.threshold)) %>% 
                   select(specificity) %>% 
                   mutate(specificity = 1 - specificity) %>% 
                   as.numeric()),
             color = "tomato",
             linetype = "dashed") +
  geom_hline(aes(yintercept = log_roc %>% 
                   filter(.threshold > .73) %>%
                   filter(.threshold == min(.threshold)) %>% 
                   select(sensitivity) %>% 
                   as.numeric()),
             color = "tomato",
             linetype = "dashed") +
  labs(title = "Curva ROC") +
  theme_bw() +
  theme(aspect.ratio = 1)
```

### Ricampionamento

```{r}
set.seed(1)
cv_folds = vfold_cv(data_train, strata = "fraud",repeats = 10)

cv_folds
```


```{r, eval=FALSE}
log_metrics = metric_set(sens, yardstick::spec, accuracy)

set.seed(1)
log_res = fit_resamples(
  log_wflw, 
  control = control_resamples(save_pred = TRUE, event_level = 'second'),
  resamples = cv_folds,
  metrics = log_metrics
)

saveRDS(log_res = "var/log_res.RData")
```

```{r, eval=TRUE}
log_res = readRDS("var/log_res.RData")

log_res_metrics = lapply(1:nrow(log_res), function(i) {
  log_res %>%
    select(.metrics) %>%
    dplyr::slice(i) %>%
    unlist() %>%
    t() %>%
    data.frame() %>%
    select_if(str_detect(names(.), ".metrics..estimate"))
}) %>% 
  purrr::map_dfr( ~ as.data.frame(.)) %>% 
  transmute(sensitivity = .metrics..estimate1,
         specificity = .metrics..estimate2,
         accuracy = .metrics..estimate3)

log_res_metrics
```



```{r}
log_res_metrics %>% 
  mutate_all(as.numeric) %>% 
  # add_row(sensitivity = 0,specificity = 0,accuracy = 0, .before = 1) %>% 
  # add_row(sensitivity = 1,specificity = 1,accuracy = 1, .after =  100) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity)) +
  geom_point(aes(color = accuracy)) +
  scale_color_gradient2(high = "springgreen",
                        low = "tomato",
                        mid = "orange",
                        midpoint = mean(as.numeric(log_res_metrics$accuracy)),
                        na.value = "gray") +
  geom_vline(aes(xintercept = log_res_metrics %>%
                   mutate_all(as.numeric) %>%
                   filter(specificity > .95) %>%
                   filter(sensitivity == max(sensitivity)) %>% 
                   select(specificity) %>% 
                   mutate(specificity = 1 - specificity) %>% 
                   as.numeric()),
             color = "tomato",
             linetype = "dashed") +
  geom_hline(aes(yintercept = log_res_metrics %>% 
                   mutate_all(as.numeric) %>%
                   filter(specificity > .95) %>%
                   filter(sensitivity == max(sensitivity)) %>% 
                   select(sensitivity) %>% 
                   as.numeric()),
             color = "tomato",
             linetype = "dashed") +
  theme_bw() 
```



```{r}
metriche = log_res %>% 
  select(.metrics) %>% 
  map_dfr(as.data.frame) %>% 
  select(.metric, .estimate, .config) %>%  
  pivot_wider(names_from = .metric, values_from = .estimate, id_cols = .config) 

metriche
```

Il modello 


```{r}
log_res %>% 
  collect_predictions() %>% 
  select(fraud, .pred_class) %>% 
  table() %>% 
  caret::confusionMatrix(positive = "Fraud")
```



## Random Forest


### Specificazione modello

```{r}
rf_mod = rand_forest() %>%
  set_engine(engine = 'ranger') %>% 
  set_args(mtry = NULL, # lo vado a inserire dopo 
           trees = 1000,
           min_n = 2,
           probability = TRUE, 
           importance = 'impurity') %>% 
  set_mode('classification') 
```

### Workflows

```{r}
rf_wflw_mtry2 = workflow() %>%
  add_model(rf_mod %>% set_args(mtry = 2)) %>%
  add_recipe(df_rec)

rf_wflw_mtry5 = workflow() %>%
  add_model(rf_mod %>% set_args(mtry = floor(sqrt(ncol(dfRid))))) %>% # 5
  add_recipe(df_rec)

rf_wflw_mtry12 = workflow() %>%
  add_model(rf_mod %>% set_args(mtry = 12)) %>%
  add_recipe(df_rec)
```

### Fit dei modelli

```{r, eval=FALSE}
set.seed(1) 
rf_fit_mtry2 = fit(rf_wflw_mtry2, data = data_train)

set.seed(1)
rf_fit_mtry12 = fit(rf_wflw_mtry12, data = data_train)

set.seed(1)
rf_fit_mtry5 = fit(rf_wflw_mtry5, data = data_train)

saveRDS(rf_fit_mtry2, "var/rf_fit_mtry2.RData")
saveRDS(rf_fit_mtry5, "var/rf_fit_mtry5.RData")
saveRDS(rf_fit_mtry12, "var/rf_fit_mtry12.RData")
```

### Valutazione del modello

```{r}
rf_fit_mtry2 = readRDS("var/rf_fit_mtry2.RData") # importo il modello nel'enviroment
rf_fit_mtry5 = readRDS("var/rf_fit_mtry5.RData") # importo il modello nel'enviroment
rf_fit_mtry12 = readRDS("var/rf_fit_mtry12.RData") # importo il modello nel'enviroment


rf_out_mtry2 = rf_fit_mtry2 %>% 
  extract_fit_engine()

rf_pred_mtry2 = predict(rf_fit_mtry2,
                    new_data = data_test,
                    type = 'prob')

rf_out_mtry5 = rf_fit_mtry5 %>% 
  extract_fit_engine()

rf_pred_mtry5 = predict(rf_fit_mtry5,
                    new_data = data_test,
                    type = 'prob')

rf_out_mtry12 = rf_fit_mtry12 %>% 
  extract_fit_engine()

rf_pred_mtry12 = predict(rf_fit_mtry12,
                    new_data = data_test,
                    type = 'prob')
```


```{r}
data_test %>% 
  select(fraud) %>% 
  add_column(rf_pred_mtry2, .name_repair = "unique") %>%
  add_column(rf_pred_mtry5, .name_repair = "unique") %>%
  add_column(rf_pred_mtry12, .name_repair = "unique") %>%
  pivot_longer(cols = c(2,4,6), values_to = "value", names_to = "mtry") %>% 
  mutate(mtry = factor(mtry, labels = c("mtry2", "mtry5", "mtry12"))) %>% 
  ggplot(aes(x = fraud, y = value)) +
  geom_violin(fill = "mediumspringgreen", alpha = .6) + 
  geom_hline(yintercept = 0.5, color='mediumvioletred') +  # try changing threshold
  geom_hline(yintercept = 0.85, color='mediumvioletred') +  # try changing threshold
  facet_grid(cols = dplyr::vars(mtry)) +
  labs(y = 'Predicted Probability of Outcome', x = 'Observed Outcome') +
  scale_y_percent() +
  theme_classic()

```
#### Matrice di confusione

```{r}

data_test %>% 
  select(fraud) %>% 
  add_column(rf_pred_mtry12 %>% 
               transmute(predicted = factor(`.pred_No fraud`<.85,
                                            labels = c("No fraud", "Fraud")))) %>% 
  table() %>% 
  caret::confusionMatrix(positive = "Fraud")
```

### Selezione variabili 

Nel grafico si può vedere come mtry12 e mtry5 hanno valori sull'importanza delle variabili decisamente superiori di mtry2, inoltre le prime variabili hanno una classifica simile tra mtry12 e mtry5, per mtry2, invece, le variabili più importanti sarebbero altre. 

procediamo quindi ha eseguire una selezione delle variabili per snellire il dataset. Verrannno selezionate solo alcune modalità di variabili categoriche, avendo già il dataset diviso in dummies non ci saranno problemi nella creazione di nuove categorie.

```{r}
rf_varImportance = rf_out_mtry2 %>% 
  vip::vi() %>% 
  right_join(rf_out_mtry5 %>% vip::vi(), by = "Variable") %>% 
  right_join(rf_out_mtry12 %>% vip::vi(), by = "Variable") 
  
colnames(rf_varImportance) = c("Variable", "mtry2", "mtry5", "mtry12")

rf_varImportance %>% 
  # mutate(Macrovariabili = creaGruppoVariabili(rf_varImportance, names(dfRid), "Variable")) %>% 
  head(25) %>% 
  pivot_longer(cols = 2:4,
               names_to = "mtry",
               values_to = "Importance") %>%
  mutate(mtry = mtry %>% 
           factor(levels = c("mtry2", "mtry5", "mtry12"), ordered = T)) %>% 
  mutate(Variable = Variable %>% 
           factor() %>% 
           reorder(Importance)) %>% 
  # mutate(mtry = mtry %>% 
  #          factor() %>% 
  #          reorder(c("mtry2", "mtry5", "mtry12"))) %>% 
  ggplot(aes(x = Variable, y = Importance)) +
  geom_col(aes(fill = mtry), position = "dodge") + 
  scale_fill_brewer(palette = "Dark2") +
  coord_flip()
```

Seleziono le 25 variabili più importanti per il modello random forest con il parametro mtry12.

```{r}
var_Importance = rf_varImportance %>% 
  dplyr::arrange(across(mtry12)) %>% 
  tail(25) %>% 
  select(Variable)

df_rec_25var = df_rec %>% 
  step_select(c(var_Importance$Variable, "fraud")) %>% 
  step_nzv(all_factor_predictors()) %>%
  step_corr(all_predictors())  

df_rec_25var

data_train_25var = df_rec_25var %>% 
  prep(training = data_train) %>% 
  bake(new_data = NULL)

data_test_25var = df_rec_25var %>% 
  prep(training = data_test)  %>% 
  bake(new_data = NULL)

rm(rf_fit_mtry2) # elimino il modello per liberare spazio
rm(rf_fit_mtry5) # elimino il modello per liberare spazio
rm(rf_fit_mtry12) # elimino il modello per liberare spazio
gc()

```


### Tuning

Svolgo il lavoro su differenti core della CPU.

```{r}
cores = parallel::detectCores() - 1
cores
```

```{r}
rf_mod_tun =
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) %>% 
  set_engine("ranger", num.threads = cores) %>% 
  set_args(probability = TRUE, 
           importance = 'impurity') %>% 
  set_mode("classification")

rf_workflow =
  workflow() %>% 
  add_model(rf_mod_tun) %>% 
  add_recipe(df_rec_25var)

rf_workflow
```

```{r, eval=FALSE}
set.seed(1)
rf_res =
  rf_workflow %>% 
  tune_grid(data_folds,
            grid = 25,
            control = control_grid(save_pred = TRUE,
                                   event_level = "second"),
            metrics = metric_set(roc_auc, 
                                 yardstick::specificity,
                                 yardstick::sensitivity))
saveRDS(rf_res, "var/rf_res.RData")
```
Alleno 250 modelli totali utilizzando la cross validation (k-fold) con il train set   diviso in 10 parti uguali. Per ogni modello k-1 gruppi del train set verranno utilizzati per allenare il modello, mentre il restante sarà usato per testare il modello.


```{r}
rf_res = readRDS("var/rf_res.RData")
metriche = rf_res %>% 
  select(.metrics) %>% 
  map_dfr(as.data.frame) %>% 
  select(mtry, min_n, .metric, .estimate, .config) 

metriche_pivot = metriche %>% 
  pivot_wider(names_from = .metric, values_from = .estimate, id_cols = .config) %>% 
  left_join(metriche %>% 
              filter(.metric == "specificity") %>% 
              select(mtry, min_n, .config),
            by = ".config")

metriche_pivot
```



```{r}
autoplot(rf_res) +
  geom_vline(aes(xintercept = 5), color = "tomato") +
  theme_bw()
```


```{r}
rf_best = metriche_pivot %>% 
  filter(specificity > .95) %>% 
  filter(sensitivity == max(sensitivity))

rf_best
```


```{r}
rf_auc =
  rf_res %>% 
  collect_predictions(parameters = rf_best) %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  mutate(model = "Random Forest")

rf_auc %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity)) + 
  geom_path(lwd = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() 
```


## Albero di decisione

### Tuning

```{r}
tree_mod = 
  decision_tree(
    cost_complexity = tune(),
    tree_depth = tune()
  ) %>% 
  set_engine("rpart") %>% 
  set_mode("classification")


tree_grid = grid_regular(cost_complexity(),
                          tree_depth(),
                          levels = 5)


tree_wflw = workflow() %>%
  add_recipe(df_rec_25var) %>%
  # add_recipe(recipe(fraud ~ ., data = data_train_25var)) %>%   
  add_model(tree_mod) 

set.seed(1)
data_folds = vfold_cv(data_train)

tree_wflw
```


```{r, eval=FALSE}
tree_res = 
  tree_wflw %>% 
  tune_grid(
    resamples = data_folds,
    grid = tree_grid
    )
saveRDS(tree_res, "var/tree_res.RData")
```



```{r}
tree_res = readRDS("var/tree_res.RData")

tree_res %>%
  collect_metrics() %>%
  mutate(tree_depth = factor(tree_depth)) %>%
  ggplot(aes(cost_complexity, mean, color = tree_depth)) +
  geom_line(linewidth = 1.5, alpha = 0.6) +
  geom_point(size = 2) +
  facet_wrap(~ .metric, scales = "free", nrow = 2) +
  scale_x_log10(labels = scales::label_number()) +
  scale_color_viridis_d(option = "plasma", begin = .9, end = 0)

```


```{r}
tree_res %>%
  show_best("accuracy")
```




```{r}
best_tree = tree_res %>%
  select_best("accuracy")

best_tree
```



```{r}
tree_final_wflw = 
  tree_wflw %>% 
  finalize_workflow(best_tree)
```

```{r}
tree_fit = 
  tree_final_wflw %>%
  last_fit(data_split) 

tree_fit %>%
  collect_metrics()
```


```{r}
tree_fit %>%
  collect_predictions() %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  autoplot()
```

```{r}
tree_fit %>%
  extract_fit_engine() %>%
  rpart.plot::rpart.plot()
  # plot()
# text()
```

```{r}
tree_fit %>%
  collect_predictions() %>% 
  mutate(pred_best = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = 0.75)) %>% 
  select(fraud,pred_best) %>% 
  table() %>% 
  caret::confusionMatrix(positive = "Fraud")
```


## GLM 


```{r}
lr_mod =
  logistic_reg(penalty = tune(), mixture = 1) %>% 
  set_engine("glmnet") %>% 
  set_mode("classification") %>% 
  set_args(family = "binomial")

lr_workflow =
  workflow() %>% 
  add_model(lr_mod) %>% 
  add_recipe(df_rec_25var)

lr_reg_grid = tibble(penalty = 10^ seq(-4, -1, length.out = 30))

lr_workflow
```


```{r, eval=FALSE}
lr_res =
  lr_workflow %>% 
  tune_grid(data_folds,
            grid = lr_reg_grid,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(roc_auc))
saveRDS(lr_res, "var/lr_res.RData")
```


```{r}
lr_res = readRDS("var/lr_res.RData")
  lr_res %>% 
  collect_metrics() %>% 
  ggplot(aes(x = penalty, y = mean)) + 
  geom_point() + 
  geom_line() + 
  ylab("Area under the ROC Curve") +
  geom_vline(aes(xintercept =   lr_res %>% 
                   collect_metrics() %>% 
                   select(penalty) %>% 
                   slice(20) %>% 
                   as.numeric()),
             color = "tomato") +
  scale_x_log10(labels = scales::label_number())
```



```{r}
top_models =
  lr_res %>% 
  show_best("roc_auc", n = 20) %>% 
  arrange(penalty) 
top_models
```

Scelgo il 20esimo modello migliore considerando l'AUC per massimizzare la penalità e di conseguenza migliorare il modello.

```{r}
lr_best =
  lr_res %>% 
  collect_metrics() %>% 
  arrange(penalty) %>% 
  slice(20)
lr_best
```

```{r}
lr_auc =  lr_res %>% 
  collect_predictions(parameters = lr_best) %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  mutate(model = "Regressione logistica")

lr_auc %>% autoplot(lr_auc)
```

```{r}
lr_res %>% 
  collect_predictions(parameters = lr_best) %>% 
  mutate(pred_best = probably::make_two_class_pred(`.pred_No fraud`, levels(fraud), threshold = 0.680)) %>% 
  select(fraud,pred_best) %>% 
  table() %>% 
  caret::confusionMatrix(positive = "Fraud")
```

```{r}
lr_final_wflw = 
  lr_workflow %>% 
  finalize_workflow(lr_best)
```

```{r}
lr_fit = 
  lr_final_wflw %>%
  last_fit(data_split)

```




# Modello finale

In questa fase si procederà ad allenare i modelli sul dataset originale con i parametri trovati nei capitoli precedenti. 
Ci si ritroverà in una situazione di imbilanciamento delle classe, sarà quindi opportuno utilizzare una tecnica di ribilanciamente. 
In questo caso è stata utilizzata ROSE. 

```{r, eval=FALSE}
set.seed(1)
df_split = initial_split(df, prop = 3/4, strata = "fraud")

df_train = training(data_split)
df_test = testing(data_split)
```


## Regressione logistica

```{r, eval=FALSE}
Sys.time()
lr_fit_finale = lr_final_wflw %>% last_fit(df_split)
Sys.time()

saveRDS(lr_fit_finale, "var/lr_fit_finale.RData")
```

### Curva ROC 

```{r}
lr_fit_finale = readRDS("var/lr_fit_finale.RData")

lr_auc_finale = lr_fit_finale %>% 
  collect_predictions() %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  mutate(model = "Regressione logistica finale senza ROSE")
```


```{r}
bind_rows(rf_auc, lr_auc, lr_auc_finale) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6)
```

## Regressione logistica con ROSE

```{r, eval=FALSE}
lr_final_wflw_rose = lr_final_wflw %>%
  update_recipe(df_rec %>% step_rose())

Sys.time()
lr_fit_finale_rose = lr_final_wflw_rose %>% last_fit(df_split)
Sys.time()

saveRDS(lr_fit_finale_rose, "var/lr_fit_finale_rose.RData")
```

### Curva ROC 

```{r, eval = FALSE}
lr_fit_finale_rose = readRDS("var/lr_fit_finale_rose.RData")
lr_auc_finale_rose = lr_fit_finale_rose %>% 
  collect_predictions() %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  mutate(model = "Regressione logistica finale con ROSE")
```


```{r}
bind_rows(rf_auc, lr_auc, lr_auc_finale, lr_auc_finale_rose) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1, alpha = 0.6) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = 1)
```

## Random forest senza ROSE

```{r}
rf_wflw_finale = rf_workflow %>% 
  finalize_workflow(rf_best)

rf_wflw_finale
```

```{r, eval=FALSE}
Sys.time()
rf_fit_finale = rf_wflw_finale %>% last_fit(df_split)
Sys.time()

saveRDS(rf_fit_finale, "var/rf_fit_finale.RData")
```


### Curva ROC

```{r, eval=FALSE}
rf_fit_finale = readRDS("var/rf_fit_finale.RData")

rf_auc_finale = rf_fit_finale %>% 
  collect_predictions() %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  mutate(model = "Random forest finale")
```

```{r}
bind_rows(rf_auc, lr_auc, lr_auc_finale, rf_auc_finale) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6)
```


### Importanza delle variabili

```{r}
rf_fit_finale %>%
  extract_fit_engine() %>% 
  vip::vi() %>% 
  head(25) %>% 
  arrange(Importance) %>% 
  ggplot(aes(x = Variable, y = Importance)) +
  geom_col() +
  coord_flip()
```

## Random forest con ROSE

```{r, eval=FALSE}
rf_wflw_finale_rose = rf_wflw_finale %>%
  update_recipe(df_rec %>% step_rose)
  

Sys.time()
rf_fit_finale_rose = rf_wflw_finale_rose %>% last_fit(df_split)
Sys.time()

saveRDS(rf_fit_finale_rose, "var/rf_fit_finale_rose.RData")
```


```{r, eval=FALSE}
rf_fit_finale_rose = readRDS("var/rf_fit_finale_rose.RData")
rf_auc_finale_rose = rf_fit_finale_rose %>% 
  collect_predictions() %>% 
  roc_curve(fraud, .pred_Fraud, event_level = "second") %>% 
  mutate(model = "Random forest finale con ROSE")
```

```{r}
bind_rows(
  rf_auc,
  lr_auc,
  lr_auc_finale,
  lr_auc_finale_rose,
  rf_auc_finale,
  rf_auc_finale_rose
) %>% 
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  theme_bw() +
  # theme(legend.position = c(.8,.23)) +
  scale_color_viridis_d(option = "plasma", end = 1)
```