---
title: "🗓️ Week 4 - Classification, Metrics & Handling Imbalanced Data"
subtitle: "Theme: Binary Classification with Real-World Challenges"
author: "Dr. Ghita Berrada"
date: "24 October 2025"
format: 
  html:
    from: markdown+emoji
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

Welcome to Week 4! This week focuses on **binary classification** with a particular emphasis on **real-world challenges** like imbalanced datasets. We'll use a heart failure dataset to predict patient survival.

**:new: What's NEW this week:**

-   Deep dive into logistic regression theory (coefficients, perfect separation, regularization)
-   K-Nearest Neighbors (K-NN) classification algorithm
-   precision, recall, F1-score
-   ROC curves and AUC
-   **PR curves** - better than ROC for imbalanced data
-   **Calibration curves** - do predicted probabilities match reality?
-   Techniques for handling class imbalance (undersampling, SMOTE, class weights) - Feature selection methods (Information Gain, RFE, Lasso) - Cross-validation for hyperparameter tuning

**:books: Quick RECAP from last week (we'll build on these):**

-   Confusion matrices (TP, TN, FP, FN)
-   Basic metrics: accuracy\
-   Why accuracy can be misleading

**:date: Coming next week:**

-   Multi-class classification (3+ classes)
-   Advanced algorithms: SVM, Random Forest, XGBoost
-   One-vs-rest and one-vs-one strategies
-   Feature/model selection

------------------------------------------------------------------------

# :gear: Set up

**Loading R libraries**

```{r import, message=FALSE, warning=FALSE}
# Tidyverse packages
library(dplyr)
library(ggplot2)
library(readr)
library(tidyr)
library(purrr)
library(tibble)
library(ggridges)

# Tidymodels packages
library(parsnip)
library(yardstick)
library(recipes)
library(workflows)
library(rsample)
library(tune)
library(themis)  # For handling imbalanced data
library(kknn)    # For K-NN

# Additional libraries
library(FSelectorRcpp)  # For information gain
library(caret)          # For RFE
library(patchwork)      # For plot layouts
library(showtext)

# Enable showtext for better font rendering
showtext_auto()
```

**Loading the Heart Failure Clinical Records dataset**

This dataset was made publicly available on the [UCI machine learning repository](https://archive.ics.uci.edu/dataset/519/heart+failure+clinical+records). [^1]

[^1]: *"The UCI Machine Learning Repository is a collection of databases, domain theories, and data generators that are used by the machine learning community for the empirical analysis of machine learning algorithms"*

This dataset contains medical records of 299 patients who had heart failure. Each patient profile has 13 clinical features, and our target is whether they died during the follow-up period. This is a **binary classification** problem.

```{r load_data, message=FALSE}
# Load heart failure data from UCI repository
heart_data <- read_csv("https://archive.ics.uci.edu/ml/machine-learning-databases/00519/heart_failure_clinical_records_dataset.csv")

# Rename target variable for clarity and convert to factors
heart_data <- heart_data %>%
  rename(died = DEATH_EVENT) %>%
  mutate(
    died = factor(died, levels = c(0, 1), labels = c("Survived", "Died")),
    anaemia = factor(anaemia, levels = c(0, 1), labels = c("No", "Yes")),
    diabetes = factor(diabetes, levels = c(0, 1), labels = c("No", "Yes")),
    high_blood_pressure = factor(high_blood_pressure, levels = c(0, 1), 
                                 labels = c("No", "Yes")),
    sex = factor(sex, levels = c(0, 1), labels = c("Female", "Male")),
    smoking = factor(smoking, levels = c(0, 1), labels = c("No", "Yes"))
  )

# Quick look at the data
glimpse(heart_data)
```

**Data dictionary**:

-   `age`: Age of patient (years)
-   `anaemia`: Decrease of red blood cells or hemoglobin
-   `creatinine_phosphokinase`: Level of CPK enzyme (mcg/L)
-   `diabetes`: Whether patient has diabetes
-   `ejection_fraction`: Percentage of blood leaving heart at each contraction (%)
-   `high_blood_pressure`: Whether patient has hypertension
-   `platelets`: Platelets in blood (kiloplatelets/mL)
-   `serum_creatinine`: Level of serum creatinine (mg/dL)
-   `serum_sodium`: Level of serum sodium (mEq/L)
-   `sex`: Male or female
-   `smoking`: Whether patient smokes
-   `time`: Follow-up period (days)
-   **`died`**: TARGET VARIABLE - Whether patient died during follow-up

------------------------------------------------------------------------

## Understanding Class Imbalance

```{r explore_imbalance}
# Check class distribution - THIS IS CRITICAL!
class_counts <- table(heart_data$died)
class_props <- prop.table(class_counts)

cat("Class Distribution:\n")
print(class_counts)
cat("\nProportions:\n")
print(round(class_props, 3))

# Visualize
ggplot(heart_data, aes(x = died, fill = died)) +
  geom_bar() +
  geom_text(stat = 'count', aes(label = after_stat(count)), vjust = -0.5, size = 5) +
  geom_text(stat = 'count', 
            aes(label = paste0("(", round(after_stat(count)/sum(after_stat(count))*100, 1), "%)")),
            vjust = 1.5, size = 4) +
  labs(title = "Class Distribution: Heart Failure Dataset",
       subtitle = "Notice the imbalance: ~68% survived, ~32% died",
       y = "Count",
       x = "Outcome") +
  scale_fill_manual(values = c("Survived" = "steelblue", "Died" = "coral")) +
  theme_minimal() +  # Use a general font
  theme(legend.position = "none", text = element_text(size = 12))
        
```

**:warning: CRITICAL OBSERVATION**: This dataset is **imbalanced**!

-   **Majority class (Survived)**: \~68% of patients
-   **Minority class (Died)**: \~32% of patients

This imbalance ratio (\~2:1) mirrors real-world healthcare data where fortunately, survival is more common than death. However, this creates challenges:

1.  Models tend to predict the majority class (easy accuracy boost)
2.  Minority class (death) is often MORE important to detect
3.  Standard metrics like accuracy become misleading
4.  Decision boundaries get biased toward majority class

------------------------------------------------------------------------

# Part 1: Logistic Regression - Theory Deep Dive

## 1.1 Why Not Linear Regression for Classification? :no_entry:️

Last week, we touched on this briefly. Let's understand it deeply now.

```{r why_not_linear, fig.width=10, fig.height=5}
# -------------------------------------------------------------------
# 1. Prepare data
# -------------------------------------------------------------------
simple_data <- heart_data %>%
  select(age, died) %>%
  mutate(died = factor(died, levels = c("Survived", "Died")))  # ensure correct order

age_grid <- tibble(age = seq(20, 100, 1))

# -------------------------------------------------------------------
# 2. WRONG MODEL: Linear regression for classification
# -------------------------------------------------------------------
# For linear regression, we need a numeric outcome
simple_data_lin <- simple_data %>%
  mutate(died_numeric = as.numeric(died) - 1)

linear_recipe <- recipe(died_numeric ~ age, data = simple_data_lin)

linear_spec <- linear_reg() %>%
  set_engine("lm")

linear_wf <- workflow() %>%
  add_model(linear_spec) %>%
  add_recipe(linear_recipe)

linear_fit <- fit(linear_wf, data = simple_data_lin)

pred_linear <- predict(linear_fit, new_data = age_grid) %>%
  bind_cols(age_grid)

# -------------------------------------------------------------------
# 3. CORRECT MODEL: Logistic regression
# -------------------------------------------------------------------
logistic_recipe <- recipe(died ~ age, data = simple_data)

logistic_spec <- logistic_reg() %>%
  set_engine("glm")

logistic_wf <- workflow() %>%
  add_model(logistic_spec) %>%
  add_recipe(logistic_recipe)

logistic_fit <- fit(logistic_wf, data = simple_data)

pred_logistic <- predict(logistic_fit, new_data = age_grid, type = "prob") %>%
  bind_cols(age_grid)

# -------------------------------------------------------------------
# 4. Plot comparisons
# -------------------------------------------------------------------
p1 <- ggplot(simple_data_lin, aes(x = age, y = died_numeric)) +
  geom_point(alpha = 0.4, position = position_jitter(height = 0.02), size = 2) +
  geom_line(data = pred_linear, aes(x = age, y = .pred),
            color = "red", linewidth = 1.2) +
  geom_hline(yintercept = c(0, 1), linetype = "dashed", alpha = 0.5) +
  annotate("rect", xmin = 20, xmax = 50, ymin = -0.2, ymax = 0,
           fill = "red", alpha = 0.2) +
  annotate("text", x = 35, y = -0.1, label = "Impossible!\nProbability < 0",
           color = "red", size = 3) +
  labs(title = "Linear Regression for Classification: THE PROBLEM",
       y = "Death (0 = Survived, 1 = Died)",
       x = "Age") +
  theme_minimal()

p2 <- ggplot(simple_data_lin, aes(x = age, y = died_numeric)) +
  geom_point(alpha = 0.4, position = position_jitter(height = 0.02), size = 2) +
  geom_line(data = pred_logistic, aes(x = age, y = .pred_Died),
            color = "steelblue", linewidth = 1.2) +
  geom_hline(yintercept = c(0, 1), linetype = "dashed", alpha = 0.5) +
  annotate("text", x = 35, y = 0.9,
           label = "Logistic regression\nalways stays in [0,1]",
           color = "steelblue", size = 3) +
  labs(title = "Logistic Regression: THE SOLUTION",
       y = "Predicted Probability of Death",
       x = "Age") +
  theme_minimal()

p1 + p2
```

**The fundamental problem with linear regression**:

-   Predictions can be \< 0 or \> 1
-   These aren't valid probabilities!
-   Makes no sense to say "probability of death = -0.15" or "1.3"

**The solution: Logistic regression**

-   Uses the **sigmoid function** to "squash" predictions into \$\$0,1\]
-   Always outputs valid probabilities
-   S-shaped curve matches reality better

------------------------------------------------------------------------

## 1.2 The Sigmoid Function - Mathematical Foundation

The **logistic (sigmoid) function** is the heart of logistic regression:

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

Where $z = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p$ is the linear combination of features.

For binary classification: $$
P(Y = 1 | X) = \sigma(\beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p)
$$

**Key properties**:

1.  **Output always in (0, 1)** - perfect for probabilities
2.  **Smooth and differentiable** - good for optimization
3.  **S-shaped curve** - models real-world binary outcomes well
4.  **Symmetric around 0.5** - when $z=0$, $P=0.5$

```{r sigmoid_viz, fig.width=10, fig.height=4}
# Create sigmoid visualization
sigmoid <- function(z) 1 / (1 + exp(-z))

p1 <- tibble(z = seq(-6, 6, 0.1)) %>%
  mutate(p = sigmoid(z)) %>%
  ggplot(aes(z, p)) +
  geom_line(linewidth = 1.5, color = "steelblue") +
  geom_hline(yintercept = c(0, 0.5, 1), linetype = "dashed", alpha = 0.3) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.3) +
  geom_point(aes(x = 0, y = 0.5), size = 4, color = "red") +
  annotate("text", x = 0, y = 0.55, label = "p = 0.5 when z = 0", 
           hjust = 0, vjust = 0, color = "red") +
  labs(title = "The Sigmoid Function σ(z)",
       subtitle = "Transforms any real number to probability in (0,1)",
       x = expression(z == beta[0] + beta[1]*X[1] + " ... + " + beta[p]*X[p]),
       y = "σ(z) = Predicted Probability") +
  theme_minimal() +
  theme(text = element_text(size = 11))

# Show derivative (for those interested)
sigmoid_derivative <- function(z) {
  s <- sigmoid(z)
  s * (1 - s)
}

p2 <- tibble(z = seq(-6, 6, 0.1)) %>%
  mutate(derivative = sigmoid_derivative(z)) %>%
  ggplot(aes(z, derivative)) +
  geom_line(linewidth = 1.5, color = "coral") +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.3) +
  annotate("text", x = 0, y = 0.27, 
           label = "Maximum change\noccurs at z = 0", 
           hjust = 0.5, vjust = 1, color = "coral", size = 3) +
  labs(title = "Sigmoid Derivative σ'(z)",
       subtitle = "Shows where probability changes most rapidly",
       x = "z",
       y = "Rate of Change") +
  theme_minimal() +
  theme(text = element_text(size = 11))

p1 + p2
```

::: callout-note
## The Sigmoid Function and Its Sensitivity

The **sigmoid function** (left plot) is what logistic regression uses to convert the model’s linear output ($z = \beta\_0 + \beta\_1 X_1 + \dots + \beta\_p X_p$) into a probability between 0 and 1:

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

This S-shaped curve is **steepest around (**$z = 0$), where the predicted probability is about 0.5.\
In this region, small changes in the predictors lead to **large changes in predicted probability** — the model is most sensitive here.

The **derivative** of the sigmoid function shows how quickly the probability changes with (z):

$$\frac{d\sigma}{dz} = \sigma(z)\big(1 - \sigma(z)\big)$$

This derivative reaches its **maximum value when (**$z = 0$) (i.e., ($P = 0.5$ where $P$ is the probability of death)), as shown in the right plot.\
That means the model’s output changes most rapidly when it is uncertain between the two classes.\
As (z) becomes very positive or very negative, the sigmoid flattens out, the derivative approaches zero, and the model’s predictions become **less sensitive** — it is already confident about the outcome.
:::

------------------------------------------------------------------------

## 1.3 Understanding Coefficients: From Linear to Log-Odds to Odds

This is where logistic regression gets interesting! Let's unpack the layers:

### Layer 1: The Linear Model (Log-Odds)

Logistic regression actually models the **log-odds**:

$$
\log\left(\frac{P(Y=1)}{1-P(Y=1)}\right) = \beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p
$$

Where:

-   $P(Y=1)$ = probability of death
-   $\frac{P(Y=1)}{1-P(Y=1)}$ = **odds** of death
-   $\log(\text{odds})$ = **log-odds** (or **logit**)

**Why log-odds?** Because the log-odds can range from $-\infty$ to $+\infty$ (just like linear regression), while probabilities are constrained to \$\$0,1\].

### Layer 2: Understanding Odds

**Odds** = ratio of probability of event happening to probability of it not happening

$$
\text{Odds} = \frac{P}{1-P}
$$

**Examples**:

-   $P = 0.5$ → Odds = 0.5/0.5 = 1 (even odds, "50-50")
-   $P = 0.75$ → Odds = 0.75/0.25 = 3 (3:1, "3 times more likely to happen")
-   $P = 0.25$ → Odds = 0.25/0.75 = 0.33 (1:3, "3 times more likely NOT to happen")

### Layer 3: Interpreting Coefficients as Odds Ratios

**Key insight**: Exponentiating a coefficient gives the **odds ratio**:

$$
\text{Odds Ratio} = e^{\beta_j}
$$

**Interpretation**:

-   OR = 1: No effect (odds don't change)
-   OR \> 1: Predictor **increases** odds of outcome (e.g., OR = 2 means odds double)
-   OR \< 1: Predictor **decreases** odds of outcome (e.g., OR = 0.5 means odds halve)

Let's see this in practice!

------------------------------------------------------------------------

## 1.4 Fitting and Interpreting a Logistic Regression Model

First, split the data:

```{r split_data}
set.seed(123)

# Stratified split to maintain class proportions
heart_split <- initial_split(heart_data, prop = 0.75, strata = died)
heart_train <- training(heart_split)
heart_test <- testing(heart_split)

# Verify proportions are maintained
cat("Training set class distribution:\n")
print(prop.table(table(heart_train$died)))

cat("\nTest set class distribution:\n")
print(prop.table(table(heart_test$died)))
```

**Good practice**: Use `strata = died` to ensure both training and test sets have similar class distributions!

Now fit a simple model with just 3 key predictors:

```{r simple_logistic}
# Create recipe with preprocessing
simple_recipe <- recipe(died ~ ejection_fraction + serum_creatinine + age, 
                        data = heart_train) %>%
  step_normalize(all_numeric_predictors())  # Standardize for easier interpretation

# Specify logistic regression model
logistic_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

# Create workflow
simple_workflow <- workflow() %>%
  add_recipe(simple_recipe) %>%
  add_model(logistic_spec)

# Fit model
simple_fit <- fit(simple_workflow, data = heart_train)

# Extract and interpret coefficients
coef_table <- simple_fit %>%
  extract_fit_engine() %>%
  broom::tidy() %>%
  mutate(
    odds_ratio = exp(estimate),
    interpretation = case_when(
      term == "(Intercept)" ~ "Baseline log-odds",
      odds_ratio > 1.5 ~ "Strong positive effect on death",
      odds_ratio > 1.1 ~ "Moderate positive effect on death",
      odds_ratio < 0.67 ~ "Strong negative effect on death (protective)",
      odds_ratio < 0.91 ~ "Moderate negative effect on death (protective)",
      TRUE ~ "Minimal effect"
    )
  ) %>%
  select(term, estimate, odds_ratio, p.value, interpretation)

print(coef_table, n = Inf)
```

::: callout-note
**Why this model?**

Before interpreting the model output, it’s important to understand **why we’re fitting this particular logistic regression**.

Our goal here is to estimate how a few key clinical variables — **ejection fraction**, **serum creatinine**, and **age** — relate to the likelihood of **death** among patients with heart failure.

Because the outcome (`died`) is binary (yes/no), we use **logistic regression**, which models the *log-odds* of the outcome as a linear combination of predictors:

$$
\text{logit}(P(\text{death})) = \beta_0 + \beta_1(\text{ejection fraction}) + \beta_2(\text{serum creatinine}) + \beta_3(\text{age})
$$

Each coefficient ($\beta_i$) represents how much the log-odds of death change for a one-unit increase in the predictor, holding the others constant.

We normalize the predictors first so that the coefficients are on a **comparable scale**, making it easier to judge which variables have stronger or weaker effects.

The model’s output will show us:

-   how strongly each variable is associated with the probability of death,\
-   whether the effect increases or decreases that probability, and\
-   how confident we are in those estimates (via p-values).

This conceptual grounding will help us interpret the coefficient table.
:::

::: callout-tip
## 🕵 How to read this coefficient table

### Understanding probabilities, odds, and odds ratios

**Starting with probabilities:**

-   **Probability** tells us how likely an event is: P(death) = 0.2 means 20% chance of death
-   Ranges from 0 to 1 (or 0% to 100%)

**Moving to odds:**

-   **Odds** compare the probability of an event *happening* vs *not happening*: $$ \text{odds} = \frac{P(\text{death})}{1 - P(\text{death})}$$
-   Example: If P(death) = 0.2, then odds = 0.2/0.8 = 0.25 (or "1 to 4")
-   Translation: "For every 1 death, there are 4 survivors" - Odds range from 0 to ∞ (unlike probabilities!)

**Understanding odds ratios:**

-   An **odds ratio (OR)** compares odds between two groups or conditions
-   It tells us how odds *change* when a predictor increases $$ \text{OR} = \frac{\text{new odds}}{\text{old odds}} \quad \text{or} \quad \text{new odds} = \text{old odds} \times \text{OR}$$

**Quick reference table:**

| Probability | Odds | Interpretation                 |
|-------------|------|--------------------------------|
| 0.10 (10%)  | 0.11 | 1 death per 9 survivors        |
| 0.20 (20%)  | 0.25 | 1 death per 4 survivors        |
| 0.50 (50%)  | 1.00 | 1 death per 1 survivor (equal) |
| 0.75 (75%)  | 3.00 | 3 deaths per 1 survivor        |
| 0.90 (90%)  | 9.00 | 9 deaths per 1 survivor        |

------------------------------------------------------------------------

### What each column means

-   **Estimate (β):** Change in the *log-odds* of death for a 1 SD increase in the predictor

    -   Log-odds are harder to interpret directly, so we exponentiate to get odds ratios!

-   **Odds Ratio (OR = e\^β):**\

    How the *odds* of death change multiplicatively for a 1 SD increase:

    -   **OR = 1:** No effect (odds unchanged)
    -   **OR \> 1:** Odds increase (risk factor)
    -   **OR \< 1:** Odds decrease (protective factor)

-   **p-value:** Statistical evidence for an association (\< 0.05 typically significant)

------------------------------------------------------------------------

### Worked examples from this model

#### Example 1: Ejection fraction (OR = 0.51) — Protective factor

**Starting scenario (baseline from the model's intercept):**

The model's **intercept** represents the log-odds of death when all standardized predictors equal zero - that is, when all predictors are at their **mean values** in the training data. This gives us our baseline:

-   Model intercept (β₀): -0.75 (example value)
-   Baseline log-odds: -0.75
-   Baseline odds: e\^(-0.75) = 0.47
-   **Baseline probability: 0.47/(1 + 0.47) = 0.32 (or 32%)**

**Interpretation of baseline:** An "average" patient

-   someone with average ejection fraction, average creatinine, average age, etc.
-   has a **32% probability of death** according to our model.

------------------------------------------------------------------------

**After 1 SD increase in ejection fraction** (all other predictors stay at their means):

-   New odds = 0.47 × 0.51 = 0.24
-   New probability = 0.24/(1 + 0.24) = 0.19 (or 19%)

**Interpretation:**

-   The odds of death are **multiplied by 0.51** (or reduced by 49%)
-   Probability dropped from 32% to 19% — a **13 percentage point** decrease
-   **Clinical meaning:** Improving ejection fraction by 1 SD (from average to above-average) is protective

------------------------------------------------------------------------

**How does baseline risk affect the interpretation?**

The odds ratio (0.51) tells us the **relative effect**, but the **absolute benefit** depends on where you start:

-   **Starting at 32% (average risk):** 32% → 19% (13 point decrease)
-   **Starting at 60% (high risk):** 60% → 43% (17 point decrease)
-   **Starting at 10% (low risk):** 10% → 5% (5 point decrease)

**Key insight:** The same intervention (improving ejection fraction by 1 SD) has the **same relative effect (OR = 0.51)** but larger **absolute benefits** for higher-risk patients.

------------------------------------------------------------------------

#### Example 2: Serum creatinine (OR = 1.61) — Risk factor

**Starting from the same baseline** (average patient with 32% risk):

-   Baseline odds: 0.47

**After 1 SD increase in serum creatinine** (all other predictors at their means):

-   New odds = 0.47 × 1.61 = 0.76
-   New probability = 0.76/(1 + 0.76) = 0.43 (or 43%)

**Interpretation:**

-   Odds multiply by 1.61 (or increase by 61%)
-   Probability: 32% → 43% (11 percentage point increase)
-   **Clinical meaning:** Kidney dysfunction (1 SD above average creatinine) substantially increases death risk

------------------------------------------------------------------------

#### Example 3: Age (OR = 1.89) — Strong risk factor

**Starting from baseline** (32% risk for average patient): - Baseline odds: 0.47

**After 1 SD increase in age:**

-   New odds = 0.47 × 1.89 = 0.89
-   New probability = 0.89/(1 + 0.89) = 0.47 (or 47%)

**Interpretation:**

-   Odds multiply by 1.89 (or increase by 89%)
-   Probability: 32% → 47% (15 percentage point increase)
-   **Clinical meaning:** Older age (1 SD above average) is a strong risk factor

------------------------------------------------------------------------

### :dart: Key takeaways

1.  **Odds ratios are multiplicative:** OR = 2 means odds *double*, OR = 0.5 means odds *halve*

2.  **Odds ≠ Probability:**

    -   Doubling the odds does NOT double the probability
    -   Effect on probability depends on baseline risk
    -   Same OR produces different absolute probability changes at different baseline risks

3.  **Converting back to probability:** $$ P = \frac{\text{odds}}{1 + \text{odds}}$$

4.  **Standardized predictors matter:**

    -   All interpretations are per **1 SD increase**, not per 1 unit
    -   Makes effects comparable across different scales

5.  **Clinical context is crucial:**

    -   OR tells us about association strength (relative effect)
    -   Clinical significance depends on baseline risk (absolute effect) and practical impact
:::

If you want to review/understand the concepts linked to logistic regression, you could have a look at [this Youtube playlist](https://www.youtube.com/watch?v=8nm0G-1uJzA&list=PLblh5JKOoLUKxzEP5HA2d-Li7IJkHfXSe)

------------------------------------------------------------------------

## :small_blue_diamond: Section 1.5 – The Problem: Perfect Separation

In logistic regression, we estimate probabilities using a smooth *S-shaped* curve that separates two outcome classes (e.g., “Yes” vs “No”).

However, sometimes a predictor (or a combination of predictors) can **perfectly separate** the classes, meaning there exists a cutoff value that divides all 0s from all 1s with no overlap.

This situation is called **Perfect Separation** (or *quasi-perfect separation* if the separation is almost perfect).

### :wrench: Example: Constructing a perfectly separated dataset

```{r perfect_sep_demo, fig.width=10, fig.height=4, message=FALSE, warning=FALSE}

# Create a toy dataset with perfect separation
perfect_sep_example <- tibble(
  x = c(1, 2, 3, 4, 5, 5.5, 6, 7, 8, 9, 10),
  y = factor(c(0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1),
             levels = c(0, 1),
             labels = c("No", "Yes"))
)

# Plot the data and the perfect separating line
p1 <- ggplot(perfect_sep_example, aes(x = x, y = as.numeric(y) - 1, color = y)) +
  geom_point(size = 4) +
  geom_vline(xintercept = 5.5, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = 5.5, y = 0.5,
           label = "Perfect separator:\nx < 5.5 → Class = No\nx ≥ 5.5 → Class = Yes",
           hjust = -0.1, color = "red", size = 3.5) +
  labs(title = "Perfect Separation Example",
       subtitle = "Predictor X perfectly divides classes",
       x = "Predictor X", y = "Observed Class (0 or 1)") +
  scale_color_manual(values = c("No" = "steelblue", "Yes" = "coral")) +
  theme_minimal() +
  theme(legend.position = "none")

p1
```

:gear: What happens when we fit logistic regression?

If we fit a standard logistic regression to this data:

``` r
# Fit logistic regression model
perf_sep_model <- glm(y ~ x, data = perfect_sep_example, family = binomial)

# View estimated coefficients
coef(perf_sep_model)
```

You’ll often see a warning such as:

> “fitted probabilities numerically 0 or 1 occurred” “algorithm did not converge”

This happens because the model can fit *perfectly* by making the slope coefficient (β₁) extremely large — approaching infinity.

### :chart_with_upwards_trend: Visualizing the unstable fitted curve

```{r perfect_sep_curve, message=FALSE, warning=FALSE}
# Generate predicted probabilities over a grid
pred_grid <- tibble(x = seq(0, 11, 0.1)) %>%
  mutate(prob = predict(perf_sep_model, newdata = ., type = "response"))
beta1_val <- coef(perf_sep_model)[2]
label_text <- sprintf("beta[1] == %.1f * ' (extremely~large!)'", beta1_val)

# Plot the logistic curve vs data
p2 <- ggplot(perfect_sep_example, aes(x = x, y = as.numeric(y) - 1, color = y)) +
  geom_point(size = 4) +
  geom_line(data = pred_grid, aes(x = x, y = prob), color = "black", linewidth = 1) +
  annotate(
  "text",
  x = 6.5, y = 0.7,
  label = label_text,
  color = "red", size = 3.5,
  parse = TRUE
) +
  labs(title = "Fitted Logistic Curve under Perfect Separation",
       subtitle = "Coefficients blow up, probabilities approach 0 or 1",
       x = "Predictor X", y = "Predicted Probability") +
  scale_color_manual(values = c("No" = "steelblue", "Yes" = "coral")) +
  theme_minimal() +
  theme(legend.position = "none")

p1 + p2
```

### :warning: Why this is problematic

| Problem | Consequence |
|-------------------------------|-----------------------------------------|
| Coefficients go to ±∞ | Model unstable, doesn’t converge |
| Standard errors inflate | Large uncertainty; unreliable inference |
| Predicted probabilities = 0 or 1 | No nuance; overconfident predictions |
| Overfitting | Great on training data, terrible on new data |
| Small data changes | Cause massive swings in coefficients |

In short, *perfect separation breaks the math behind logistic regression*.

------------------------------------------------------------------------

## :bulb: Key Takeaway

> When predictors perfectly (or almost perfectly) separate outcomes, logistic regression can’t find finite, stable parameter estimates.
>
> **Solution:** Use *regularization* to control coefficient growth.

## :small_blue_diamond: Section 1.6 – The Solution: Regularization

Regularization stabilizes logistic regression by **penalizing large coefficients**, preventing them from shooting off to infinity (as in perfect separation).

### :abacus: Basic idea

We modify the objective function:

$$ \text{Loss} = \text{Logistic Loss} + \lambda \times \text{Penalty}(\boldsymbol{\beta})$$

-   **λ** controls the strength of the penalty (larger → more shrinkage)\
-   The **Penalty** depends on the chosen regularization type

### :1234: Common types of regularization

| Type | Penalty | Effect on Coefficients | Best Used When… |
|-----------------|-----------------------|-----------------|-----------------|
| **Ridge (L2)** | ($\lambda \sum \beta_j^2$) | Shrinks all toward zero (none exactly zero) | All predictors might be relevant |
| **Lasso (L1)** | ($\lambda \sum | \beta\_j |$) |  |  |
| **Elastic Net** | ($\lambda (\alpha \sum \|\beta_j\| + (1-\alpha) \sum \beta_j^2)$) | Blend of L1 and L2 | You want balance between sparsity and stability |

### :toolbox: Demonstration: Comparing models with and without regularization

```{r regularization_demo, fig.width=10, fig.height=5, message=FALSE, warning=FALSE}

# Preprocessing: standardize numeric predictors
full_recipe <- recipe(died ~ ., data = heart_train) %>%
  step_dummy(all_nominal_predictors()) %>% #the glmnet engine expects numeric values!
  step_normalize(all_numeric_predictors())

# 1. Unregularized Logistic Regression
unreg_spec <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification")

unreg_fit <- workflow() %>%
  add_recipe(full_recipe) %>%
  add_model(unreg_spec) %>%
  fit(data = heart_train)

# 2. Ridge (L2)
ridge_spec <- logistic_reg(penalty = 0.1, mixture = 0) %>%
  set_engine("glmnet") %>%
  set_mode("classification")

ridge_fit <- workflow() %>%
  add_recipe(full_recipe) %>%
  add_model(ridge_spec) %>%
  fit(data = heart_train)

# 3. Lasso (L1)
lasso_spec <- logistic_reg(penalty = 0.1, mixture = 1) %>%
  set_engine("glmnet") %>%
  set_mode("classification")

lasso_fit <- workflow() %>%
  add_recipe(full_recipe) %>%
  add_model(lasso_spec) %>%
  fit(data = heart_train)

# 4. Elastic Net (α = 0.5)
elastic_spec <- logistic_reg(penalty = 0.1, mixture = 0.5) %>%
  set_engine("glmnet") %>%
  set_mode("classification")

elastic_fit <- workflow() %>%
  add_recipe(full_recipe) %>%
  add_model(elastic_spec) %>%
  fit(data = heart_train)

# Helper function to extract coefficients
extract_coefs <- function(fit, model_name) {
  if (model_name == "Unregularized") {
    tidy(extract_fit_engine(fit)) %>%
      mutate(model = model_name) %>%
      select(model, term, estimate)
  } else {
    coef_mat <- extract_fit_engine(fit) %>% coef(s = 0.1)
    tibble(term = rownames(as.matrix(coef_mat)),
           estimate = as.numeric(coef_mat)) %>%
      mutate(model = model_name)
  }
}

# Combine all coefficients
all_coefs <- bind_rows(
  extract_coefs(unreg_fit, "Unregularized"),
  extract_coefs(ridge_fit, "Ridge (L2)"),
  extract_coefs(lasso_fit, "Lasso (L1)"),
  extract_coefs(elastic_fit, "Elastic Net")
) %>%
  filter(term != "(Intercept)")

# Visualize coefficient shrinkage
ggplot(all_coefs, aes(x = term, y = estimate, fill = model)) +
  geom_col(position = "dodge") +
  coord_flip() +
  labs(title = "Effect of Regularization on Coefficient Estimates",
       subtitle = "Regularization shrinks coefficients; Lasso can zero out some predictors",
       x = "Predictor Variable", y = "Coefficient Value", fill = "Model Type") +
  theme_minimal() +
  theme(text = element_text(size = 10))
```

### :bar_chart: Interpretation of results

-   **Unregularized (GLM):**\
    Coefficients can take **large or unstable values**, especially when predictors are correlated or the sample size is limited.\
    This model fits the training data most closely but risks **overfitting**.

-   **Ridge (L2 penalty):**\
    Shrinks all coefficients smoothly toward zero but never fully to zero.\
    This reduces variance and improves model **stability**, while keeping all predictors in the model.

-   **Lasso (L1 penalty):**\
    Shrinks some coefficients exactly to zero, effectively **removing less important variables**.\
    This yields a simpler, more **interpretable** model by performing **automatic feature selection**.

-   **Elastic Net (α = 0.5):**\
    Combines Ridge’s stability and Lasso’s sparsity.\
    Useful when predictors are correlated — it can select groups of related variables rather than just one.

**In the plot above:**

Regularization gradually pulls coefficient magnitudes toward zero.\
Notice how **Lasso** has dropped some variables entirely, while **Ridge** only compresses their values.\
**Elastic Net** lies in between, showing both shrinkage and sparsity.

### :compass: Practical guidelines

| Situation                             | Recommended Regularization     |
|---------------------------------------|--------------------------------|
| Many correlated predictors            | **Ridge** or **Elastic Net**   |
| Only some predictors likely important | **Lasso**                      |
| Unsure which scenario applies         | **Elastic Net** (often safest) |

### :white_check_mark: Key Takeaway

> Regularization stabilizes logistic regression under perfect or near-perfect separation by preventing coefficient explosion. It improves generalization, reduces variance, and can even perform automatic feature selection.

### How does our dataset do in terms of separation?

```{r visualize_imbalance}
# Visualize class imbalance across key features
heart_train %>%
  select(died, age, ejection_fraction, serum_creatinine) %>%
  pivot_longer(-died, names_to = "feature", values_to = "value") %>%
  ggplot(aes(x = value, y = died, fill = died)) +
  geom_density_ridges(alpha = 0.7) +
  facet_wrap(~feature, scales = "free_x") +
  labs(title = "Feature Distributions by Outcome",
       subtitle = "Notice the overlap - no perfect separation") +
  theme_minimal() +
  theme(legend.position = "top")
```

------------------------------------------------------------------------

# Part 2: :books: Recap: Confusion Matrix & Accuracy

Last week, we introduced the **confusion matrix** and **accuracy** as our first evaluation metrics. We also discussed why accuracy can be misleading for imbalanced data. Let's quickly recap before expanding our toolkit.

## 2.1 Confusion Matrix (recap)

```{r confusion_matrix_recap}
# Get predictions from simple model
simple_pred <- predict(simple_fit, heart_test, type = "class") %>%
  bind_cols(heart_test %>% select(died))

# Create confusion matrix
cm <- conf_mat(simple_pred, truth = died, estimate = .pred_class)
print(cm)

# Visualize
cm %>%
  autoplot(type = "heatmap") +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(
    title = "Confusion Matrix",
    subtitle = "Rows = Actual (Truth), Columns = Predicted"
  ) +
  theme_minimal() %>%
  print() #In a Quarto document, however, only objects explicitly printed are shown. Wrapping the final plot in `print()` ensures it always appears in the rendered output.
```

**Anatomy** (as a reminder):

|   | **Predicted: Survived** | **Predicted: Died** |
|------------------------|--------------------------|-----------------------|
| **Actual: Survived** | True Negative (TN) :white_check_mark: | False Positive (FP) :x: |
| **Actual: Died** | False Negative (FN) :x: | True Positive (TP) :white_check_mark: |

-   **TN**: Correctly predicted survived :white_check_mark:
-   **TP**: Correctly predicted died :white_check_mark:
-   **FP**: Predicted died but actually survived :x: (Type I error, "false alarm")
-   **FN**: Predicted survived but actually died :x: (Type II error, "missed detection")

**Which errors matter more?**

-   Depends on the application!
-   In healthcare: Missing a death (FN) is usually worse than a false alarm (FP)
-   In spam detection: False alarms (FP) might be worse than missing spam (FN)

::: callout-note
## How to read this confusion matrix

-   **Diagonal cells** (48 and 12) show correct predictions.\
-   **Off-diagonal cells** (12 and 3) show errors.\
-   The model correctly identifies most survivors,\
    but **some patients who actually died were predicted to survive** —\
    suggesting a tendency to err on the side of predicting survival.
:::

**Interpreting the cells**

| Cell | Meaning | What it implies |
|-------------------|---------------------------|--------------------------|
| **48 (top-left)** | People who *survived* and were correctly predicted to survive. | The model is good at identifying survivors. |
| **12 (top-right)** | People who *died* but were predicted to survive. | These are **missed deaths** — the model failed to flag them. |
| **3 (bottom-left)** | People who *survived* but were predicted to die. | A few **false alarms**, where the model was overly cautious. |
| **12 (bottom-right)** | People who *died* and were correctly predicted to die. | The model did catch some true deaths. |

**:speech_balloon: Summary interpretation**

> The model predicts survival correctly most of the time, but it **misses several patients who actually died**. It’s more confident about recognizing those who survive than those who die. This pattern is common when the dataset is **imbalanced**, with more survivors than deaths.

------------------------------------------------------------------------

## 2.2 Accuracy and why it fails (recap)

Last week, we learned about **accuracy**:

$$
\text{Accuracy} = \frac{TP + TN}{TP + TN + FP + FN} = \frac{\text{Correct predictions}}{\text{Total predictions}}
$$

Accuracy seems intuitive, just the proportion of correct predictions. But let's see the **critical problem** with imbalanced data:

```{r accuracy_problem}
# Calculate our model's accuracy
model_acc <- accuracy(simple_pred, truth = died, estimate = .pred_class)$.estimate

# Create a "stupid" baseline that ALWAYS predicts "Survived"
stupid_baseline <- tibble(
  died = heart_test$died,
  .pred_class = factor(rep("Survived", nrow(heart_test)), 
                      levels = c("Survived", "Died"))
)

stupid_acc <- accuracy(stupid_baseline, truth = died, estimate = .pred_class)$.estimate

# Compare
cat("Our logistic regression accuracy:", sprintf("%.1f%%\n", model_acc * 100))
cat("'Always predict Survived' baseline:", sprintf("%.1f%%\n\n", stupid_acc * 100))

cat("⚠️ CRITICAL PROBLEM:\n")
cat("A model that NEVER detects death still achieves ~68% accuracy!\n")
cat("This is because 68% of patients survived in our dataset.\n\n")

cat("Why this is dangerous:\n")
cat("- High accuracy doesn't mean the model is useful\n")
cat("- It could be ignoring the minority class (deaths) completely\n")
cat("- We need better metrics that don't reward lazy predictions!\n")
```

::: callout-warning
**:warning: Critical takeaway: Accuracy ≠ Usefulnes**

-   The model achieves **80% accuracy**, but most of that comes from correctly predicting the majority class (*Survived*).\
-   A trivial model that always predicts *Survived* would already reach **68% accuracy**.\
-   The confusion matrix reveals the problem: **many deaths are missed** (half of them!).

Accuracy alone can hide serious weaknesses — especially when one class (e.g., survival) dominates the data.
:::

**The lesson**: With imbalanced data, a model can achieve high accuracy by simply always predicting the majority class. This makes accuracy **dangerously misleading** when classes are imbalanced.

**We need better metrics!** This is why today we'll learn about precision, recall, F1-score, and metrics specifically designed for imbalanced data.

------------------------------------------------------------------------

# Part 3: :new: Beyond Accuracy - Better Metrics for Classification

## 3.1 :new: Precision, Recall, and F1-Score

Since accuracy fails us with imbalanced data, we need metrics that better capture performance on both classes, especially the minority class.

```{r precision_recall_f1}
# Calculate comprehensive metrics
comprehensive_metrics <- yardstick::metric_set(yardstick::precision, yardstick::recall, f_meas, bal_accuracy)
results <- simple_pred %>%
  comprehensive_metrics(truth = died, estimate = .pred_class, event_level = "second")

print(results)
```

### Precision - "When we predict positive, are we usually right?"

$$
\text{Precision} = \frac{TP}{TP + FP} = \frac{TP}{\text{all predicted positives}}
$$

**What it tells us:** Of all patients we predicted would die, what percentage actually died?

**When it matters most:**

-   When **false positives are costly**
-   Example: Spam detection (don't want to lose important emails)
-   In our case: Unnecessary anxiety/treatment for patients incorrectly flagged as high-risk

------------------------------------------------------------------------

### Recall (Sensitivity) - "Of all actual positives, did we find them?"

$$
\text{Recall} = \frac{TP}{TP + FN} = \frac{TP}{\text{all actual positives}}
$$

**What it tells us:** Of all patients who actually died, what percentage did we correctly identify?

**When it matters most:**

-   When **false negatives are costly**
-   Example: Disease detection (don't want to miss cases!)
-   In our case: **Missing a death is worse than a false alarm**
-   this is our priority!

------------------------------------------------------------------------

### F1-Score - "Give me one number that balances both"

$$
\text{F1} = 2 \times \frac{\text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}}
$$

**What it is:** The harmonic mean of precision and recall

**When to use it:** - When you need to balance both precision and recall - When both false positives AND false negatives matter - Good single metric for imbalanced data

**Note:** F1 = 0 if either precision or recall is 0 (can't hide bad performance on one!)

------------------------------------------------------------------------

### Balanced Accuracy - "Accuracy that treats both classes equally"

$$
\text{Balanced Accuracy} = \frac{\text{Sensitivity} + \text{Specificity}}{2}
$$

Where **Specificity** = $\frac{TN}{TN + FP}$ (proportion of negatives correctly identified)

**Why it's better than regular accuracy:**

-   Treats both classes equally regardless of their frequency
-   For imbalanced data, gives a more honest assessment
-   If model only predicts majority class: Balanced Acc ≈ 0.5 (reveals the problem!)

```{r compare_accuracy_types}
# Compare regular accuracy vs balanced accuracy for our "stupid" baseline
stupid_bal_acc <- bal_accuracy(stupid_baseline, truth = died, estimate = .pred_class, 
                               event_level = "second")$.estimate

cat("'Always predict Survived' baseline:\n")
cat(sprintf("  Regular accuracy: %.1f%%\n", stupid_acc * 100))
cat(sprintf("  Balanced accuracy: %.1f%%\n", stupid_bal_acc * 100))
cat("\nBalanced accuracy reveals the model is useless (≈50% = random guessing)!\n")
```

------------------------------------------------------------------------

### :bar_chart: What our metrics tell us about our actual model

| Metric | Value | Meaning |
|-------------------|-------------------|----------------------------------|
| **Precision** | 0.80 | When the model predicts *Died*, it’s right 80% of the time. |
| **Recall (Sensitivity)** | 0.50 | It only identifies 50% of the actual deaths — **misses half of them**. |
| **F₁ Score** | 0.62 | The balance between precision and recall — reflects that recall is low. |
| **Balanced Accuracy** | 0.72 | Average of sensitivity (0.50) and specificity (≈0.94) — accounts for class imbalance. |

#### :speech_balloon: Interpretation

-   The model is **cautious** in predicting deaths — it only calls “death” when it’s fairly confident (high precision).
-   However, that caution means it **fails to flag many true deaths** (low recall).
-   The **F₁ score** summarizes this trade-off: decent precision, weak recall → moderate overall performance.
-   The **balanced accuracy (72%)** is lower than the raw accuracy (80%) because it accounts for the imbalance — it penalizes the model for missing so many deaths.

#### :brain: Connecting back to the confusion matrix

| From the matrix | Metric connection |
|-------------------------------------------|-----------------------------|
| 12 false negatives (actual deaths predicted as survived) | → hurts **recall** |
| 3 false positives (survivors predicted as died) | → slightly lowers **precision** |
| 48 true negatives + 12 true positives | → together drive **balanced accuracy** |

So, the confusion matrix **explains the story behind the numbers**:

> The model is “playing it safe” — rarely cries wolf (high precision), but sometimes fails to bark when it should (low recall).

------------------------------------------------------------------------

### The Precision-Recall Tradeoff

There's an inherent tension between precision and recall:

```{r precision_recall_tradeoff, fig.width=10, fig.height=5}
# Show how changing threshold affects precision and recall
thresholds <- seq(0.1, 0.9, by = 0.05)

# Get probability predictions from the fitted logistic model
simple_pred_prob <- predict(simple_fit, heart_test, type = "prob") %>%
  bind_cols(heart_test %>% select(died))

# Compute precision, recall, and F1 for each threshold value
threshold_results <- map_df(thresholds, function(thresh) {
  pred_at_threshold <- simple_pred_prob %>%
    mutate(
      .pred_class = factor(
        if_else(.pred_Died >= thresh, "Died", "Survived"),
        levels = c("Survived", "Died")
      )
    )
  
  tibble(
    threshold = thresh,
    precision = precision_vec(pred_at_threshold$died, pred_at_threshold$.pred_class, 
                              event_level = "second"),
    recall = recall_vec(pred_at_threshold$died, pred_at_threshold$.pred_class, 
                        event_level = "second"),
    f1 = f_meas_vec(pred_at_threshold$died, pred_at_threshold$.pred_class, 
                    event_level = "second")
  )
})

# --- Plot with shaded regions highlighting the precision–recall tradeoff ---
threshold_results %>%
  # Convert metrics to long format so ggplot can map color to "metric"
  pivot_longer(-threshold, names_to = "metric", values_to = "value") %>%
  ggplot(aes(x = threshold, y = value, color = metric)) +
  
  # --- Shaded regions to visually emphasize tradeoff extremes ---
  # Low thresholds → liberal model → high recall, low precision
  annotate("rect", xmin = 0.1, xmax = 0.4, ymin = 0, ymax = 1,
           alpha = 0.1, fill = "steelblue") +
  # High thresholds → conservative model → high precision, low recall
  annotate("rect", xmin = 0.6, xmax = 0.9, ymin = 0, ymax = 1,
           alpha = 0.1, fill = "darkgreen") +
  
  # --- Metric curves ---
  geom_line(linewidth = 1.2) +
  
  # --- Default 0.5 threshold line for reference ---
  geom_vline(xintercept = 0.5, linetype = "dashed", alpha = 0.6) +
  
  # --- Text labels for regions ---
  annotate("text", x = 0.25, y = 0.15, 
           label = "High recall\nLow precision", 
           color = "blue", size = 3.5, fontface = "italic") +
  annotate("text", x = 0.75, y = 0.15, 
           label = "High precision\nLow recall", 
           color = "darkgreen", size = 3.5, fontface = "italic") +
  
  # --- Titles, axis labels, legend ---
  labs(
    title = "Precision–Recall Tradeoff",
    subtitle = "Shaded regions highlight threshold extremes",
    x = "Classification Threshold",
    y = "Score",
    color = "Metric"
  ) +
  
  # --- Clean minimalist theme ---
  theme_minimal(base_size = 11)
```

#### :compass: Reading the plot

-   **X-axis:** the classification threshold — how confident the model must be before predicting *“Died.”*
-   **Y-axis:** the value of each metric (precision, recall, F₁).
-   **Vertical dashed line:** the **default threshold (0.5)** used by logistic regression.
-   **Colored lines:** how each metric changes as you make the model more or less “cautious.”

#### :bar_chart: What the lines tell you

| Trend | What it means |
|----------------------------|--------------------------------------------|
| **Blue line (recall)** falls as the threshold increases | As you require higher confidence to predict *“Died,”* you miss more true deaths. |
| **Green line (precision)** rises as the threshold increases | The model makes fewer false alarms — it’s more confident when it does predict *“Died.”* |
| **Red line (F₁)** peaks at an intermediate threshold | That’s the “sweet spot” balancing catching deaths (recall) and avoiding false alarms (precision). |

#### Interpretation

> When the threshold is **low**, the model calls “death” too easily — it catches most deaths (high recall) but makes lots of false alarms (low precision). When the threshold is **high**, the model calls “death” only when very sure — precision rises, but recall collapses. Around **0.4–0.5**, both are moderate, and the **F₁ score** reaches its maximum — this is usually the most balanced choice.

**For our heart failure problem:** We probably want to prioritize **recall** (don't miss deaths) even if it means more false alarms!

------------------------------------------------------------------------

## 3.2 :new: ROC and PR Curves - Visualizing Performance

### ROC Curve (Receiver Operating Characteristic)

The **ROC curve** plots True Positive Rate (TPR = Recall) vs False Positive Rate (FPR = 1 - Specificity) across all possible classification thresholds.

```{r roc_curve}
# ROC curve
roc_data <- roc_curve(simple_pred_prob, truth = died, .pred_Died, 
                      event_level = "second")

autoplot(roc_data) +
  labs(title = "ROC Curve",
       subtitle = "Closer to top-left corner = better model") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray50") +
  annotate("text", x = 0.6, y = 0.4, 
           label = "Random classifier\n(diagonal line)",
           color = "gray50", size = 3)

# Calculate AUC
roc_auc_value <- roc_auc(simple_pred_prob, truth = died, .pred_Died, 
                         event_level = "second")
cat(sprintf("ROC-AUC: %.3f\n", roc_auc_value$.estimate))
```

#### :compass: Interpreting the ROC curve

-   The **x-axis (1 − specificity)** measures the **false positive rate** — how often the model incorrectly predicts “Died” when the patient actually survived.
-   The **y-axis (sensitivity)** measures the **true positive rate** — how often the model correctly identifies patients who died.

Each point on the curve corresponds to a **different threshold** for classifying someone as “Died.” By varying that threshold, we trace how the tradeoff between *sensitivity* and *specificity* changes.

#### :bulb: Reading the plot visually

| Region | Meaning |
|-------------------------|-----------------------------------------------|
| **Top-left corner (ideal)** | Model catches nearly all deaths (high sensitivity) while rarely flagging survivors (low false positive rate). |
| **Diagonal line** | Random guessing — sensitivity rises only in proportion to false positives. |
| **Curve above diagonal** | Model performs better than chance — the farther above, the better. |

Our curve hugs the upper-left corner quite closely, showing the model is doing a **good job separating deaths from survivors** across all thresholds.

-   **AUC = 0.85** means there’s an **85% chance** the model will rank a randomly chosen patient who died **higher** (i.e., more likely to die) than a randomly chosen patient who survived.
-   AUC values:
    -   0.5 → random guessing (diagonal line)

    -   0.7–0.8 → fair

    -   0.8–0.9 → **good**

    -   

        > 0.9 → excellent (but check for overfitting!)

So, this model has **good discriminative power**, meaning it’s quite effective at distinguishing high-risk patients.

> The precision–recall curve focuses on performance for the *positive class* (especially useful for imbalanced data), while the ROC curve summarizes performance **across all thresholds** using both classes. Both visualize the same principle — threshold tuning — but from different perspectives.

::: callout-warning
Check out [this Medium post](https://medium.com/data-science/imbalanced-data-stop-using-roc-auc-and-use-auprc-instead-46af4910a494) and [this one](https://juandelacalle.medium.com/how-and-why-i-switched-from-the-roc-curve-to-the-precision-recall-curve-to-analyze-my-imbalanced-6171da91c6b8) to see why ROC curves and AUC are not necessarily suitable for data with class imbalance. You can also check [this post](https://towardsdatascience.com/precision-recall-curve-is-more-informative-than-roc-in-imbalanced-data-4c95250242f6/)
:::

------------------------------------------------------------------------

## :brain: Why ROC curves look too good on imbalanced data

While ROC curves are widely used to assess classification performance, they can appear **overly optimistic on imbalanced datasets**.

This happens because the ROC curve’s *False Positive Rate (FPR)* = FP / (TN + FP) is **artificially low** when negatives dominate — the large number of True Negatives (TN) in the denominator makes the FPR tiny, even if the model produces many false positives relative to the few actual positives.

As a result, the ROC curve can look excellent even when the model performs poorly on the minority (positive) class.

> Let’s recall the definitions: $$
> \text{True Negative Rate (TNR)} = \frac{TN}{TN + FP}, \quad
> \text{False Positive Rate (FPR)} = 1 - \text{TNR} = \frac{FP}{TN + FP}
> $$ Now suppose negatives vastly outnumber positives — e.g., 990 negatives and 10 positives.

> Even if your model **misclassifies many of the 10 positives** (bad for Recall), the denominator of FPR — ( TN + FP ) — is dominated by **the huge number of true negatives (TN)**. That means:

> -   ( FP ) is small relative to ( TN + FP )
> -   → ( FPR = FP / (TN + FP) ) is tiny
> -   → ( TNR = TN / (TN + FP) ) is very high

> So the ROC curve — which uses TPR vs. FPR — shows a curve that *hugs the top-left corner*, suggesting great performance, even when **the model barely detects positives**.

**Precision–Recall (PR) curves**, by contrast, focus only on the positive class — plotting **Precision** (how many predicted positives are correct) against **Recall** (how many true positives are found). This makes PR curves **much more informative** for evaluating performance on **imbalanced data** such as fraud detection or rare disease prediction.

```{r pr_curve_new}
# Calculate PR curve
pr_data <- pr_curve(simple_pred_prob, truth = died, .pred_Died, 
                    event_level = "second")

# Plot PR curve
baseline_precision <- sum(heart_test$died == "Died") / nrow(heart_test)

autoplot(pr_data) +
  geom_hline(yintercept = baseline_precision,
             linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = 0.5, y = baseline_precision + 0.05,
           label = sprintf("Random baseline\n(%.1f%% prevalence)", 
                          baseline_precision * 100),
           color = "red", size = 3) +
  labs(title = "Precision-Recall Curve",
       subtitle = "Better than ROC for imbalanced data - focuses on minority class") +
  theme_minimal() 

# Calculate PR-AUC
pr_auc_value <- pr_auc(simple_pred_prob, truth = died, .pred_Died, 
                       event_level = "second")
cat(sprintf("PR-AUC: %.3f\n", pr_auc_value$.estimate))
cat(sprintf("Baseline (random): %.3f\n", baseline_precision))
```

**Interpreting our results**

PR-AUC: 0.760 Baseline (random): 0.320 (the positive class prevalence)

*What this means:*

-   A random classifier (pure guessing) would have a flat PR curve at the (minority positive class) prevalence line (here, 32%).

-   Our model’s PR-AUC = 0.76 is well above baseline, showing that:

    -   It identifies positives much better than random chance.

    -   Its predictions maintain high precision even as recall increases.

-   The curve’s early high precision (near 1.0) means the top-ranked predictions are highly reliable — few false positives.

-   As recall rises, precision drops, meaning that to find more true positives, the model starts including more false alarms (typical trade-off).

**Why PR is better for imbalanced data:**

1.  **Focuses on minority class** (positive class) - what we care about!
2.  **More sensitive to improvements** in minority class performance
3.  **Baseline is class prevalence** (not 0.5 like ROC)
4.  **Shows true cost of errors** in the positive class
5.  **Can't hide poor performance** with lots of true negatives

**When to use which:**

-   **Balanced data**: ROC curve is fine
-   **Imbalanced data**: PR curve is more informative and honest
-   **Best practice**: Plot both for complete picture!

**Comparison:**

| Aspect                     | ROC Curve            | PR Curve                  |
|--------------------------|---------------------|-------------------------|
| Focus                      | Both classes equally | Minority (positive) class |
| Baseline                   | 0.5 (diagonal)       | Class prevalence          |
| Best for                   | Balanced data        | Imbalanced data           |
| Can hide poor performance? | Yes (with many TNs)  | No                        |

------------------------------------------------------------------------

## 3.3 :new: Calibration Curves - Are Probabilities Reliable?

**Calibration** measures whether predicted probabilities match actual frequencies.

**A well-calibrated model**: If the model predicts 70% probability of death for 100 patients, approximately 70 of them should actually die.

```{r calibration_curve}
# Create calibration data
calibration_data <- simple_pred_prob %>%
  mutate(
    prob_bin = cut(.pred_Died, 
                   breaks = seq(0, 1, 0.1),
                   include.lowest = TRUE),
    actual_died = as.numeric(died == "Died")
  ) %>%
  group_by(prob_bin) %>%
  summarise(
    mean_predicted = mean(.pred_Died),
    mean_actual = mean(actual_died),
    n = n(),
    .groups = "drop"
  ) %>%
  filter(n > 0)  # Only bins with data

# Plot calibration curve
ggplot(calibration_data, aes(x = mean_predicted, y = mean_actual)) +
  geom_point(aes(size = n), alpha = 0.7, color = "steelblue") +
  geom_line(color = "steelblue") +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  annotate("text", x = 0.7, y = 0.3, 
           label = "Perfect calibration\n(diagonal line)", 
           color = "red", hjust = 0) +
  labs(title = "Calibration Curve",
       subtitle = "Do predicted probabilities match actual frequencies?",
       x = "Predicted Probability of Death",
       y = "Actual Proportion Who Died",
       size = "# Patients in bin") +
  coord_fixed() +
  xlim(0, 1) + ylim(0, 1) +
  theme_minimal()
```

**How to read the calibration curve:**

The calibration curve compares predicted probabilities to observed outcomes.

A perfectly calibrated model would have all points on the red diagonal line — for example, among patients with a predicted 30% probability of death, about 30% would actually die.

In this model, most points lie close to the diagonal, suggesting the probabilities are reasonably accurate. Some deviations around 0.4–0.6 indicate mild underestimation, meaning the model predicts lower risk than actually observed for those patients.

Calibration complements discrimination metrics (like ROC or PR curves) by assessing whether the predicted probabilities mean what they say.

**Why calibration matters:**

-   In healthcare, we use probabilities to make decisions (e.g., "treat if risk \> 60%")
-   Miscalibrated probabilities lead to poor decisions
-   A model can have high AUC but poor calibration!
-   Well-calibrated models → trustworthy risk predictions for individual patients

**Example:** If your model is well-calibrated and predicts a patient has 80% risk of death, you can confidently say "4 out of 5 similar patients would die" - this is crucial for informed decision-making!

------------------------------------------------------------------------

# Part 4: K-Nearest Neighbors (K-NN) Classification

## 4.1 :new: The K-NN Algorithm

**K-Nearest Neighbors (K-NN)** is a simple, intuitive classification algorithm:

**Basic idea**:

1\. To classify a new data point, look at its $K$ nearest neighbors in the training data

2\. Count how many neighbors belong to each class

3\. Assign the new point to the majority class

**Key hyperparameter:** $K$

-   $K$ = number of neighbors to consider
-   $K$ = 1: Very flexible, prone to overfitting
-   $K$ = large: Very smooth, may underfit
-   $K$ is typically odd (to avoid ties)

**Quick example (images and example from [here](https://www.analyticsvidhya.com/blog/2018/03/introduction-k-neighbours-algorithm-clustering/))**

Look at this example.

![](images/clipboard-3976936043.png)

Say you want to predict the class of the blue star (BS) and see whether it's a red circle (RC) or a green square (GS).

![](images/clipboard-3904064819.png)

$K$-NN will take "votes" from the nearest neighbours (nearest defined based on pre-specified distance usually Euclidian distance) and will only look at $K$ neighbours (hence the $K$ in $K$-NN).

In this case, if we take $K=3$, we see that the blue star is surrounded by 3 red circles and no green squares, so the star would be classified as a red circle.

### KNN on the heart failure example

```{r knn_illustration, fig.width=10, fig.height=5}
set.seed(456)

# ------------------------------------------------------------------------------
# Data subset for visualization
# ------------------------------------------------------------------------------
viz_data <- heart_train %>%
  select(ejection_fraction, serum_creatinine, died) %>%
  slice_sample(n = 100)

# Test point to classify
test_point <- tibble(
  ejection_fraction = 35,
  serum_creatinine = 1.5
)

# ------------------------------------------------------------------------------
# Recipe + model specs (consistent with tidymodels style)
# ------------------------------------------------------------------------------
viz_recipe <- recipe(died ~ ejection_fraction + serum_creatinine, data = viz_data) %>%
  step_normalize(all_numeric_predictors())

knn_spec_k3 <- nearest_neighbor(neighbors = 3) %>%
  set_engine("kknn") %>%
  set_mode("classification")

knn_spec_k7 <- nearest_neighbor(neighbors = 7) %>%
  set_engine("kknn") %>%
  set_mode("classification")

# ------------------------------------------------------------------------------
# Fit models using workflows
# ------------------------------------------------------------------------------
knn_fit_k3 <- workflow() %>%
  add_recipe(viz_recipe) %>%
  add_model(knn_spec_k3) %>%
  fit(data = viz_data)

knn_fit_k7 <- workflow() %>%
  add_recipe(viz_recipe) %>%
  add_model(knn_spec_k7) %>%
  fit(data = viz_data)

# ------------------------------------------------------------------------------
# Generate grid for visualization
# ------------------------------------------------------------------------------
grid <- expand.grid(
  ejection_fraction = seq(min(viz_data$ejection_fraction), max(viz_data$ejection_fraction), length.out = 100),
  serum_creatinine  = seq(min(viz_data$serum_creatinine), max(viz_data$serum_creatinine), length.out = 100)
)

# Use augment() to get predictions for each grid point
grid_pred_k3 <- augment(knn_fit_k3, new_data = grid)
grid_pred_k7 <- augment(knn_fit_k7, new_data = grid)

# ------------------------------------------------------------------------------
# Visualization: decision regions + sample points + test point
# ------------------------------------------------------------------------------
p1 <- ggplot() +
  geom_tile(data = grid_pred_k3,
            aes(x = ejection_fraction, y = serum_creatinine, fill = .pred_class),
            alpha = 0.25) +
  geom_point(data = viz_data,
             aes(x = ejection_fraction, y = serum_creatinine, color = died),
             size = 3) +
  geom_point(data = test_point,
             aes(x = ejection_fraction, y = serum_creatinine),
             shape = 8, size = 5, color = "black") +
  labs(
    title = "KNN with K = 3",
    subtitle = "Star = new test point; shaded regions = predicted class",
    x = "Ejection Fraction",
    y = "Serum Creatinine"
  ) +
  scale_fill_manual(values = c("Lived" = "steelblue", "Died" = "coral")) +
  theme_minimal() +
  guides(fill = "none")

p2 <- ggplot() +
  geom_tile(data = grid_pred_k7,
            aes(x = ejection_fraction, y = serum_creatinine, fill = .pred_class),
            alpha = 0.25) +
  geom_point(data = viz_data,
             aes(x = ejection_fraction, y = serum_creatinine, color = died),
             size = 3) +
  geom_point(data = test_point,
             aes(x = ejection_fraction, y = serum_creatinine),
             shape = 8, size = 5, color = "black") +
  labs(
    title = "KNN with K = 7",
    subtitle = "Larger K = smoother decision boundaries",
    x = "Ejection Fraction",
    y = "Serum Creatinine"
  ) +
  scale_fill_manual(values = c("Lived" = "steelblue", "Died" = "coral")) +
  theme_minimal() +
  guides(fill = "none")

p1 + p2
```

Perfect — that plot looks great! :dart: You’ve successfully used **tidymodels** to visualize **K-Nearest Neighbors (KNN) decision boundaries** with different values of ( K ) (3 vs. 7).

Let’s briefly interpret your results, point out a few improvements you could make, and explain how to take it further (e.g., showing class probabilities or CV tuning).

------------------------------------------------------------------------

## :brain: Interpretation of your current figure

| Plot | What it shows | Key takeaway |
|-------------------|---------------------------------|---------------------|
| **Left (K = 3)** | Highly detailed boundaries — the shaded “Died” (gray) vs “Survived” (red) regions change rapidly around training points. | Model is **very flexible**, capturing local noise; prone to **overfitting**. The decision boundaries are more irregular: the model closely follows the training data. |
| **Right (K = 7)** | Boundaries are smoother and more stable. | Larger ( K ) → **more generalization**, less sensitivity to outliers. |

The **black star** marks your *test point* (EF = 35, SC = 1.5). Its surrounding region determines how KNN would classify it.

→ If it’s in the gray area, predicted = **“Died”**; if red, **“Survived.”**

**The K-NN Algorithm Steps**:

1.  **Choose K** (number of neighbors)
2.  **Calculate distance** between new point and all training points
    -   Usually Euclidean distance: $d = \sqrt{\sum_{i=1}^p (x_i - y_i)^2}$
3.  **Find K nearest neighbors** (smallest distances)
4.  **Take majority vote** among these K neighbors
5.  **Assign class** based on majority

------------------------------------------------------------------------

## 4.2 :new: Fitting K-NN in tidymodels (review)

**Important preprocessing for K-NN**:

-   K-NN uses **distance**, so scale matters!
-   **Always normalize/standardize features** before applying K-NN
-   Otherwise, features with large ranges dominate

```{r knn_fit}
# Create recipe with normalization (CRITICAL for K-NN!)
knn_recipe <- recipe(died ~ ejection_fraction + serum_creatinine + age + 
                            serum_sodium + time, 
                     data = heart_train) %>%
  step_normalize(all_numeric_predictors())  # Essential!

# Specify K-NN model with K=5
knn_spec <- nearest_neighbor(neighbors = 5) %>%
  set_engine("kknn") %>%
  set_mode("classification")

# Create workflow
knn_workflow <- workflow() %>%
  add_recipe(knn_recipe) %>%
  add_model(knn_spec)

# Fit
knn_fit <- fit(knn_workflow, data = heart_train)

cat("K-NN model with K=5 trained successfully!\n")
```

------------------------------------------------------------------------

## 4.3 :new: Evaluating K-NN

```{r knn_evaluate}
# Get predictions (both class and prob)
log_pred <- predict(unreg_fit, heart_test, type = "prob") %>%
  bind_cols(predict(unreg_fit, heart_test, type = "class")) %>%
  bind_cols(heart_test %>% select(died))

# Define two metric sets
label_metrics <- metric_set(bal_accuracy, yardstick::precision, yardstick::recall, f_meas)
prob_metrics  <- metric_set(roc_auc, pr_auc)

# Compute class-label metrics
log_label_metrics <- log_pred %>%
  label_metrics(truth = died, estimate = .pred_class, event_level = "second")

# Compute probability-based metrics
log_prob_metrics <- log_pred %>%
  prob_metrics(truth = died, .pred_Died, event_level = "second")

# Combine both
log_metrics <- bind_rows(log_label_metrics, log_prob_metrics) %>%
  mutate(model = "Logistic Regression")

# Do the same for KNN
knn_label_metrics <- knn_pred %>%
  label_metrics(truth = died, estimate = .pred_class, event_level = "second")

knn_prob_metrics <- knn_pred %>%
  prob_metrics(truth = died, .pred_Died, event_level = "second")

knn_metrics <- bind_rows(knn_label_metrics, knn_prob_metrics) %>%
  mutate(model = "K-NN (K=5)")

# Combine for comparison
comparison <- bind_rows(log_metrics, knn_metrics) %>%
  select(model, .metric, .estimate) %>%
  pivot_wider(names_from = .metric, values_from = .estimate)

# Display
cat("\nComparison:\n")
print(comparison)
```

#### :compass: Interpretation

| Metric | Logistic Regression | KNN (K=5) | What it tells us |
|:----------------|:----------------|:----------------|:-------------------|
| **Balanced Accuracy** | **0.888** | 0.794 | Logistic regression correctly identifies both classes more evenly (good for imbalance). |
| **Precision** | 0.808 | 0.800 | Both models are similar — when they predict “Died,” they’re correct \~80% of the time. |
| **Recall (Sensitivity)** | **0.875** | 0.667 | Logistic regression finds a higher proportion of actual deaths (fewer false negatives). |
| **F1-score** | **0.84** | 0.73 | Logistic regression balances precision and recall better. |
| **ROC-AUC** | **0.936** | 0.880 | Logistic regression has stronger class separation overall. |
| **PR-AUC** | **0.874** | 0.754 | Logistic regression maintains better precision at all recall levels, especially important with imbalance. |

**Summary:**

> Logistic Regression outperforms KNN across all key metrics, especially Recall, F1, and PR-AUC — meaning it’s better at detecting the minority “Died” class while keeping false positives low.

KNN still performs decently, but its lower recall and PR-AUC show it’s less reliable when positives are rare.

```{r}
comparison_long <- comparison %>%
  pivot_longer(cols = -model, names_to = "metric", values_to = "value")

ggplot(comparison_long, aes(x = metric, y = value, fill = model)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.6) +
  coord_flip() +
  geom_text(aes(label = round(value, 3)), 
            position = position_dodge(width = 0.7), 
            hjust = -0.2, size = 3.2) +
  labs(
    title = "Model Comparison on Imbalanced Classification Metrics",
    subtitle = "Balanced Accuracy, Precision, Recall, F1, ROC-AUC, and PR-AUC",
    x = NULL, y = "Metric Value", fill = "Model"
  ) +
  ylim(0, 1.1) +
  theme_minimal(base_size = 12)
```

```{r}
# Helper function for calibration data
make_calibration <- function(pred_df, model_name) {
  pred_df %>%
    mutate(bin = ntile(.pred_Died, 10)) %>%  # deciles of predicted probability
    group_by(bin) %>%
    summarise(
      mean_pred = mean(.pred_Died),
      obs_rate  = mean(died == "Died"),
      n = n()
    ) %>%
    mutate(model = model_name)
}

# Get calibration data for both models
cal_data <- bind_rows(
  make_calibration(log_pred, "Logistic Regression"),
  make_calibration(knn_pred, "K-NN (K=5)")
)

# Plot calibration curves
ggplot(cal_data, aes(x = mean_pred, y = obs_rate, color = model)) +
  geom_point(size = 3) +
  geom_line(linewidth = 1) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "gray40") +
  labs(
    title = "Calibration Curves",
    subtitle = "Perfect calibration = dashed diagonal line",
    x = "Predicted Probability of Death",
    y = "Observed Proportion Who Died",
    color = "Model"
  ) +
  theme_minimal(base_size = 12)
```

#### What the calibration curves show

##### :small_blue_diamond: Logistic Regression (teal)

-   Closely tracks the **dashed diagonal** — meaning predicted probabilities correspond closely to actual observed rates.
-   For example, among patients predicted to have a 0.6 risk of death, about 60% actually died.
-   **Conclusion:** Logistic regression is *well-calibrated* — its probabilities are reliable as risk estimates.

##### :small_blue_diamond: KNN (red)

-   Deviates more from the diagonal — especially at mid-to-high probabilities.
-   Sometimes **overestimates risk** (curve above diagonal), sometimes **underestimates** (below).
-   This is typical for KNN: it produces more "step-like" probability behavior that can overfit local structure.
-   **Conclusion:** KNN has reasonable discrimination (AUC) but poorer calibration.

> *"Both models can classify patients fairly well, but only logistic regression produces probabilities you can actually trust as predicted risks. KNN might say ‘70% risk of death,’ but in reality, only \~50–60% of such patients die — it’s not well calibrated."*

> To summarize: → **AUC ≠ Calibration.** → Calibration reflects *probability reliability*, not *ranking ability*.

------------------------------------------------------------------------

## 4.4 :new: The Effect of K (Bias-Variance Tradeoff)

$K$ controls the complexity of the model:

```{r k_effect, fig.width=10, fig.height=6}
 # Define values of K to try
k_values <- c(1, 3, 5, 7, 11, 15, 21, 31, 51)

# Evaluate each K
k_results <- map_df(k_values, function(k) {
  
  # Define model
  knn_k_spec <- nearest_neighbor(neighbors = k) %>%
    set_engine("kknn") %>%
    set_mode("classification")
  
  knn_k_workflow <- workflow() %>%
    add_recipe(knn_recipe) %>%
    add_model(knn_k_spec)
  
  # Fit on training data
  knn_k_fit <- fit(knn_k_workflow, data = heart_train)
  
  # --- TRAIN performance ---
  train_pred <- predict(knn_k_fit, heart_train, type = "prob") %>%
    bind_cols(predict(knn_k_fit, heart_train, type = "class")) %>%
    bind_cols(heart_train %>% select(died))
  
  train_f1 <- f_meas(train_pred, truth = died, estimate = .pred_class,
                     event_level = "second")$.estimate
  
  # --- TEST performance ---
  test_pred <- predict(knn_k_fit, heart_test, type = "prob") %>%
    bind_cols(predict(knn_k_fit, heart_test, type = "class")) %>%
    bind_cols(heart_test %>% select(died))
  
  test_f1 <- f_meas(test_pred, truth = died, estimate = .pred_class,
                    event_level = "second")$.estimate
  
  # Return results
  tibble(
    K = k,
    train_f1 = train_f1,
    test_f1 = test_f1
  )
})

# --- Visualization ---
k_results %>%
  pivot_longer(cols = starts_with("train_") | starts_with("test_"),
               names_to = "set", values_to = "F1") %>%
  ggplot(aes(x = K, y = F1, color = set)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  scale_color_manual(values = c("train_f1" = "steelblue", "test_f1" = "coral"),
                     labels = c("Training", "Test")) +
  labs(
    title = "Effect of K on K-NN F1-score",
    subtitle = "Small K = overfitting (high train, low test); large K = underfitting",
    x = "K (Number of Neighbors)",
    y = "F1-score",
    color = "Dataset"
  ) +
  theme_minimal(base_size = 12)

# --- Find optimal K based on test F1 ---
optimal_k <- k_results %>%
  arrange(desc(test_f1)) %>%
  slice(1)

cat(sprintf("\nOptimal K based on test F1-score: %d\n", optimal_k$K))
cat(sprintf("Test F1-score: %.3f\n", optimal_k$test_f1))
```

> “For imbalanced data, F1 is a better optimization target than accuracy — it rewards models that maintain both precision and recall.”

**Key observations:**

-   **Small K (e.g., K = 1–3)**

    -   Very high training F1 (memorizes training data → *overfitting*)
    -   Much lower test F1 (poor generalization)
    -   **High variance**, **low bias**: predictions depend heavily on local noise

-   **Large K (e.g., K ≥ 30)**

    -   Training and test F1 both decrease — the model becomes too smooth
    -   Fails to capture local structure or rare patterns
    -   **Low variance**, **high bias**: underfitting due to excessive averaging

-   **Optimal K**

    -   Found at the **balance point** between bias and variance
    -   Typically chosen via **cross-validation** (covered next!)
    -   Gives stable, generalizable performance on unseen data

**This motivates cross-validation**: We need a principled way to choose K!

------------------------------------------------------------------------

# Part 5: :new: Cross-Validation for Hyperparameter Tuning

## 5.1 The Problem with Single Train-Test Split

So far, we've used a single train-test split. But this has problems:

1.  **Randomness**: Results depend on which data points end up in test set
2.  **Data waste**: Test set isn't used for training
3.  **No hyperparameter tuning**: How do we choose $K$ for K-NN? Or penalty for regularization?

**The solution: Cross-validation!**

------------------------------------------------------------------------

## 5.2 K-Fold Cross-Validation

**K-Fold CV** (don't confuse with K in K-NN!):

1.  Split data into **K equal-sized folds** (typically K=5 or K=10)
2.  For each fold:
    -   Train on K-1 folds
    -   Test on remaining fold
3.  Average performance across all K folds
4.  All data gets used for both training and testing!

![](https://upload.wikimedia.org/wikipedia/commons/4/4b/KfoldCV.gif){fig-alt="Animated illustration of K-fold cross-validation" width="70%"}

```{r cv_explanation, fig.width=10, fig.height=6}
# Visual representation of 5-fold CV
cv_viz <- tibble(
  fold = rep(1:5, each = 5),
  iteration = rep(1:5, 5),
  role = case_when(
    fold == iteration ~ "Test",
    TRUE ~ "Train"
  )
)

ggplot(cv_viz, aes(x = fold, y = iteration, fill = role)) +
  geom_tile(color = "white", linewidth = 1) +
  scale_fill_manual(values = c("Train" = "steelblue", "Test" = "coral")) +
  labs(title = "5-Fold Cross-Validation",
       subtitle = "Each fold serves as test set exactly once",
       x = "Fold",
       y = "CV Iteration",
       fill = "Role") +
  theme_minimal() +
  theme(panel.grid = element_blank(),
        text = element_text(size = 12)) %>%
  print()
```

------------------------------------------------------------------------

## 5.3 Using CV to Tune K in K-NN

Let's use cross-validation to find the optimal K!

```{r cv_tune_k}
set.seed(789)

# ------------------------------------------------------------------------------
# 1. Create stratified 5-fold CV splits
# ------------------------------------------------------------------------------

cv_folds <- vfold_cv(heart_train, v = 5, strata = died)

cat("Created 5 CV folds with", nrow(heart_train), "training observations\n")
cat("Each fold has ~", round(nrow(heart_train) / 5), "observations\n")

# ------------------------------------------------------------------------------
# 2. Define tunable KNN model
# ------------------------------------------------------------------------------

knn_tune_spec <- nearest_neighbor(neighbors = tune()) %>%
  set_engine("kknn") %>%
  set_mode("classification")

knn_tune_workflow <- workflow() %>%
  add_recipe(knn_recipe) %>%
  add_model(knn_tune_spec)

# ------------------------------------------------------------------------------
# 3. Define grid of K values
# ------------------------------------------------------------------------------

k_grid <- tibble(neighbors = c(1, 3, 5, 7, 11, 15, 21, 31, 51, 71))

# ------------------------------------------------------------------------------
# 4. Define metrics (F1 and a few reference metrics)
# ------------------------------------------------------------------------------

f1_metrics <- metric_set(f_meas, yardstick::precision, yardstick::recall, bal_accuracy, pr_auc, roc_auc)

# ------------------------------------------------------------------------------
# 5. Tune across CV folds
# ------------------------------------------------------------------------------

knn_tuning <- tune_grid(
  knn_tune_workflow,
  resamples = cv_folds,
  grid = k_grid,
  metrics = f1_metrics
)

# ------------------------------------------------------------------------------
# 6. Review results
# ------------------------------------------------------------------------------

cv_results <- collect_metrics(knn_tuning)
print(cv_results)

# ------------------------------------------------------------------------------
# 7. Visualize cross-validation results
# ------------------------------------------------------------------------------

cv_results %>%
  ggplot(aes(x = neighbors, y = mean, color = .metric)) +
  geom_line(linewidth = 1.1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = mean - std_err, ymax = mean + std_err),
                width = 2, alpha = 0.5) +
  facet_wrap(~.metric, scales = "free_y", ncol = 2) +
  labs(
    title = "Cross-Validation Results for K-NN (F1-focused)",
    subtitle = "Error bars = standard error across 5 folds",
    x = "K (Number of Neighbors)",
    y = "Mean Performance",
    color = "Metric"
  ) +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none")

# ------------------------------------------------------------------------------
# 8. Select best K (based on F1)
# ------------------------------------------------------------------------------

best_k <- select_best(knn_tuning, metric = "f_meas")
cat(sprintf("\nBest K based on F1-score: %d\n", best_k$neighbors))

# ------------------------------------------------------------------------------
# 9. Finalize and fit model on full training data
# ------------------------------------------------------------------------------

final_knn_workflow <- finalize_workflow(knn_tune_workflow, best_k)
final_knn_fit <- fit(final_knn_workflow, data = heart_train)

# ------------------------------------------------------------------------------
# 10. Evaluate final model on test set
# ------------------------------------------------------------------------------

final_knn_pred <- predict(final_knn_fit, heart_test, type = "prob") %>%
  bind_cols(predict(final_knn_fit, heart_test, type = "class")) %>%
  bind_cols(heart_test %>% select(died))

# Class-based metrics
label_metrics <- metric_set(f_meas, yardstick::precision, yardstick::recall, bal_accuracy)
# Probabilistic metrics
prob_metrics <- metric_set(pr_auc, roc_auc)

# Compute both
final_metrics <- bind_rows(
  label_metrics(final_knn_pred, truth = died, estimate = .pred_class, event_level = "second"),
  prob_metrics(final_knn_pred, truth = died, .pred_Died, event_level = "second")
)

cat("\nFinal K-NN (tuned) performance on test set:\n")
print(final_metrics)
```

**Benefits of cross-validation**:

1.  **More reliable estimates**: Average over multiple splits
2.  **Hyperparameter tuning**: Find optimal K (or penalty, etc.)
3.  **Detect overfitting**: Large gap between folds = unstable model
4.  **Uses all data**: No data wasted

------------------------------------------------------------------------

## 5.4 Detecting Overfitting with CV

```{r detect_overfitting}
# --- K=1 (likely overfit) ---
knn_k1_spec <- nearest_neighbor(neighbors = 1) %>%
  set_engine("kknn") %>%
  set_mode("classification")

knn_k1_workflow <- workflow() %>%
  add_recipe(knn_recipe) %>%
  add_model(knn_k1_spec)

# Fit and evaluate on training (F1)
knn_k1_fit <- fit(knn_k1_workflow, heart_train)
k1_train_pred <- predict(knn_k1_fit, heart_train, type = "class") %>%
  bind_cols(heart_train %>% select(died))
k1_train_f1 <- f_meas(k1_train_pred, truth = died, estimate = .pred_class, event_level = "second")$.estimate

# Get CV performance for K=1 (F1)
k1_cv_f1 <- cv_results %>%
  filter(.metric == "f_meas", neighbors == 1) %>%
  pull(mean)

# --- Optimal K ---
best_k_cv_f1 <- cv_results %>%
  filter(.metric == "f_meas", neighbors == best_k$neighbors) %>%
  pull(mean)

optk_train_pred <- predict(final_knn_fit, heart_train, type = "class") %>%
  bind_cols(heart_train %>% select(died))
optk_train_f1 <- f_meas(optk_train_pred, truth = died, estimate = .pred_class, event_level = "second")$.estimate

# --- Combine results ---
comparison_overfit <- tibble(
  K = c(1, 1, best_k$neighbors, best_k$neighbors),
  Set = c("Training", "CV (5-fold)", "Training", "CV (5-fold)"),
  F1 = c(k1_train_f1, k1_cv_f1, optk_train_f1, best_k_cv_f1)
)

# --- Plot ---
ggplot(comparison_overfit, aes(x = factor(K), y = F1, fill = Set)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = sprintf("%.3f", F1)),
            position = position_dodge(width = 0.9), vjust = -0.5) +
  labs(title = "Detecting Overfitting with Cross-Validation (F1-based)",
       subtitle = "Large gap between Training and CV = overfitting",
       x = "K",
       y = "F1-score") +
  scale_fill_manual(values = c("Training" = "steelblue", "CV (5-fold)" = "coral")) +
  ylim(0, 1) +
  theme_minimal(base_size = 12)

# --- Print summary ---
cat("\n📊 Overfitting Analysis (F1-based):\n")
cat(sprintf("K=1: Train=%.3f, CV=%.3f, Gap=%.3f (%.1f%%) → OVERFITTING!\n",
            k1_train_f1, k1_cv_f1, k1_train_f1 - k1_cv_f1,
            100 * (k1_train_f1 - k1_cv_f1) / k1_train_f1))
cat(sprintf("K=%d: Train=%.3f, CV=%.3f, Gap=%.3f (%.1f%%) → Good generalization\n",
            best_k$neighbors, optk_train_f1, best_k_cv_f1, 
            optk_train_f1 - best_k_cv_f1,
            100 * (optk_train_f1 - best_k_cv_f1) / optk_train_f1))

```

A large positive gap between training and test performance indicates overfitting and conversely, a large negative gap is a sign of underfitting.

#### :jigsaw: 1. Think conceptually first

| Size of | Pattern | Typical interpretation | Model behaviour |
|:------------------|:-----------------|:-----------------|:-----------------|
| **Very small gap** (\< ≈ 0.05 F1) | Train ≈ CV | Good generalization | Bias–variance in balance |
| **Moderate positive gap** (≈ 0.05 – 0.15) | Train \> CV | Mild overfitting | A bit too flexible |
| **Large positive gap** (\> ≈ 0.15) | Train ≫ CV | Strong overfitting | High variance, memorizing |
| **Moderate negative gap** (≈ −0.05 – −0.15) | CV \> Train | Slight underfitting | Slightly too simple |
| **Large negative gap** (\< −0.15 to −0.20) | CV ≫ Train | Clear underfitting | Too smooth / high bias |

> In short: **\|Gap\| \> 0.15 (≈ 15 percentage points)** is usually a *red flag* either way.

These thresholds aren’t strict rules — they’re guideposts that work well for most F1-based or accuracy-based evaluations.

#### :brain: 2. Context matters

-   **Imbalanced data** → expect *larger metric variance* across folds; a “large” gap might start around 0.10 rather than 0.15.
-   **Tiny datasets** → CV can fluctuate more, so interpret gaps qualitatively.
-   **High-variance metrics** (like F1 or PR-AUC) → use standard errors to judge whether the gap is statistically meaningful.

A quick rule of thumb:

> If the gap exceeds **twice the cross-validation standard error**, treat it as significant.

> -   **Small gap → good generalization**
> -   **Large positive gap → overfitting**
> -   **Large negative gap → underfitting**
> -   The “large” threshold depends on the metric’s scale and variability, but ± 0.10 – 0.15 is a practical guide.

Cross-validation tells us which model performs best on average, but it doesn’t guarantee that the model is perfectly balanced. In this case, our best K (51) generalizes well but is slightly too smooth — it underfits the data.

# Part 6: Handling Imbalanced Data

Now we tackle the elephant in the room: **how to handle class imbalance**!

**Why is imbalance a problem?**

1.  Models learn to predict the majority class (easy way to minimize error)
2.  Minority class often more important (e.g., detecting disease)
3.  Standard evaluation metrics misleading
4.  Decision boundary biased toward majority class

------------------------------------------------------------------------

## 6.1 Solution 1: Class Weights

**Idea**: Penalize misclassifications of the minority class more heavily during training.

```{r class_weights}
# --- Compute class weights ---
class_weights <- heart_train %>%
  count(died) %>%
  mutate(weight = max(n) / n) %>%
  pull(weight, name = died)

heart_train_wt <- heart_train %>%
  mutate(weight = importance_weights(class_weights[died]))

# --- Define weighted recipe ---
weighted_recipe <- recipe(died ~ ., data = heart_train_wt) %>%
  step_dummy(all_nominal_predictors()) %>%       # <-- convert categorical vars to numeric
  step_normalize(all_numeric_predictors())

# --- Weighted logistic regression spec ---
weighted_spec <- logistic_reg(penalty = 0, mixture = 0) %>%  # unregularized logistic regression
  set_engine("glmnet") %>%
  set_mode("classification")

# --- Workflow ---
weighted_workflow <- workflow() %>%
  add_recipe(weighted_recipe) %>%
  add_model(weighted_spec)

# --- Fit model ---
weighted_fit <- fit(weighted_workflow, data = heart_train_wt)

cat("\n✅ Weighted logistic regression fit complete!\n")

```

:bulb: **Weighted Logistic Regression for Imbalanced Data**

This code builds a logistic regression model that gives *more importance* to the minority class ("Died").\
We first calculate **class weights** inversely proportional to class frequency—so rarer outcomes receive higher weights.\
These weights are stored in a special column (`importance_weights()`), which tells tidymodels to use them when fitting the model.

The recipe then: 1. Converts categorical predictors into numeric form with `step_dummy()`,\
2. Standardizes numeric features with `step_normalize()` (important for glmnet).

Finally, the weighted logistic regression is trained using the `glmnet` engine,\
which supports observation-level weights in the model’s loss function.

:point_right: The result: the model penalizes misclassifications of the minority class more heavily,\
helping correct bias from class imbalance *without resampling the data*.

```{r}
class_weights
```

**How it works**:

-   Misclassifying a "Died" case costs \~2x more than misclassifying a "Survived" case
-   Model learns to pay more attention to minority class

------------------------------------------------------------------------

## 6.2 Solution 2: Undersampling the Majority Class

**Idea**: Randomly remove instances from the majority class to balance the dataset.

```{r undersampling}
# Create recipe with undersampling
undersample_recipe <- recipe(died ~ ., data = heart_train) %>%
  step_normalize(all_numeric_predictors()) %>%
  themis::step_downsample(died, under_ratio = 1)  # 1:1 ratio

# Workflow
undersample_workflow <- workflow() %>%
  add_recipe(undersample_recipe) %>%
  add_model(logistic_spec)

# Fit
undersample_fit <- fit(undersample_workflow, data = heart_train)

# Check what happened to the training data
undersampled_data <- prep(undersample_recipe) %>% bake(new_data = NULL)
cat("Original training set:\n")
table(heart_train$died)

cat("\nAfter undersampling:\n")
table(undersampled_data$died)
```

**Pros**:

-   Balanced classes
-   Faster training (less data)

**Cons**:

-   Loss of information (discards majority class data)
-   May reduce model performance

------------------------------------------------------------------------

## 6.3 Solution 3: Oversampling the Minority Class (SMOTE)

**SMOTE (Synthetic Minority Over-sampling Technique)**: Creates synthetic examples of the minority class by interpolating between existing examples.

```{r smote}
# Create recipe with SMOTE
smote_recipe <- recipe(died ~ ., data = heart_train) %>%
  step_dummy(all_nominal_predictors()) %>%   # convert categorical → numeric
  themis::step_smote(died, over_ratio = 1) %>%  # perform SMOTE on numeric data
  step_normalize(all_numeric_predictors())      # normalize after resampling

# Workflow
smote_workflow <- workflow() %>%
  add_recipe(smote_recipe) %>%
  add_model(logistic_spec)

# Fit
smote_fit <- fit(smote_workflow, data = heart_train)

# Check what happened
smote_data <- prep(smote_recipe) %>% bake(new_data = NULL)
cat("After SMOTE:\n")
table(smote_data$died)
```

**How SMOTE works**:

1.  Select a minority class instance
2.  Find its k nearest neighbors (also minority class)
3.  Create synthetic instance by interpolating between them

**Pros**:

-   No information loss (keeps all original data)
-   Creates realistic synthetic examples

**Cons**:

-   Can create noisy examples if classes overlap
-   Increases training time
-   If the dataset contains categorical features, they need to be transformed into dummy variables as SMOTE only works on numeric data

:bulb: **SMOTE-NC: Handling mixed data**

The standard SMOTE algorithm only works on numeric data because it interpolates between samples.

SMOTE-NC (Nominal + Continuous) extends this to mixed data types:

-   Numeric features are interpolated,
-   Categorical features are copied from the nearest neighbor.

Here's how it applies to logistic regression:

```{r}

smotenc_recipe <- recipe(died ~ ., data = heart_train) %>%
  themis::step_smotenc(died, over_ratio = 1) %>%   # 1) SMOTE-NC on mixed data
  step_dummy(all_nominal_predictors()) %>%         # 2) one-hot for glmnet
  step_normalize(all_numeric_predictors())         # 3) standardize for glmnet

log_spec <- logistic_reg(penalty = 0.0, mixture = 0) %>%
  set_engine("glmnet") %>%
  set_mode("classification")

wf <- workflow() %>%
  add_recipe(smotenc_recipe) %>%
  add_model(log_spec)

fit_smotenc <- fit(wf, data = heart_train)

# Check what happened
smotenc_data <- prep(smotenc_recipe) %>% bake(new_data = NULL)
cat("After SMOTENC:\n")
table(smote_data$died)
```

:bulb: Note that we still used a `step_dummy` after `step_smotenc`: that's because logistic regression with `glmnet` can't be applied directly on categorical features and they need to be transformed to dummy variables. It’s not all logistic regression models, only certain engines (implementations) — like `glmnet` — that require fully numeric input. For models that handle categorical features (e.g tree models), the `step_dummy` step would be unnecessary.

------------------------------------------------------------------------

## 6.4 Comparing All Approaches

Let's compare all our approaches on the test set!

```{r compare_approaches}
# Helper to get predictions and truth
get_predictions <- function(fit, test_data, model_name) {
  predict(fit, test_data, type = "prob") %>%
    bind_cols(predict(fit, test_data, type = "class")) %>%
    bind_cols(test_data %>% select(died)) %>%
    mutate(model = model_name)
}

# Get test predictions
pred_undersample <- get_predictions(undersample_fit, heart_test, "Undersampling")
pred_smote       <- get_predictions(smote_fit, heart_test, "SMOTE")
pred_weighted    <- get_predictions(weighted_fit, heart_test, "Class Weights")

# Combine into one dataset
all_preds <- bind_rows(pred_undersample, pred_smote, pred_weighted)

# Detect predicted probability columns automatically
prob_cols <- grep("^\\.pred_", names(all_preds), value = TRUE)
print(prob_cols)
# e.g. c(".pred_Survived", ".pred_Died")

# Use the positive class (usually the second level)
pos_class <- prob_cols[2]

# Compute metrics safely using dynamic column selection
all_metrics <- all_preds %>%
  group_by(model) %>%
  summarise(
    bal_accuracy = bal_accuracy_vec(truth = died, estimate = .pred_class),
    precision    = precision_vec(truth = died, estimate = .pred_class, event_level = "second"),
    recall       = recall_vec(truth = died, estimate = .pred_class, event_level = "second"),
    f_meas       = f_meas_vec(truth = died, estimate = .pred_class, event_level = "second"),
    roc_auc      = roc_auc_vec(truth = died, estimate = .data[[pos_class]], event_level = "second"),
    pr_auc       = pr_auc_vec(truth = died, estimate = .data[[pos_class]], event_level = "second")
  ) %>%
  arrange(desc(f_meas))

all_metrics
```

Let's visualise the results:

```{r}
all_metrics %>%
  pivot_longer(-model, names_to = "metric", values_to = "value") %>%
  group_by(metric) %>%
  mutate(mean_val = mean(value)) %>%  # for ordering metrics top-to-bottom
  ggplot(aes(x = reorder(metric, mean_val), y = value, fill = model)) +
  geom_col(position = position_dodge(width = 0.7), width = 0.7) +
  geom_text(
    aes(label = sprintf("%.2f", value)),
    position = position_dodge(width = 0.7),
    hjust = -0.1,
    size = 3.5,
    color = "black"
  ) +
  coord_flip() +
  scale_y_continuous(limits = c(0, 1.05), expand = c(0, 0)) +
  scale_fill_manual(values = c("Class Weights" = "#E74C3C",
                               "SMOTE" = "#27AE60",
                               "Undersampling" = "#3498DB")) +
  labs(
    title = "Comparison of Imbalance Handling Strategies",
    subtitle = "Trade-offs between precision and recall across methods",
    x = NULL,
    y = "Metric Value",
    fill = "Model"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "top",
    legend.title = element_text(face = "bold"),
    plot.title = element_text(face = "bold")
  ) %>%
  print()
```

> :bulb: **Model comparison under class imbalance** Each strategy shows a distinct trade-off between *precision* and *recall*:
>
> -   **Undersampling** slightly increases recall (more sensitivity to the minority class) but may lose some precision because it discards majority samples.
> -   **SMOTE** achieves a good balance by enriching minority examples synthetically.
> -   **Class weighting** maintains all data but can be less stable when the imbalance is large.
>
> Overall, the **F1-score and PR-AUC** often give the most realistic picture of performance on imbalanced data, while **ROC-AUC** can look deceptively high.
