---
title: "Predictive Analytics Unit 4: Logistic Classifiers (revised)"
author: "Ken Arnold<br>Calvin University"
format: 
  revealjs:
    smaller: true
---


```{r setup, echo=FALSE}
include_graphics <- knitr::include_graphics
kable <- knitr::kable
glue <- glue::glue

library(conflicted)
conflict_prefer("filter", "dplyr", quiet = TRUE)
conflict_prefer("cols", "vroom", quiet = TRUE)
conflict_prefer("col_number", "vroom", quiet = TRUE)

knitr::opts_chunk$set(echo = TRUE, comment = "",
                      fig.width = 8,
                      fig.asp = 0.618,
                      out.width = "90%",
                      fig.align = "center",
                      dpi = 300,
                      message = FALSE)

ggplot2::theme_set(ggplot2::theme_bw(base_size = 16))

```

<style>
.col2 { column-count: 2; }
.tiny { font-size: 10px; }
</style>

```{r include=FALSE}
library(tidyverse)
library(tidymodels)
```

## Two Types of Supervised Learning

- **Regression**: fill in a *number*
  - How many customers?
  - What price?
- **Classification**: fill in a *choice*
  - How likely is this transaction to be fraudulent?
  - What is the *risk of default* for this credit application?


---

## Classification example: Blood test for autism?

::: {.col2}

We'll use an example from a [2017 PLOS Computational Biology paper](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005385)


```{r download-data, include=FALSE}
data_filename <- "data/autism.csv"
if (!file.exists(data_filename)) {
  dir.create("data")
  download.file("https://doi.org/10.1371/journal.pcbi.1005385.s001", data_filename)
}
```

```{r read-data, include=FALSE}
col_names <- names(read_csv(data_filename, n_max = 1, col_types = cols(.default = col_character())))
autism <- read_csv(data_filename, skip = 2, col_names = col_names, col_types = cols(
  .default = col_double(),
  Group = col_character()
))
```

```{r class.output="no-overflow-x"}
#| class-output: "no-overflow-x"
skimr::skim(autism)
```

:::

---

We have 3 kinds of data about `r nrow(autism)` children:

1. The outcome (`Group`): ASD (diagnosed with ASD), SIB (sibling not diagnosed with ASD), and NEU (age-matched neurotypical children, for control)

```{r group-counts}
autism %>% group_by(Group) %>% summarize(n = n()) %>% knitr::kable()
```

---

1. The outcome (`Group`): ASD, SIB, NEU
2. Concentrations of various metabolites in a blood sample:

```{r show-metabolites}
autism %>% select(-1, -last_col()) %>% names()
```

---

1. The outcome (`Group`): ASD, SIB, NEU
2. Concentrations of various metabolites in a blood sample
3. For the ASD children only, a measure of life skills ("Vineland ABC")

```{r show-behavior-score, fig.asp=.25}
autism %>% 
  ggplot(aes(x = `Vineland ABC`, y = Group)) + geom_boxplot()
```

---

## Exploratory Data Analysis (EDA)

What do these metabolites look like?

```{r bad-metabolite-plot, echo=FALSE, out.width="80%"}
library(ggridges)
autism %>%
  select(-`Vineland ABC`) %>% 
  pivot_longer(-Group, names_to = "Measure") %>% 
  ggplot(aes(x = value, y = Measure)) +
  geom_density_ridges() + 
  facet_wrap(vars(Group), scales = "free_x")
```

---

code for the previous plot:
```{r ref.label="bad-metabolite-plot", eval=FALSE}
```

---

## EDA

Better question: **Can these metabolites help us distinguish autism?**

---

```{r ridgeplot, fig.width=10, fig.asp=1, echo=FALSE, out.width="80%"}
autism %>%
  select(-`Vineland ABC`) %>% 
  pivot_longer(-Group, names_to = "Measure") %>% 
  ggplot(aes(x = value, y = Group)) +
  geom_density_ridges() +
  facet_wrap(vars(Measure), scales = "free_x") + 
  theme_ridges()
```

---

code for previous plot:

```{r ref.label="ridgeplot", eval=FALSE}
```

---

## Can we predict ASD vs non-ASD from metabolites?

* Let's start by (1) ignoring the behavior scores (that's an *outcome*) and comparing
just ASD and NEU.
* We need to drop SIB and encode `Group` as a factor. (Ensure that "ASD" is the *first* level.)

```{r}
data <- autism %>% 
  select(-`Vineland ABC`) %>% 
  filter(Group != "SIB") %>% 
  mutate(Group = as_factor(Group) %>% fct_relevel("ASD"))
```

```{r}
data$Group
```

---

## Linear Model for *Classification*

```{r fit-model}
model <- logistic_reg() %>% #<<
  fit(Group ~ `tGSH/GSSG` + `fGSH` + `Methion.`, data = data)
```


```{r show-coef}
model %>% tidy() %>% select(term, estimate) %>% kable()
```

```{r echo=FALSE, include=FALSE}
equatiomatic::extract_eq(extract_fit_engine(model), use_coefs = T)
```


$$
\log\left[ \frac { \widehat{P( \operatorname{..y} = \operatorname{NEU} )} }{ 1 - \widehat{P( \operatorname{..y} = \operatorname{NEU} )} } \right] = -15.34 + 0.14(\operatorname{tGSH/GSSG}) + 2.17(\operatorname{fGSH}) + 0.29(\operatorname{Methion.})
$$


---

### What do the *predictions* look like?

```{r predicted-probs-tree}
model %>% predict(data, type = "prob")
```

---

## Intuition: *Risk*

- First, linear model predicts a *risk score* (negative = low risk, positive = high risk).
  - Technically: *log-odds*: log(prob of *yes* / prob of *no*)
- Then, *logistic* function turns risk score into probability

Other examples:

- Risk of a customer churning
- Risk of a transaction being fraudulent
- Chance of one team winning: difference in skill scores (Elo)

---

### Were those predictions *good*?

::: {.small-code}

```{r plot-predicted-probs-tree, fig.asp=.4, out.width="90%"}
model %>%
  predict(data, type = "prob") %>%
  bind_cols(data) %>% 
  mutate(idx = row_number()) %>% 
  ggplot(aes(x = idx, y = .pred_ASD, color = Group, shape = Group)) +
  geom_hline(yintercept = .5) +
  geom_point() 
```

:::

---

### Quantifying Classification Quality

```{r}
decision_metrics <- yardstick::metric_set(accuracy, sensitivity,
                                 specificity, yardstick::f_meas)
model %>% 
  predict(data, type = "class") %>% 
  bind_cols(data) %>% 
  decision_metrics(truth = Group, estimate = .pred_class)
```

Note: these are all measures of the quality of the *decision*. Other metrics can evaluate the *scores*:

```{r}
model %>% 
  augment(data) %>% 
  roc_auc(truth = Group, .pred_ASD)
```


---

## Classification Metrics  {.smaller}

Note: this is *transposed* of the Wikipedia table <!-- TODO -->

|                      | Event happened              | No event happened           |
|----------------------|-------------------------------|-------------------------------|
| Event predicted    | True positive                 | False positive (Type 1 error) |
| No event predicted | False negative (Type 2 error) | True negative                 |

. . .

- **Accuracy** (% correct) = (TP + TN) / (# episodes)
- **False negative** ("miss") **rate** = FN / (# actual events)
- **False positive** ("false alarm") **rate** = FP / (# true non-events)

. . .

- **Sensitivity** ("true positive rate") = TP / (# actual events)
  - Sensitivity = 1 − False negative rate
- **Specificity** ("true negative rate") = TN / (# actual events)
  - Specificity = 1 − False positive rate
- [Wikipedia article](https://en.wikipedia.org/wiki/Sensitivity_and_specificity)

---

## Predicted probabilities tell confidence {.smaller}

- **Comparing models**: It's better to be *confident* when you're right, it's better to be *un-confident* when you're wrong.
- **Setting thresholds**
  - Model is usually trained as if false negative and false positive are equally bad.
  - But if false negative is worse than false positive, we can adjust the confidence threshold (default: 50%)

```{r}
metrics <- yardstick::metric_set(mn_log_loss, roc_auc)
model %>% predict(data, type = "prob") %>% 
  bind_cols(data) %>% 
  metrics(truth = Group, .pred_ASD)
```
