---
title: "🗓️ Week 03 - From Regression to Classification"
subtitle: "DS202: Data Science for Social Scientists"
author: "Theme: Supervised Learning"
date: "16 October 2025"
format: 
  html:
    page-layout: full
    toc: true
    toc-depth: 2
self-contained: true
---

# Welcome to Week 03! 🎉

**Last week**, we started working through a case study on predicting national happiness. We covered: - Data loading and exploration - Understanding missingness mechanisms (MCAR, MAR, MNAR)

**We stopped at the imputation step.** Today we'll: 1. Complete the regression workflow (imputation → modeling → evaluation) 2. Transition to a new problem: **classification** 3. Learn how to evaluate classification models

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

## ⚙️ Setup

```{r setup, message=FALSE, warning=FALSE}
# Core packages
library(tidyverse)
library(tidymodels)

# Data handling
library(janitor)
library(naniar)
library(visdat)
library(corrplot)

# Feature selection
library(FSelectorRcpp)
library(caret)

# Visualization
library(patchwork)

# Country name normalization
library(countrycode)

#EDA
library(skimr)

theme_set(theme_minimal(base_size = 12))
```

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

# 🎯 Learning Objectives

**Hour 1: Completing the Regression Workflow**

-   Handle missing data with imputation
-   Understand the recipe/prep/bake workflow
-   Compare multiple regression models
-   Apply regularization (Ridge, Lasso, Elastic Net)

**Hour 2: Introduction to Classification**

-   Distinguish regression from classification
-   Build logistic regression models
-   Evaluate with confusion matrices and metrics
-   Perform feature selection

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

# Part 1: Completing the Regression Case Study 📊

## 1.1 Picking Up Where We Left Off

Last week we loaded the happiness data and discussed missingness mechanisms. Let's load it again:

```{r load_data, message=FALSE}
# Load data
wb_data <- read_csv("../data/wb_data_clean_2020.csv")  #adapt filepath to your actual setup!
whr_data <- read_csv("../data/whr_2020.csv") %>% #adapt filepath to your actual setup!
  clean_names() %>%
  rename(happiness_score = ladder_score)

# PROBLEM: Country names differ between World Bank and World Happiness Report datasets
# SOLUTION: Standardize country names to match World Bank conventions before merging
# This prevents losing countries during the join operation

whr_final <- whr_data %>%
  mutate(
    # STEP 0: Clean up any extra spaces at start/end of country names
    country_name = trimws(country_name),
    
    # STEP 1: Manual mapping for known mismatches
    # We explicitly map WHR country names → World Bank country names
    # These are cases where countries have different official names in the two datasets
    country_clean = case_when(
      # European countries with naming differences
      country_name == "Czech Republic" ~ "Czechia",                    # WB uses new official name
      country_name == "Slovakia" ~ "Slovak Republic",                  # WB uses full name
      
      # Asian countries
      country_name == "South Korea" ~ "Korea, Rep.",                   # WB uses formal abbreviation
      country_name == "Russia" ~ "Russian Federation",                 # WB uses full formal name
      country_name == "Kyrgyzstan" ~ "Kyrgyz Republic",               # WB uses alternative spelling
      country_name == "Hong Kong S.A.R. of China" ~ "Hong Kong SAR, China",  # Different SAR format
      country_name == "Vietnam" ~ "Viet Nam",                         # WB uses two-word spelling
      country_name == "Laos" ~ "Lao PDR",                             # WB uses official abbreviation
      
      # Middle Eastern countries
      country_name == "Iran" ~ "Iran, Islamic Rep.",                  # WB includes "Islamic Rep."
      country_name == "Palestinian Territories" ~ "West Bank and Gaza", # WB uses territorial description
      country_name == "Yemen" ~ "Yemen, Rep.",                        # WB includes "Rep."
      country_name == "Turkey" ~ "Turkiye",                           # Updated to new official name (2022)
      
      # African countries
      country_name == "Ivory Coast" ~ "Cote d'Ivoire",               # WB uses French name
      country_name == "Congo (Brazzaville)" ~ "Congo, Rep.",         # WB distinguishes by "Rep."
      country_name == "Congo (Kinshasa)" ~ "Congo, Dem. Rep.",       # WB distinguishes by "Dem. Rep."
      country_name == "Gambia" ~ "Gambia, The",                      # WB includes "The"
      country_name == "Swaziland" ~ "Eswatini",                      # Updated to new official name (2018)
      country_name == "Egypt" ~ "Egypt, Arab Rep.",                   # WB includes "Arab Rep."
      
      # Countries with multiple name variations in the data
      country_name %in% c("Bosnia & Herzegovina", "Bosnia and Herzegovina") ~ "Bosnia and Herzegovina",  # Standardize ampersand
      country_name %in% c("Trinidad & Tobago", "Trinidad and Tobago") ~ "Trinidad and Tobago",          # Standardize ampersand
      country_name %in% c("Myanmar (Burma)", "Burma") ~ "Myanmar",    # Handle historical name
      
      # European countries
      country_name == "Macedonia" ~ "North Macedonia",                # Updated to new official name (2019)
      
      # South American countries
      country_name == "Venezuela" ~ "Venezuela, RB",                  # WB includes "RB" (Bolivarian Republic)
      
      # Disputed territories / special cases
      country_name == "Kosovo" ~ "Kosovo",                            # Keep as is (not in WB but handle explicitly)
      
      # DEFAULT: If no manual mapping found, set to NA (will be handled in Step 2)
      TRUE ~ NA_character_
    ),
    
    # STEP 2: For countries NOT manually mapped (NA values), try automatic standardization
    # Use the countrycode package to standardize country names
    # This catches countries that don't need special handling
    country_clean = if_else(
      is.na(country_clean),                                           # Only apply to unmapped countries
      countrycode(country_name, "country.name", "country.name", warn = FALSE),  # Auto-standardize
      country_clean                                                   # Keep manual mapping if it exists
    )
  )

# RESULT: country_clean column now contains standardized names that will match World Bank data
# Any countries that couldn't be matched will have NA in country_clean

# Merge World Bank and standardised World Happiness datasets
survey_results <- wb_data %>%
  inner_join(whr_final %>% select(country_name, happiness_score),
             by = c("country" = "country_name"))
```

::: callout-notes
I've updated the World Bank a bit and have replaced the urban population variable by the percentage of urban population (wrt the total population).

It makes more sense to use percentages when comparing countries rather than sticking to absolute values (using absolute values would let countries with large populations dominate)!
:::

Let's count missing values per variable

```{r}
# Count total missing values
sum(is.na(survey_results))

# Count missing values per column
colSums(is.na(survey_results))

# Generate a tidy summary
survey_results %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  pivot_longer(
    everything(),
    names_to = "variable",
    values_to = "missing_count"
  ) %>%
  mutate(
    total_rows = nrow(survey_results),
    missing_pct = round((missing_count / total_rows) * 100, 1),
    missing_pct = paste0(missing_pct, "%")
  ) %>%
  arrange(desc(missing_count)) %>%
  knitr::kable()
```

Let's visualise missingness:

```{r}
# Optional: sample rows if too many
max_rows <- 200

df_long <- survey_results %>%
  mutate(.row = row_number()) %>%
  mutate(across(-.row, ~ is.na(.))) %>%  # convert to logical
  { if(nrow(.) > max_rows) filter(., .row %in% sample(unique(.$.row), max_rows)) else . } %>%
  pivot_longer(-.row, names_to = "variable", values_to = "missing") %>%
  group_by(variable) %>%
  mutate(missing_prop = mean(missing)) %>%
  ungroup() %>%
  mutate(variable = fct_reorder(variable, missing_prop, .desc = TRUE))

# 1️⃣ Matrix plot
matrix_plot <- ggplot(df_long, aes(x = .row, y = variable, fill = missing)) +
  geom_tile(width = 1, height = 0.9) +
  scale_fill_manual(values = c("TRUE" = "#007BA7", "FALSE" = "#FEF5D0"), name = "Missing") +
  labs(title = "Missingness Matrix") +
  theme_minimal(base_size = 10) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    panel.grid = element_blank()
  )

# 2️⃣ Side bar plot with percentages
bar_plot <- df_long %>%
  distinct(variable, missing_prop) %>%
  ggplot(aes(x = missing_prop, y = variable)) +
  geom_col(fill = "#007BA7") +
  geom_text(aes(label = paste0(round(missing_prop*100,1), "%")),
            hjust = -0.05, size = 3) +  # annotate outside
  scale_x_continuous(labels = percent_format(accuracy = 1),
                     expand = expansion(mult = c(0, 0.2))) +  # space for text
  labs(x = "% Missing", y = NULL, title = "Missing % per Variable") +
  theme_minimal(base_size = 10) +
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    panel.grid = element_blank()
  )

# 3️⃣ Combine side by side
combined_plot <- matrix_plot + bar_plot + plot_layout(widths = c(5, 1.5))

# Display
combined_plot

# Optional: save wide figure
# ggsave("missingness_matrix_annotated.png", combined_plot, width = 18, height = 12)
```

::: callout-tip
If your variable names are not too long, you could use the `vis_miss` function from `naniar` to produce the missingness matrix:

``` r
vis_miss(survey_results)
```

This doesn't quite work here (i.e the name of variables in the plot appear too high above the missingness matrix and are truncated) since the variable names are quite long (e.g `community_health_workers_per_1000`)
:::

We'll drop the variables that have too many missing values i.e:

-   `community_health_workers_per_1000`
-   `social_safety_net`
-   `adult_literacy_pct`
-   `cpi_gender_eq` (this one is borderline but it is highly likely that imputing will add more noise than actual signal)

```{r}
# Select key variables (remove those with too much missingness)
survey_final <- survey_results |> select(-community_health_workers_per_1000,-social_safety_net,-adult_literacy_pct,-cpi_gender_eq)
```

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

## 1.2 Exploratory Data Analysis (EDA)

Before jumping into modeling, let's understand our data better:

```{r eda_summary}
# Summary statistics
skim(survey_final)
```

The dataset `survey_final` contains 131 observations, each representing a different country, and 26 variables. Of these, four are character variables: `country` and `iso3c` identify the country by name and 3-letter code, while `region` indicates the geographic region and `income` specifies the country's income group. The remaining 22 variables are numeric and measure a range of socio-economic, environmental, and health indicators. Most numeric variables are nearly complete, although some have substantial missingness, notably `gini_index`, `education_upper_secondary_pct`, and `nurses_midwives_per_1000`.

Among the numeric variables, several exhibit notable patterns. Variables such as `rule_of_law`, `happiness_score`, and `life_expectancy` are roughly normally distributed, while others show strong skewness or extreme values. For instance, `population_density` and `gni_ppp` are highly right-skewed, reflecting very dense or very wealthy countries relative to the majority. `ghg_per_capita` is concentrated near zero, but a few countries have substantially higher emissions. Some health-related variables, such as `maternal_death_risk_pct` and `suicide_rate`, have low median values across countries, but certain countries experience very high rates, creating notable outliers. Percent-based variables like `electricity_access_pct` and `internet_users_pct` also vary widely, highlighting disparities in infrastructure and technology access across countries.

**Particularly skewed variables** that may require attention in regression analyses include:

-   `population_density`
-   `gni_ppp`
-   `ghg_per_capita`
-   `maternal_death_risk_pct`

When performing regression, we might consider transformations to reduce skewness. For strictly positive variables, a log-transformation can be applied. For variables that include zeros or negative values, the Yeo-Johnson transformation is more appropriate. In the first instance, however, we will use the raw variables directly while being aware of their distributions.

Overall, the dataset is well-suited for cross-country comparisons and exploratory analyses. The variation and extreme values in economic, demographic, and health indicators make it a rich resource for examining global inequalities and country-level patterns.

### Distribution of Happiness Scores

```{r happiness_distribution, fig.width=10, fig.height=5}
p1 <- ggplot(survey_final, aes(x = happiness_score)) +
  geom_histogram(bins = 30, fill = "steelblue", alpha = 0.7) +
  geom_vline(xintercept = median(survey_final$happiness_score, na.rm = TRUE),
             linetype = "dashed", color = "red", size = 1) +
  labs(title = "Distribution of Happiness Scores",
       subtitle = "Red line shows median - we'll use this to create binary outcome later",
       x = "Happiness Score", y = "Count")

p2 <- ggplot(survey_final, aes(x = region, y = happiness_score, fill = region)) +
  geom_boxplot(show.legend = FALSE) +
  coord_flip() +
  labs(title = "Happiness by Region",
       x = "", y = "Happiness Score")

p1 / p2
```

The plots provide a detailed picture of how happiness scores are distributed across all countries in your dataset and how these scores vary by global region.

**Overall Distribution of Happiness Scores**

-   The top histogram displays a moderately spread distribution of happiness scores, with most values falling between roughly 4 and 7.
-   The distribution is somewhat bimodal, showing clusters of countries with both lower and higher happiness scores.
-   The red dashed line marks the median happiness score (just above 6), which serves as a natural threshold for dividing countries into "higher" and "lower" happiness groups for potential future binary classification.

**Regional Differences in Happiness**

-   The boxplots below the histogram reveal pronounced regional variation:
    -   North America and East Asia & Pacific have notably higher median happiness scores with low variability—most countries in these regions experience consistently high life satisfaction.
    -   Europe & Central Asia and Latin America & Caribbean show a broader spread of scores but still tend toward relatively high medians.
    -   Sub-Saharan Africa, South Asia, and Middle East & North Africa are distinguished by much lower median happiness scores and a wider spread, reflecting more mixed and generally lower wellbeing within these regions.

**Hypothesis for Regional Patterns**

-   These observed differences likely reflect well-established global patterns: regions with higher economic development, stronger institutions, improved health and social services, and stable governance see higher and more homogeneous happiness scores.\[1\]\[2\]\[3\]\[4\]
-   In contrast, regions facing greater economic hardship, political instability, health crises, or inequality exhibit both lower average happiness and greater disparity in life satisfaction among their populations.
-   The broader spread in some regions (such as Latin America & Caribbean) may also signal the impact of cultural or social strengths (like family and community ties) that partially offset economic disadvantages.

**Conclusion**

Together, the plots suggest that happiness is not evenly distributed globally but instead follows clear regional patterns driven by a mix of economic, social, and institutional factors. The median value serves as an effective reference point for analysis, and the boxplots visually reinforce that geography—and the factors that come with it—play a major role in shaping how people report their life satisfaction. Outliers within regions may represent effects of unique national policies, cultural norms, or recent events. This rich, multi-layered view of the data offers a strong starting point for further analysis, predictive modeling, or research into the roots of global happiness disparities.

For more on the analysis of global happiness patterns, see:

-   [1](https://www.worldhappiness.report/ed/2025/caring-and-sharing-global-analysis-of-happiness-and-kindness/)
-   [2](https://www.visualcapitalist.com/mapped-countries-happiest-least-happy/)
-   [3](https://www.happierlivesinstitute.org/2025/03/27/6-lessons-about-world-happiness-from-the-world-happiness-report/)
-   [4](https://www.voronoiapp.com/maps/Mapping-Global-Happiness-By-Country-in-2025-4434)

### Relationships Between Key Variables

```{r correlations, fig.width=10, fig.height=10}
# Select key numeric variables for correlation plot
key_vars <- survey_final %>%
  select(happiness_score, life_expectancy, gni_ppp, rule_of_law, 
         unemployment_pct, education_upper_secondary_pct) %>%
  drop_na()

# Correlation matrix
cor_matrix <- cor(key_vars)

# Visualize
corrplot(cor_matrix, method = "color", type = "upper", 
         addCoef.col = "black", number.cex = 0.7,
         tl.col = "black", tl.srt = 45,
         title = "Correlation Matrix: Happiness and Key Predictors",
         mar = c(0,0,2,0))
```

**What we see:**

-   Strong positive correlation: `life_expectancy` ↔ `happiness_score` (r ≈ 0.75)
-   Strong positive correlation: `gni_ppp` ↔ `happiness_score` (r ≈ 0.77)
-   Multicollinearity alert: `rule_of_law` ↔ `gni_ppp` (r ≈ 0.88) - these are highly correlated with each other!

This suggests regularization (Ridge/Lasso/Elastic Net) will be important!

### Exploring Missingness Patterns

```{r missingness_exploration}
# Which variables have most missing data?
miss_var_summary(survey_final) %>%
  filter(n_miss > 0) %>%
  ggplot(aes(x = reorder(variable, n_miss), y = n_miss)) +
  geom_col(fill = "coral") +
  coord_flip() +
  labs(title = "Missing Data by Variable",
       x = "Variable", y = "Number of Missing Values")
```

**Recap: Missingness Mechanisms**

Last week we discussed three types:

1.  **MCAR (Missing Completely at Random)**: Missingness unrelated to anything
    -   Example: Random survey technical failures
    -   ✅ Safe to ignore or impute
2.  **MAR (Missing at Random)**: Missingness depends on *observed* variables
    -   Example: Low-income countries less likely to report healthcare data
    -   ✅ Can impute using other variables
3.  **MNAR (Missing Not at Random)**: Missingness depends on the *unobserved* value
    -   Example: Countries with high unemployment deliberately don't report it
    -   ⚠️ Most problematic - can bias results

**Today:** We'll handle this missing data properly using imputation!

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

## 1.3 Train/Test Split FIRST (Preventing Data Leakage)

> 🚨 **CRITICAL:** Always split BEFORE any preprocessing!

```{r train_test_split}
set.seed(123)
data_split <- initial_split(survey_final, prop = 0.7, strata = region)

train_data <- training(data_split)
test_data <- testing(data_split)

cat("Training set:", nrow(train_data), "countries\n")
cat("Test set:", nrow(test_data), "countries\n")
```

**Why stratify by region?**

-   Ensures both sets have proportional representation from each region
-   Prevents one set from being "easier" than the other

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

## 1.4 The Recipe/Prep/Bake Workflow

Let's now pre-process our data (e.g imputing missing values).

This is where tidymodels shines! Let's create a **preprocessing recipe**:

```{r create_recipe}
# Define preprocessing steps
survey_recipe <- recipe(happiness_score ~ ., data = train_data) %>%
  
  # 1️⃣ Identify ID variables (don't use as predictors)
  update_role(country, iso3c, region, year, new_role = "id") %>%
  
  # 2️⃣ Handle categorical variables
  step_novel(all_nominal_predictors(), new_level = "novel") %>%  # Handle new levels in test
  step_unknown(all_nominal_predictors(), new_level = "unknown") %>%  # Replace NA with "unknown"
  step_other(all_nominal_predictors(), threshold = 0.01) %>%  # Collapse rare levels
  step_dummy(all_nominal_predictors()) %>%  # Convert to dummy variables
  
  # 3️⃣ Handle missing numeric values
  step_impute_bag(all_numeric_predictors()) %>%  # Impute using bagged trees
  
  # 4️⃣ remove zero-variance variables
  
  step_zv(all_numeric_predictors()) %>%
  
  # 5️⃣ Normalize numeric predictors (CRITICAL for regularization!)
  step_normalize(all_numeric_predictors()) 

# View the recipe
survey_recipe
```

**Code intuition: the cake analogy 🍰**

Think of preprocessing a dataset like **baking a cake**:

-   Each column in the dataset is an ingredient:
    -   **IDs** = ingredients you don’t actually eat 😅
    -   **Categorical variables** = flavors
    -   **Numeric variables** = quantities
-   A **recipe** tells you **how to prep your ingredients** so the cake (model) turns out well.

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

**What's `step_impute_bag()`?**

> It fills missing values by looking at patterns in other variables. Think of it as asking: "What values would make sense here based on similar countries?"
>
> More sophisticated than mean/median imputation because it considers relationships between variables (and less likely to distort the data distribution!).
>
> It also works on datasets that contain a mix of categorical and numeric variables.

**What about `update_role()`?**

> It's mostly used to change the role of variables to a non-predictor role so that they are no longer predictors and are no longer used to build the model.

**What do the steps for handling categorical variables do?**

Sometimes, the test data contains variable categories that did not exist in the training data. This causes errors if not dealt with.

> `step_novel(all_nominal_predictors(), new_level = "novel")` ensures that any new category seen by the model in the test data is labeled "novel"
>
> `step_unknown(all_nominal_predictors(), new_level = "unknown")` replaces missing categories (NAs) with "unknown"
>
> `step_other(all_nominal_predictors(), threshold = 0.01)` collapses rare categories (less than 1% of observations) into "other"
>
> `step_dummy(all_nominal_predictors())` convert categorical variables into dummy/one-hot columns for modeling

**What about `step_normalize()`?**

> It normalizes numeric variables i.e it centers and scales numeric columns to a mean of 0 and a standard deviation of 1. This ensures fair weighting of the variables, which is especially important for regularized models.

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

## 1.5 Understanding Recipe/Prep/Bake

```{r prep_bake}
# 🔧 PREP: Learn parameters from training data ONLY
survey_recipe_prep <- prep(survey_recipe, training = train_data)

# 🍞 BAKE: Apply transformations
train_imputed <- bake(survey_recipe_prep, new_data = NULL)  # NULL = training data
test_imputed <- bake(survey_recipe_prep, new_data = test_data)

glimpse(train_imputed)
```

**The Workflow:**

```         
┌─────────────────────────────────────────────────────────┐
│ recipe()  →  Define WHAT to do                          │
│              (impute, normalize, dummy encode, etc.)    │
└─────────────────────────────────────────────────────────┘
                          ↓
┌─────────────────────────────────────────────────────────┐
│ prep()    →  LEARN parameters from TRAINING data        │
│              (means, modes, imputation models, etc.)    │
└─────────────────────────────────────────────────────────┘
                          ↓
┌─────────────────────────────────────────────────────────┐
│ bake()    →  APPLY transformations to any dataset       │
│              (train, test, or new data)                 │
└─────────────────────────────────────────────────────────┘
```

**Why this prevents leakage:**

-   `prep()` only sees training data
-   `bake()` applies the same transformations to test
-   Test data NEVER influences the preprocessing parameters!

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

## 1.6 Building Regression Models

Now let's fit several models:

```{r fit_models}
# Remove ID columns for modeling
train_modeling <- train_imputed %>% select(-country, -iso3c, -region, -year)
test_modeling <- test_imputed %>% select(-country, -iso3c, -region, -year)

# Model 1: Ordinary Least Squares (OLS)
ols_model <- linear_reg() %>%
  set_engine("lm") %>%
  fit(happiness_score ~ ., data = train_modeling)

# Model 2: Ridge Regression (L2 penalty)
ridge_model <- linear_reg(penalty = 0.1, mixture = 0) %>%
  set_engine("glmnet") %>%
  fit(happiness_score ~ ., data = train_modeling)

# Model 3: Lasso Regression (L1 penalty)
lasso_model <- linear_reg(penalty = 0.1, mixture = 1) %>%
  set_engine("glmnet") %>%
  fit(happiness_score ~ ., data = train_modeling)

# Model 4: Elastic Net (combination)
elastic_model <- linear_reg(penalty = 0.1, mixture = 0.5) %>%
  set_engine("glmnet") %>%
  fit(happiness_score ~ ., data = train_modeling)
```

**Reminder: What's the difference?**

| Model           | Penalty                  | Effect                            |
|------------------|-----------------------|-------------------------------|
| **OLS**         | None                     | Can overfit with many predictors  |
| **Ridge**       | $\lambda \sum \beta_j^2$ | Shrinks all coefficients smoothly |
| **Lasso**       | $\lambda \sum |\beta_j|$ | Shrinks AND selects (some → 0)    |
| **Elastic Net** | Both L1 & L2             | Best of both worlds               |

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

## 1.7 Evaluating Models

```{r evaluate_models}
# Function to evaluate on test set with Adjusted R2
evaluate_regression <- function(model, test_data, model_name) {
  preds <- predict(model, test_data) %>%
    bind_cols(test_data %>% select(happiness_score))
  
  # Basic R2
  r2 <- rsq_vec(preds$happiness_score, preds$.pred)
  
  # Number of observations and predictors
  n <- nrow(test_data)
  p <- ncol(test_data) - 1  # subtract 1 for target variable
  
  # Adjusted R2
  adj_r2 <- 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
  
  tibble(
    Model = model_name,
    RMSE = rmse_vec(preds$happiness_score, preds$.pred),
    MAE = mae_vec(preds$happiness_score, preds$.pred),
    R_squared = r2,
    Adjusted_R2 = adj_r2
  )
}

# Evaluate all models
results <- bind_rows(
  evaluate_regression(ols_model, test_modeling, "OLS"),
  evaluate_regression(ridge_model, test_modeling, "Ridge"),
  evaluate_regression(lasso_model, test_modeling, "Lasso"),
  evaluate_regression(elastic_model, test_modeling, "Elastic Net")
)

results %>%
  arrange(RMSE) %>%
  mutate(across(where(is.numeric), ~round(., 3)))

```

**Interpreting the metrics:**

-   **RMSE**: Average prediction error (same units as happiness score)
    -   Lower is better
    -   Penalizes large errors more
-   **MAE**: Average absolute error
    -   Lower is better
    -   More robust to outliers
-   **R²**: Proportion of variance explained
    -   Higher is better (0 to 1)

**Question:** Which model performs best? Why might regularization help here?

**1. RMSE (Root Mean Squared Error)**

-   Measures typical size of our prediction errors in **the same units as happiness scores**.
-   The happiness scores range from 2.57 to 7.81 → a range of \~5.2.
-   RMSE values: 0.598–0.660 → errors are around **0.6 points on average**, which is roughly **\~11–12% of the total range**.
-   **Interpretation:** Ridge produces the smallest typical error, so its predictions are closer to actual happiness values.

**2. MAE (Mean Absolute Error)**

-   Average absolute difference between predicted and actual scores.
-   Values: 0.441–0.478 → about **0.44–0.48 points**, slightly smaller than RMSE because RMSE penalizes large errors more heavily.
-   Elastic Net has the lowest MAE → fewer large deviations.
-   **Interpretation:** On average, predicted happiness is off by \~0.44–0.48 points, which is reasonably small relative to a 5.2-point scale.

**3. R² (Coefficient of Determination)**

-   Proportion of variance in happiness explained by the model.
-   Values: 0.677–0.750 → **67–75% of variability in happiness scores is explained**.
-   Elastic Net leads slightly (0.750) → explains the most variance.
-   OLS lags behind (0.677) → regularization improves performance, because predictors are correlated and/or numerous.

**4. Adjusted R²**

-   Adjusted R² accounts for **number of predictors** (`p`) relative to observations (`n`) and penalizes unnecessary complexity.
-   Values: 0.241–0.413 → much lower than R²!
-   **Why?** Our models have many predictors, and some contribute little to prediction. Regularization (Ridge, Elastic Net, Lasso) shrinks coefficients, which increases adjusted R² compared to OLS.
-   Elastic Net has the highest adjusted R² (0.413) → best balance of predictive power **and model simplicity**.
-   OLS has a big drop from R² to Adjusted R² (0.677 → 0.241) → adding many variables without shrinkage reduces the model’s effective explanatory power.

**5. Key Takeaways**

1.  **Prediction accuracy (RMSE, MAE):** Ridge is slightly better in RMSE; Elastic Net is better in MAE. Both are better than OLS or Lasso.

2.  **Variance explained (R²):** Elastic Net \> Ridge \> Lasso \> OLS.

3.  **Model parsimony (Adjusted R²):** Elastic Net best balances prediction and complexity; OLS has overfitting concerns.

4.  **Scale context:**

    -   Errors (\~0.44–0.60) are modest relative to the 5.2-point range.
    -   These models give reasonably precise predictions for happiness, but there’s still **\~25–30% of variance unexplained** even with regularization.

**Summary Judgment:**

-   **Best predictive model overall:** Ridge (lowest RMSE) or Elastic Net (best combination of MAE, R², adjusted R²).
-   **OLS:** simplest, but loses explanatory efficiency once many predictors are included.
-   **Lasso:** performs slightly worse here, possibly because fully zeroing coefficients removed useful variables.

To confirm the presence/absence of overfitting, let's calculate the metrics on the training data too!

```{r evaluate_models_training}
# Function to evaluate on test set with Adjusted R2
evaluate_regression <- function(model, train_data, model_name) {
  preds <- predict(model, train_data) %>%
    bind_cols(train_data %>% select(happiness_score))
  
  # Basic R2
  r2 <- rsq_vec(preds$happiness_score, preds$.pred)
  
  # Number of observations and predictors
  n <- nrow(train_data)
  p <- ncol(train_data) - 1  # subtract 1 for target variable
  
  # Adjusted R2
  adj_r2 <- 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
  
  tibble(
    Model = model_name,
    RMSE = rmse_vec(preds$happiness_score, preds$.pred),
    MAE = mae_vec(preds$happiness_score, preds$.pred),
    R_squared = r2,
    Adjusted_R2 = adj_r2
  )
}

# Evaluate all models
results <- bind_rows(
  evaluate_regression(ols_model, train_modeling, "OLS"),
  evaluate_regression(ridge_model, train_modeling, "Ridge"),
  evaluate_regression(lasso_model, train_modeling, "Lasso"),
  evaluate_regression(elastic_model, train_modeling, "Elastic Net")
)

results %>%
  arrange(RMSE) %>%
  mutate(across(where(is.numeric), ~round(., 3)))
```

**1. RMSE and MAE**

-   Errors on the training set are **smaller than on the test set**:

    -   RMSE: 0.489–0.584 (train) vs 0.598–0.660 (test)
    -   MAE: 0.399–0.483 (train) vs 0.441–0.478 (test)

-   This is typical — models perform better on the data they were trained on.

-   **Interpretation:** On average, predictions are off by **\~0.4–0.58 points**, which is \~8–11% of the total happiness range.

**2. R² and Adjusted R²**

-   **R² (variance explained)** is higher on training than test: 0.753–0.818 (train) vs 0.677–0.750 (test).
-   **Adjusted R²** accounts for number of predictors: 0.667–0.755 (train) vs 0.241–0.413 (test).
-   The **drop from training to test** is especially dramatic for Adjusted R² on OLS (0.755 → 0.241).

**3. What this tells us**

-   **OLS:** Very high training R² and Adjusted R² → fits training data extremely well.

-   **Drop on test set:** Indicates **overfitting** — OLS captures noise in training set.

-   **Ridge & Elastic Net:**

    -   Slightly higher RMSE/MAE on train vs OLS, but much **smaller drop on test**.
    -   Regularization reduces overfitting and improves generalization.

-   **Lasso:** Higher training error → has aggressively zeroed coefficients, loses some predictive power.

**4. Key Takeaways**

1.  **Overfitting detection:** OLS performs best on training but poorly on test → classic overfit.
2.  **Generalization:** Ridge or Elastic Net perform slightly worse on train but **much better on test** → these are our best choices.
3.  **Adjusted R² is very telling:** Large drop for OLS highlights that adding predictors without regularization inflates R² artificially.
4.  **Scale context:** RMSE/MAE errors are relatively small (\~0.4–0.6 points) compared to total happiness score range (\~5.2 points).

**Summary Judgment:**

-   **Predictive focus:** Ridge or Elastic Net are best for test set predictions.
-   **Model simplicity vs. performance:** Elastic Net balances variance explained and prediction accuracy, especially if we want to shrink unimportant coefficients.
-   **OLS:** Good for insight if we want full variable coefficients, but beware overfitting.

Let's visualise our train/test metrics:

```{r}
# Combine train and test results dynamically
train_results <- bind_rows(
  evaluate_regression(ols_model, train_modeling, "OLS"),
  evaluate_regression(ridge_model, train_modeling, "Ridge"),
  evaluate_regression(lasso_model, train_modeling, "Lasso"),
  evaluate_regression(elastic_model, train_modeling, "Elastic Net")
) %>% mutate(Dataset = 'Train')


test_results <- bind_rows(
  evaluate_regression(ols_model, test_modeling, "OLS"),
  evaluate_regression(ridge_model, test_modeling, "Ridge"),
  evaluate_regression(lasso_model, test_modeling, "Lasso"),
  evaluate_regression(elastic_model, test_modeling, "Elastic Net")
) %>% mutate(Dataset = 'Test')

# Combine train and test
all_results <- bind_rows(train_results, test_results)

# Reshape for ggplot
all_results_long <- all_results %>%
  pivot_longer(cols = c(RMSE, MAE, R_squared, Adjusted_R2),
               names_to = 'Metric',
               values_to = 'Value')

# Plot
ggplot(all_results_long, aes(x = Model, y = Value, fill = Dataset)) +
  geom_col(position = position_dodge(width = 0.8)) +
  facet_wrap(~Metric, scales = 'free_y') +
  labs(title = 'Train vs Test Performance Across Models',
       y = 'Metric Value', x = 'Model') +
  theme_minimal() +
  theme(text = element_text(size = 12))

```

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

## 1.8 Understanding Regularization Effects

Let's visualize what our regularization methods (i.e Lasso, Ridge and ElasticNet) did.

Let's first visualise this with a dot plot:

```{r}
# Get union of top features from each method
top_features <- all_coefs %>%
  group_by(model) %>%
  arrange(desc(abs(estimate))) %>%
  slice(1:10) %>%
  ungroup() %>%
  pull(term) %>%
  unique()

# Filter and plot
all_coefs %>%
  filter(term %in% top_features) %>%
  ggplot(aes(x = estimate, y = reorder(term, estimate), color = model, shape = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray40", size = 0.8) +
  geom_point(size = 4, alpha = 0.8) +
  geom_line(aes(group = term), color = "gray70", size = 0.3) +
  labs(
    title = "Feature Coefficients Across Regularization Methods",
    subtitle = "Connected dots show how each method treats the same feature",
    x = "Coefficient Value",
    y = NULL,
    color = "Method",
    shape = "Method"
  ) +
  scale_color_manual(values = c("Ridge" = "#E69F00", "Lasso" = "#56B4E9", "Elastic Net" = "#009E73")) +
  scale_shape_manual(values = c("Ridge" = 16, "Lasso" = 17, "Elastic Net" = 15)) +
  theme_minimal(base_size = 12) +
  theme(
    legend.position = "top",
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_line(color = "gray90"),
    axis.text.y = element_text(size = 11)
  )
```

A second possible visualisation is this heatmap:

```{r coefficients_heatmap, fig.width=10, fig.height=6}
# Create coefficient matrix for top features
coef_matrix <- all_coefs %>%
  group_by(model) %>%
  arrange(desc(abs(estimate))) %>%
  slice(1:12) %>%
  ungroup() %>%
  pull(term) %>%
  unique() %>%
  {filter(all_coefs, term %in% .)} %>%
  select(term, model, estimate) %>%
  pivot_wider(names_from = model, values_from = estimate, values_fill = 0)

# Plot heatmap
coef_matrix %>%
  pivot_longer(cols = -term, names_to = "model", values_to = "estimate") %>%
  mutate(model = factor(model, levels = c("Ridge", "Lasso", "Elastic Net"))) %>%
  ggplot(aes(x = model, y = reorder(term, estimate), fill = estimate)) +
  geom_tile(color = "white", size = 0.5) +
  geom_text(aes(label = sprintf("%.2f", estimate)), size = 3) +
  scale_fill_gradient2(
    low = "coral", 
    mid = "white", 
    high = "steelblue",
    midpoint = 0,
    name = "Coefficient"
  ) +
  labs(
    title = "Coefficient Heatmap: Comparing Regularization Methods",
    subtitle = "White = zero coefficient (feature eliminated)",
    x = NULL,
    y = NULL
  ) +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, face = "bold"),
    axis.text.y = element_text(size = 10),
    panel.grid = element_blank(),
    legend.position = "right"
  )
```

#### Understanding Regularization Effects: A Tale of Three Approaches

Looking at the coefficient plots, we see a fascinating story unfold - three different regularization methods looking at the same data and reaching different conclusions about what matters.

**The Dominant Signal: GNI Per Capita**

All three methods agree on one thing: **wealth (`gni_ppp`) is the strongest predictor of national happiness.** Notice how all three dots cluster far to the right - Ridge (orange circle), Elastic Net (green square), and Lasso (blue triangle) all give this variable a substantial positive coefficient. When methods disagree on almost everything else but unanimously agree on this, you can be confident: economic prosperity is paramount.

**Where the Methods Diverge**

Now look what happens with everything else - this is where the story gets interesting:

**Infrastructure variables** (`electricity_access_pct`, `drinking_water_pct`, `urban_population_pct`) show the first signs of disagreement. All three methods keep them, but notice the spacing: - Ridge (orange) tends to give slightly larger coefficients - Lasso (blue) is more aggressive, shrinking them more - Elastic Net (green) sits in between, playing diplomat

These aren't just technical differences - they represent fundamentally different philosophies about how infrastructure contributes to happiness. Ridge says: "Each infrastructure element contributes independently." Lasso says: "They're redundant - electricity access captures most of the infrastructure story." Elastic Net says: "They overlap somewhat, but each adds something unique."

**Healthcare workforce** (`nurses_midwives_per_1000`, `physicians_per_1000`) reveals the classic multicollinearity problem. These variables likely correlate with wealth and each other. Look at the spread: - Ridge keeps both with moderate coefficients - Lasso aggressively shrinks `physicians_per_1000` toward zero (small blue triangle near zero) - Elastic Net keeps both but smaller than Ridge

The connected line between these dots stretches considerably - a visual signature of multicollinearity. The methods can't agree because both variables are telling similar stories, and there's no "right" answer to which one matters more.

**The Casualties of Lasso's Aggression**

Now we come to the heart of why Lasso performed worse than the other two methods. Look at these features clustering near zero:

**`gini_index`** (income inequality):

-   Ridge (orange) says: "This matters! Give it a positive coefficient around 0.27"
-   Elastic Net (green) agrees but is more cautious: around 0.05
-   Lasso (blue) nearly eliminates it: basically zero

Wait - *positive* coefficient? Countries with *more* inequality are happier? This seems counterintuitive until you realize we're looking at a model that already accounts for wealth. The Ridge interpretation might be: "Among countries of similar wealth, some inequality can coexist with happiness if other factors (rule of law, social mobility) are strong." Or it could simply be that `gini_index` is so entangled with other variables that its coefficient is unstable. The fact that Elastic Net keeps it (albeit small) while Lasso eliminates it entirely shows this uncertainty.

**`rule_of_law`**:

-   Ridge keeps it with a small positive effect
-   Elastic Net also keeps it (very small)
-   Lasso zeros it out completely

This is a significant loss. Rule of law represents governance quality, corruption levels, and institutional trust - all plausibly important for happiness independent of wealth. Lasso's decision to eliminate it entirely suggests it saw `rule_of_law` as redundant with `gni_ppp` (wealthier countries tend to have better rule of law). But Ridge and Elastic Net preserved it, suggesting there's signal beyond just wealth correlation.

**Environmental variables** (`ghg_per_capita`, `renewable_energy_pct`):

-   Ridge gives them small coefficients
-   Elastic Net shrinks them further
-   Lasso essentially eliminates them

All three methods agree these matter less than economic and infrastructure factors. But Ridge and Elastic Net say: "Not zero, though." This is the difference between saying "environment is less important" versus "environment doesn't matter at all."

**The Unemployment Story**

The only negative coefficient in the entire plot belongs to `unemployment_pct` (far left). Interestingly:

-   All three methods agree unemployment hurts happiness
-   But they disagree on magnitude: Ridge thinks it's more important (coefficient ≈ -0.25) than Lasso does (≈ -0.13)
-   Elastic Net splits the difference (≈ -0.18)

The connected line here is moderately long - suggesting some uncertainty about unemployment's independent effect once you account for wealth, infrastructure, and other factors. But unlike some other variables, all three methods kept it with a negative sign. This is strong evidence that unemployment has a real, happiness-reducing effect beyond just its correlation with poverty.

**The Population Density Puzzle**

Look at `population_density` - all three methods give it a *negative* coefficient, though Ridge's is notably more negative than the others. This suggests that controlling for everything else (wealth, infrastructure, urbanization), higher population density might slightly reduce happiness. Perhaps reflecting overcrowding effects or loss of space? The fact that all three agree on the direction gives this finding more credibility.

**Income Categories: A Study in Disagreement**

The `income_Lower.middle.income` category shows one of the longest connected lines in the plot:

-   Ridge: Strong negative effect (≈ -0.15)
-   Elastic Net: Moderate negative effect (≈ -0.08)
-   Lasso: Similar to Elastic Net (≈ -0.08)

This long line screams multicollinearity. Income categories obviously correlate with `gni_ppp` - they're both measures of wealth! The question is: do income *categories* add information beyond continuous GNI? Ridge says yes (and penalizes lower-middle-income countries heavily). Lasso and Elastic Net are more skeptical, shrinking this effect substantially.

**Why Ridge and Elastic Net Outperformed Lasso**

Our test set results showed Ridge and Elastic Net performed better than Lasso:

-   **Ridge**: RMSE = 0.598, R² = 0.740
-   **Elastic Net**: RMSE = 0.610, R² = 0.750 ✓ Best
-   **Lasso**: RMSE = 0.660, R² = 0.702

The plots explains why:

**Lasso was too aggressive.** Look at all those blue triangles clustered near zero - `gini_index`, `rule_of_law`, `physicians_per_1000`, environmental variables. Lasso decided these features were redundant noise and eliminated them. But our test results suggest they contained useful signal, even if weak or correlated.

Think about what Lasso threw away:

-   Governance quality (`rule_of_law`)
-   Income inequality effects (`gini_index`)
-   Some healthcare nuance (`physicians_per_1000`)
-   Environmental factors

These variables might be correlated with wealth, but they capture *different aspects* of development. A wealthy country with terrible governance isn't the same as a wealthy country with strong institutions. Lasso's assumption that you can represent all these dimensions with `gni_ppp` alone proved too simplistic.

**Ridge kept everything.** Notice how the orange circles span the widest range and none are truly at zero. Ridge said: "I'll shrink everything smoothly, but I won't eliminate anything." For our highly correlated development indicators, this democratic approach worked well.

Ridge's philosophy: "Yes, these variables are correlated, but each one might capture a slightly different facet of what makes countries happy. Better to include them all with smaller coefficients than to arbitrarily pick winners."

This proved wise. The small coefficients Ridge assigned to "minor" features like `rule_of_law` and `renewable_energy_pct` apparently did help prediction, even if they weren't stars of the show.

**Elastic Net found the middle path.** The green squares typically sit between Ridge and Lasso. It eliminated some features (or nearly so), but kept others that Lasso discarded. This balanced approach delivered the best test performance - aggressive enough to simplify, gentle enough to retain useful correlated information.

Elastic Net essentially said: "Lasso, you're right that we should eliminate *some* redundancy. But you went too far. Ridge, you're right that most features contain *some* signal. But your model is too complex." The result was a Goldilocks model - just right.

**The Multicollinearity Lesson**

The connected lines tell the real story. When you see a horizontal line (all three dots at similar values), the methods agree. When you see a long line stretching across the plot, you're looking at a feature where multicollinearity created disagreement:

-   **`gini_index`**: Very long line → High disagreement → Definitely correlated with other features
-   **`gni_ppp`**: Short line → High agreement → Clear, strong, independent signal
-   **`income_Lower.middle.income`**: Long line → Methods can't agree → Redundant with GNI
-   **`unemployment_pct`**: Moderate line → Some disagreement but directional consensus

These line lengths are diagnostic. Short lines = robust findings. Long lines = be cautious in interpretation.

When features are highly correlated (like our development indicators), Lasso's "pick one winner" strategy backfires. Ridge's "keep everyone" and Elastic Net's "keep most" strategies prove superior because they acknowledge that correlation ≠ perfect redundancy.

**What This Means for Understanding Happiness**

Each method tells a slightly different story about national happiness:

**Lasso's narrative** (the minimalist): \> "Happiness is primarily about wealth (`gni_ppp`). Beyond that, it's about basic infrastructure (electricity, water, urbanization), some healthcare access (nurses), and employment. Everything else - governance, inequality, environment, physician density - is either noise or redundant with wealth."

**Ridge's narrative** (the pluralist): \> "Happiness is primarily about wealth, but it's a complex web. Infrastructure matters significantly. Healthcare matters (both nurses AND doctors). Governance quality matters. Even environmental factors and inequality play small roles. Population density creates modest challenges. Lower-income countries face compounded disadvantages beyond just GDP."

**Elastic Net's narrative** (the pragmatist): \> "Happiness is primarily about wealth. Infrastructure and healthcare are clearly important. Governance and inequality matter, but less than infrastructure. Environmental factors are marginal. Employment clearly matters. We should acknowledge these various dimensions without overcomplicating the story."

**Which narrative is "true"?** The test set performance suggests Elastic Net's balanced view is most accurate for *prediction*. But if your goal is *understanding*, Ridge's more inclusive view might be preferable - it acknowledges that multiple factors matter even if some effects are small.

For *policy*, you might actually prefer something between Elastic Net and Ridge:

-   **Focus resources** on the big levers: economic development, infrastructure, healthcare
-   **Don't ignore** governance, employment policy, and social equity just because their effects are smaller
-   **Recognize** that these dimensions interact - good governance might help translate wealth into happiness more effectively

**The Philosophical Dimension**

These plots reveal a deep tension in statistical modeling:

**Occam's Razor** (Lasso's philosophy): Prefer simpler explanations. If wealth explains most of the happiness variance, and other factors are correlated with wealth, why complicate things?

**Pluralism** (Ridge's philosophy): Social phenomena are complex. Just because factors correlate doesn't mean they're the same thing. Wealth, governance, inequality, environment - these are distinct dimensions that all contribute.

**Pragmatism** (Elastic Net's philosophy): Simplicity is valuable, but not at the cost of accuracy. Keep what matters, even if slightly redundant, but don't include everything.

The fact that Elastic Net achieved the best test performance suggests pragmatism wins for prediction. But remember: prediction and explanation aren't the same thing. A model can predict well with a simplified story, while the true causal structure remains more complex.

**The Practical Takeaway**

These plots beautifully illustrate why regularization isn't one-size-fits-all:

-   **Use Ridge** when features are correlated and you suspect most matter (even if weakly). When you want to preserve all dimensions of a complex phenomenon. When interpretation of many features is valuable.

-   **Use Lasso** when you have many truly irrelevant features and need aggressive selection. When you want a sparse, simple model and are willing to sacrifice some accuracy. When you believe most features are just different measurements of a few underlying factors.

-   **Use Elastic Net** when you're uncertain - it hedges your bets and often wins. When you have correlated features but don't want to eliminate entire dimensions. When you want balance between interpretability and accuracy.

For predicting national happiness with its web of correlated development indicators, Elastic Net's diplomatic balance of shrinkage and selection proved optimal. It was smart enough to simplify (eliminating near-redundancy), but wise enough not to over-simplify (preserving distinct dimensions like governance and inequality that Lasso discarded).

**Final Thought: The Danger of Algorithmic Overconfidence**

Lasso made decisive, confident choices: eliminate this, keep that. Its test performance suggests some of those decisions were wrong. Ridge made cautious, inclusive choices: keep everything but shrink it. This proved more robust.

The lesson? When features are correlated (as they often are with real-world data), be wary of methods that make aggressive elimination decisions. Correlation doesn't mean perfect redundancy. Small signals from multiple correlated features might collectively matter more than you think.

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

## 1.9 Residuals Analysis

```{r residuals_plot, fig.width=10, fig.height=5}
# Compare residuals for best models
plot_residuals <- function(model, data, title) {
  preds <- predict(model, data) %>%
    bind_cols(data %>% select(happiness_score))
  
  ggplot(preds, aes(x = happiness_score, y = happiness_score - .pred)) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
    geom_point(alpha = 0.5) +
    geom_smooth(method = "loess", se = FALSE, color = "blue") +
    labs(
      title = title,
      x = "Actual Happiness Score",
      y = "Residuals"
    )
}

p1 <- plot_residuals(ols_model, test_modeling, "OLS")
p2 <- plot_residuals(elastic_model, test_modeling, "Elastic Net")

p1 + p2
```

#### Interpretation of the residual plots

Both residual plots reveal non-random patterns, indicating model misspecification. The OLS model (left) shows a pronounced U-shaped trend—under-predicting at both low and high happiness levels—suggesting it fails to capture non-linear relationships in the data. The Elastic Net model (right) improves overall by flattening this curve, but a **clear upward trend remains at happiness scores above \~6.5**, meaning the model consistently **under-predicts** for the happiest individuals.

This residual pattern is highly suggestive of **non-linear effects driven by skewed predictors**—common in well-being data (e.g., GDP, social support, or life expectancy). When such predictors are right-skewed and included in a linear model without transformation, their diminishing marginal effects on happiness are mis-modeled, leading to systematic errors at the outcome extremes.

Importantly, many predictors may contain zeros or negative values (e.g., standardized scores), making log or Box-Cox transformations inappropriate. A more robust approach is to apply the **Yeo-Johnson transformation**, which accommodates any real-valued data. In R with **tidymodels**, this is cleanly implemented within a preprocessing recipe:

``` r
rec <- recipe(happiness_score ~ ., data = train_data) %>%
  step_YeoJohnson(all_numeric_predictors())
```

Adding this step before fitting Elastic Net can help linearize relationships, reduce skew-induced bias, and potentially eliminate the residual trend at high happiness scores. If the pattern persists even after transformation, it may signal the need for more flexible modeling approaches (e.g., GAMs or tree-based methods) to capture complex non-linearities.

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

# Part 2: Introduction to Classification 🎯

## 2.1 From Regression to Classification: Same Data, Different Question

**Last hour** we predicted: *"What is a country's happiness score?"* (continuous: 3.5, 6.2, 7.8, etc.)

**This hour** we'll predict: *"Is a country happy or unhappy?"* (binary: Yes/No)

**The key insight:** Same dataset, same predictors, but fundamentally different type of problem!

```{r create_classification_outcome}
# Create binary classification outcome
happy_classif <- survey_final %>%
  drop_na(happiness_score) %>%
  mutate(
    # Countries above median = "Happy"
    is_happy = factor(
      if_else(happiness_score >= median(happiness_score), "Yes", "No"),
      levels = c("No", "Yes")  # "No" is reference level
    )
  )

# Check the split
happy_classif %>%
  count(is_happy) %>%
  mutate(pct = round(n / sum(n) * 100, 1))
```

Perfect! We have a balanced dataset (roughly 50/50 split).

```{r visualize_split, fig.width=10, fig.height=5}
# Visualize what we just did
ggplot(happy_classif, aes(x = happiness_score, fill = is_happy)) +
  geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  geom_vline(xintercept = median(happy_classif$happiness_score),
             linetype = "dashed", size = 1, color = "black") +
  scale_fill_manual(values = c("No" = "coral", "Yes" = "lightgreen"),
                    name = "Happy Country?") +
  labs(
    title = "Creating a Binary Classification from Continuous Happiness Score",
    subtitle = "Countries above median (dashed line) = 'Happy'",
    x = "Happiness Score",
    y = "Count"
  )
```

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

## 2.2 Regression vs Classification: What's the Difference?

```{r compare_tasks, echo=FALSE}
tibble(
  Aspect = c("Outcome Type", "Example Output", "Prediction Range", "Algorithm Used"),
  Regression = c(
    "Continuous number",
    "Happiness score = 6.347",
    "Any real number (e.g., 2.5 to 7.8)",
    "Linear regression (OLS, Ridge, Lasso)"
  ),
  Classification = c(
    "Category/Label",
    "Is Happy = Yes",
    "Fixed categories (Yes/No)",
    "Logistic regression"
  )
) %>%
  knitr::kable()
```

**Key difference:**

-   **Regression**: "How much?" (predicting quantity)
-   **Classification**: "Which group?" (predicting category)

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

## 2.3 Why Linear Regression Fails

Let's see what happens if we try to use linear regression for our binary outcome:

```{r bad_linear, fig.width=10, fig.height=5}
# Convert to numeric (DON'T DO THIS!)
happy_bad <- happy_classif %>%
  mutate(is_happy_numeric = if_else(is_happy == "Yes", 1, 0))

# Fit linear regression
bad_model <- lm(is_happy_numeric ~ life_expectancy + gni_ppp, data = happy_bad)

# Plot predictions
ggplot(happy_bad, aes(x = life_expectancy, y = is_happy_numeric, color = is_happy)) +
  geom_point(alpha = 0.5, size = 3) +
  geom_smooth(method = "lm", se = FALSE, color = "black", size = 1.5) +
  geom_hline(yintercept = c(0, 1), linetype = "dashed", color = "red") +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = 1, ymax = Inf,
           fill = "red", alpha = 0.1) +
  annotate("rect", xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = 0,
           fill = "red", alpha = 0.1) +
  annotate("text", x = 70, y = 1.15, 
           label = "Probabilities > 1?\nNonsense!", 
           color = "red", size = 5, fontface = "bold") +
  annotate("text", x = 55, y = -0.15, 
           label = "Negative probabilities?\nAlso nonsense!", 
           color = "red", size = 5, fontface = "bold") +
  labs(
    title = "Why Linear Regression Fails for Classification",
    subtitle = "Predictions escape [0,1] - what's a probability of -0.2 or 1.3?",
    x = "Life Expectancy (years)",
    y = "Predicted 'Probability'",
    color = "Actually Happy?"
  ) +
  scale_color_manual(values = c("No" = "coral", "Yes" = "lightgreen"))
```

**Problems:**

1.  Predictions can be \< 0 or \> 1 (nonsensical probabilities!)
2.  Assumes linear relationship between predictors and probability
3.  Treats categories as numbers (0, 1) which doesn't make mathematical sense

We need a better approach! ➡️

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

## 2.4 Enter Logistic Regression

Logistic regression models the probability that a given input belongs to class 1 by applying the **logistic (sigmoid) function** to a linear combination of features:

$$
P(Y = 1 \mid X) = \frac{1}{1 + e^{-(\beta_0 + \beta_1 X_1 + \ldots + \beta_p X_p)}}
$$

This sigmoid function transforms any real input into a value between 0 and 1, making it ideal for predicting probabilities.

### Key Properties of the Sigmoid Function:

-   Outputs values strictly between 0 and 1
-   Smooth, S-shaped curve
-   $$ P = 0.5 $$ when the input to the sigmoid is 0

```{r}
# Visualize the sigmoid function
x <- seq(-6, 6, 0.1)
sigmoid <- function(x) 1 / (1 + exp(-x))

tibble(x = x, y = sigmoid(x)) %>%
  ggplot(aes(x, y)) +
  geom_line(size = 2, color = "steelblue") +
  geom_hline(yintercept = c(0, 0.5, 1), linetype = "dashed", alpha = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) +
  annotate("text", x = 0, y = 0.5, label = "P = 0.5 when input = 0",
           hjust = -0.1, vjust = -0.5, size = 5) +
  labs(
    title = "The Sigmoid Function",
    subtitle = "Always outputs between 0 and 1 - perfect for probabilities!",
    x = "Linear Combination (β0 + β1X1 + ...)",
    y = "Predicted Probability"
  )
```

### The Shared Principle: Minimization of an Appropriate Loss Function

At the core of logistic regression, as with many other supervised learning models (e.g linear regression), lies the principle of **minimizing a loss function**. The loss function measures how well the model’s predictions align with the actual observed outcomes in the training data. The goal during training is to find the parameters $\beta$ that minimize this loss.

In logistic regression, the appropriate loss function is the **logistic loss** (also called log loss or binary cross-entropy loss), which is defined as:

$$
\text{Log Loss} = -\frac{1}{N} \sum_{i=1}^N \left[ y_i \log(p_i) + (1 - y_i) \log(1 - p_i) \right]
$$

where

-   $N$ is the number of data points,
-   $y_i \in \{0,1\}$ is the true label for data point $i$,
-   $p_i$ is the predicted probability by the model that $y_i = 1$.

The logistic loss penalizes predictions that are confidently wrong more heavily than those that are less confident. For example, predicting a probability close to 0 when the true label is 1 results in a large loss. Minimizing this loss ensures predicted probabilities converge towards true class labels as closely as possible.

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

## 2.5 Preparing Data for Classification

```{r prep_classification}
# Split data (using stratified sampling)
set.seed(456)
happy_split <- initial_split(happy_classif, prop = 0.7, strata = is_happy)

happy_train <- training(happy_split)
happy_test <- testing(happy_split)

# Create preprocessing recipe (reusing what we learned!)
happy_recipe <- recipe(is_happy ~ ., data = happy_train) %>%
  # Exclude ID / metadata variables
  update_role(country, iso3c, region, happiness_score, new_role = "id") %>%
  
  # Handle categorical predictors robustly
  step_novel(all_nominal_predictors()) %>%        # Handle new factor levels at prediction time
  step_unknown(all_nominal_predictors()) %>%      # Recode missing factor levels
  step_other(all_nominal_predictors(), threshold = 0.01) %>%  # Collapse rare levels
  step_dummy(all_nominal_predictors(), one_hot = FALSE) %>%   # Create dummy variables
  
  # Preprocess numeric predictors
  step_impute_bag(all_numeric_predictors()) %>%   # Impute missing values (bagging works well)
  step_nzv(all_numeric_predictors()) %>%          # Remove near-zero variance predictors (more flexible than step_zv)
  step_YeoJohnson(all_numeric_predictors()) %>%   # Transform skewed numeric predictors (handles 0/negatives)
  step_normalize(all_numeric_predictors())        # Center & scale (required for penalized logistic regression)

# Prep and bake
happy_prep <- prep(happy_recipe)
happy_train_processed <- bake(happy_prep, new_data = NULL)
happy_test_processed <- bake(happy_prep, new_data = happy_test)
```

**Note:**

We're reusing almost the same preprocessing workflow we learned for regression! The recipe/prep/bake pattern works for both. The only modifications we are making are actually based on what we learnt from the regression models! For example, removing predictors with almost no variation -`step_nzv()`- (instead of simply zero-variance predictors) as these add noise and can destabilize penalized fits. Or applying the Yeo-Johnson transformation on the numeric variables (`step_YeoJohnson()`) to deal with their skewness.

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

## 2.6 Logistic Regression & Coefficient Interpretation

### 2.6.1 What We Are Doing

So far, we predicted **happiness scores** as continuous numbers using linear regression. Now we want to predict a **category**: is a country “happy” or “unhappy”?

**Key difference:**

| Regression            | Classification               |
|-----------------------|------------------------------|
| Predict a number      | Predict a category           |
| Linear relationship   | Probability of class         |
| Can predict any value | Must predict between 0 and 1 |

We’ll use **logistic regression**.

**In plain words:** Logistic regression finds a **relationship between predictors (X)** and the **probability of being happy (p)**. It ensures the predicted probability is always between 0 and 1.

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

### 2.6.2 Fit the Logistic Regression Model

```{r}
logistic_model <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification") %>%
  fit(is_happy ~ ., 
      data = happy_train_processed %>% 
        select(-country, -iso3c, -region, -happiness_score))
```

-   This fits a **generalized linear model (GLM) with a logit link**, i.e., logistic regression.
-   The model predicts **log-odds of being happy (Yes)**, not probability directly:

$$\text{log-odds}(\textrm{is_happy = Yes}) = \beta_0 + \beta_1 X_1 + \dots + \beta_p X_p$$

-   $(X_j)$ = predictors (life expectancy, GNI, etc.)
-   $(\beta_j)$ = coefficients
-   $(\beta_0)$ = intercept

**Reference level note:**

-   The first factor level (`No`) is the **reference category**.
-   The model predicts the probability of the **second level** (`Yes`).
-   Coefficients describe the effect of predictors on the **log-odds of being “Yes” vs “No”**.
-   Changing the reference level changes the baseline; coefficients adjust to reflect the new reference, but interpretation is still relative to that baseline.

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

### 2.6.3 Interpreting Coefficients

#### Step 1: Standardization

All numeric predictors were **centered and scaled**:

-   **Centering**: subtract the mean
-   **Scaling**: divide by standard deviation (SD)

> One unit = **1 SD above or below average**.

**Example:**

-   Life expectancy: mean = 72 years, SD = 7 years → 1 SD above mean = 79 years
-   GNI per capita: mean = 15,000 USD, SD = 10,000 USD → 1 SD above mean = 25,000 USD

Why standardize?

-   Makes coefficients **comparable across features**
-   Keeps predicted log-odds on a reasonable scale

#### Step 2: Coefficients

```{r}
tidy(logistic_model) %>%
  mutate(
    odds_ratio = exp(estimate),
    interpretation = case_when(
      term == "(Intercept)" ~ "Baseline",
      odds_ratio > 2 ~ "Strong positive effect",
      odds_ratio > 1.2 ~ "Moderate positive effect",
      odds_ratio < 0.5 ~ "Strong negative effect",
      odds_ratio < 0.83 ~ "Moderate negative effect",
      TRUE ~ "Minimal effect"
    )
  ) %>%
  arrange(desc(odds_ratio)) %>%
  select(term, estimate, odds_ratio, interpretation) %>%
  head(10)
```

**How to read this table:**

-   $β_j$ (estimate) = change in log-odds of being happy for a 1 SD increase in predictor X_j (with all other predictors held constant at their means). To get the multiplicative effect on odds, calculate $exp(β_j)$.
-   **Odds ratio (**$exp(β_j)$) = how many times the **odds of being happy** multiply for a 1 SD increase

**Perfect separation and instability:**

-   If predictors can **perfectly separate the two classes**, e.g., all countries with `life_expectancy > 75` are happy, the model tries to assign **very large coefficients** to fit the data exactly.

-   Consequences:

    -   Coefficients become **extremely large**
    -   Predicted probabilities approach **0 or 1**
    -   Small changes in data can cause **large changes in predictions or coefficients**
    -   Predictions become **unstable** — tiny changes in the data could flip some results.

-   **Regularization** (penalized logistic regression) can shrink large coefficients and reduce instability.

**What about the model we trained?**

*The Problem: Perfect Separation and Unstable Coefficients*

In our model, some predictors like `drinking_water_pct` and `electricity_access_pct` show extremely large estimated coefficients (e.g., 9.03 and 7.35) and correspondingly huge odds ratios (over 8,300 and 1,560, respectively). This indicates the predictors almost perfectly separate the two outcome classes. Such **perfect separation** means that for certain predictor values, the model can classify the outcome with near certainty.

When logistic regression encounters perfect separation, it tries to fit the data exactly by inflating certain coefficients toward very large magnitudes. As a result, the model’s predictions become extreme (probabilities very close to 0 or 1). This causes instability: small changes or noise in the input data can lead to drastically different coefficients or predictions, reducing the model’s reliability and generalizability.

*The Proposed Fix: Regularization to Stabilize the Model*

To address this instability, **regularization** adds a penalty to the logistic regression's loss function, discouraging excessively large coefficients, and improving stability and generalization. The new loss function takes the form:

$$
\text{Loss} = \text{Logistic Loss} + \lambda \times \text{Penalty}
$$

where $\lambda$ controls the strength of penalization.

Different regularization penalties include:

-   **Lasso (L1)**: Adds $\lambda \sum |\beta_j|$. This can shrink some coefficients exactly to zero, effectively selecting important features and simplifying the model.
-   **Ridge (L2)**: Adds $\lambda \sum \beta_j^2$. This shrinks coefficients smoothly toward zero but rarely makes them exactly zero.
-   **Elastic Net**: Combines both L1 and L2 penalties, balancing feature selection and coefficient shrinkage.

### 4. Practical Implication

In essence, the current logistic regression is characterized by very large coefficients and odds ratios and could benefit from replacing the simple model with a **penalized logistic regression**. This approach controls coefficient magnitudes, stabilizes predictions, and enhances the model’s ability to generalize well to new data.

For more details, see:

-   ["FAQ What is complete or quasi-complete separation in logistic regression and what are some strategies to deal with the issue?"](https://stats.oarc.ucla.edu/other/mult-pkg/faq/general/faqwhat-is-complete-or-quasi-complete-separation-in-logistic-regression-and-what-are-some-strategies-to-deal-with-the-issue/)
-   ["Perfect separation, Quasi-perfect separation, and Linear dependence"](https://www.graphpad.com/guides/prism/latest/curve-fitting/reg_mult_logistic_dependence_and_separation.htm)
-   [Robert G. Clark, Wade Blanchard, Francis K.C. Hui, Ran Tian, Haruka Woods, *Dealing with complete separation and quasi-complete separation in logistic regression for linguistic data*, Research Methods in Applied Linguistics, Volume 2, Issue 1, 2023,100044, ISSN 2772-7661](https://www.sciencedirect.com/science/article/pii/S2772766123000046)
-   [*Separation and Convergence Issues in Logistic Regression*, Cornell Statistical Consulting Unit](https://cscu.cornell.edu/wp-content/uploads/lgsbias.pdf)
-   [Mohammad Ali Mansournia, Angelika Geroldinger, Sander Greenland, Georg Heinze, *Separation in Logistic Regression: Causes, Consequences, and Control*, American Journal of Epidemiology, Volume 187, Issue 4, April 2018, Pages 864–870](https://academic.oup.com/aje/article/187/4/864/4084405)
-   [Šinkovec H, Geroldinger A, Heinze G. *Bring More Data!-A Good Advice? Removing Separation in Logistic Regression by Increasing Sample Size*. Int J Environ Res Public Health. 2019;16(23):4658](https://pmc.ncbi.nlm.nih.gov/articles/PMC6926877/)

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

### 2.6.4 Predicting Probabilities

Logistic regression predicts the **probability of being happy (Yes)**:

$$
p = P(\textrm{is_happy = Yes} \mid X_1, \dots, X_p)
$$

-   Output is always between 0 and 1.
-   Log-odds can be converted to probability:

$$
p = \frac{e^{\textrm{log-odds}}}{1 + e^{\textrm{log-odds}}}
$$

-   Probabilities are **non-linear in log-odds**:

    -   Small changes near 0.5 → large probability change
    -   Near 0 or 1 → even large log-odds changes produce probabilities close to 0 or 1 (i.e., **saturate**)

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

### 2.6.5 Standardization

$$
X_\textrm{scaled} = \frac{X - \text{mean}(X)}{\text{SD}(X)}
$$

-   **SD = standard deviation**, shows spread of values
-   1 unit = 1 SD above or below mean

**Why normalize?**

-   Makes coefficients **comparable across variables on different scales**
-   Prevents optimization issues due to very large/small predictors
-   Ensures **interpretation**: β = change in log-odds per 1 SD increase

**Example revisited:**

-   Life expectancy SD = 7 years → β applies to a 7-year increase
-   GNI SD = 10,000 USD → β applies to 10,000 USD increase

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

### 2.6.6 Coefficients, Log-Odds, Odds, and Odds Ratios

We will use the **same simulated model for 2.6.6 and 2.6.7**, to make examples consistent:

$$
\textrm{log-odds} = \beta_0 + \beta_1 (\textrm{life_expectancy SD}) + \beta_2 (\textrm{GNI SD})
$$

-   Simulated coefficients: β₀ = 0, β₁ = 1.2, β₂ = 0.5
-   Standardized predictors (mean = 0, SD = 1)

**Step 1: Log-Odds**

$$
\textrm{log-odds} = \beta_0 + \beta_1 X_1 + \dots + \beta_p X_p
$$

-   $β$ = effect of predictor (1 SD increase) on log-odds

**Step 2: Odds**

-   **Odds = probability of success / probability of failure** = ( $\frac{p}{1-p}$ )
-   Complementary view: odds = $e^(log-odds)$
-   Example: p = 0.75 → odds = 0.75 / 0.25 = 3 → “3 happy per 1 unhappy”

**Step 3: Odds Ratio (OR)**

-   **Definition:** multiplicative change in odds per 1 SD increase in predictor
-   OR = $exp(β)$
-   OR \> 1 → increases odds
-   OR \< 1 → decreases odds
-   OR = 1 → no effect

**Illustrative Example Table (simulated model)**

| Predictor | $β$ (log-odds) | OR = $exp(β)$ | Interpretation |
|----|----|----|----|
| life_expectancy | 1.2 | 3.3 | 1 SD increase multiplies odds by 3.3× |
| gni_ppp | 0.5 | 1.65 | 1 SD increase multiplies odds by 1.65× |

> Notes:
>
> -   $β$ values are **simulated** for explanations purposes. They do not correspond to the $\beta$ coefficients our trained model estimated
> -   OR refers to **odds**, not probability.

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

### 2.6.7 From Log-Odds to Probability

We now use the **same simulated logistic regression model** as in 2.6.6.

**Step-by-step calculations:**

-   Fix GNI at mean (0), vary life expectancy ±1 or ±2 SDs

1.  **Log-odds**: log-odds = 0 + 1.2 × (SD units of life expectancy) + 0.5 × 0

2.  **Odds**: odds = $e^(\textrm{log-odds})$

3.  **Probability**: $p = \frac{\textrm{odds}}{(1 + \textrm{odds})}$

**Table: Life Expectancy effect**

| Life Expectancy (SD units) | Log-Odds | Odds  | Probability ($p$) |
|----------------------------|----------|-------|-------------------|
| -1                         | -1.2     | 0.30  | 0.23              |
| 0                          | 0        | 1     | 0.50              |
| +1                         | 1.2      | 3.32  | 0.77              |
| +2                         | 2.4      | 11.02 | 0.92              |

> “The probability saturate near extremes (i.e 0 and 1)” means:
>
> -   Large log-odds changes when $p$ is already near 0 or 1 → **small change in probability**
> -   Near $p$ = 0.5 → small log-odds change → **large probability change**

> Other predictors = 0 (mean in standardized units). These numbers are **illustrative**, not actual model coefficients.

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

## 2.7 Making Predictions

Let's try and make predictions on the test set based on our trained model (we are aware obviously that they might be unstable due to the perfect separability problem!)

We can use the same code to make predictions on the training set obviously (will be useful if we want to later on calculate training set metrics too and compare them with test set metrics to check for overfitting).

```{r make_predictions}
# Get predictions
predictions <- logistic_model %>%
  augment(happy_test_processed)

# Show different outputs
predictions %>%
  select(country, is_happy, .pred_class, .pred_No, .pred_Yes) %>%
  head(10)
```

Logistic regression outputs probabilities, not class labels.

-   `.pred_Yes` → probability that a country is happy

-   `.pred_No` → probability that a country is not happy

-   `.pred_class` → predicted class using a **threshold**

*What is the threshold?*

-   It's the probability $p$ (equivalent to `.pred_Yes`) above which the data point is classified as belonging to the `happy` class

-   Default = 0.5

-   `.pred_Yes >= 0.5` → classify as “Yes”

-   `.pred_Yes < 0.5` → classify as “No”

Note: `.pred_No + .pred_Yes = 1.0` (they're complementary!)

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

# Part 3: Evaluating Classification Models 📊

Now, we know how to construct and train a logistic regression model and get predictions out of it. It's time to evaluate the model!

## 3.1 The Confusion Matrix

First up is what we call a confusion matrix.

```{r confusion_matrix}
# Create confusion matrix
cm <- conf_mat(predictions, truth = is_happy, estimate = .pred_class)
print(cm)

# Visualize
autoplot(cm, type = "heatmap") +
  scale_fill_gradient(low = "white", high = "steelblue") +
  labs(title = "Confusion Matrix: Happy Country Classification",
       subtitle = "Rows = Actual, Columns = Predicted")
```

**Anatomy:**

|                | Predicted No        | Predicted Yes       |
|----------------|---------------------|---------------------|
| **Actual No**  | True Negative (TN)  | False Positive (FP) |
| **Actual Yes** | False Negative (FN) | True Positive (TP)  |

-   **TN**: Correctly predicted unhappy countries ✅
-   **TP**: Correctly predicted happy countries ✅
-   **FP**: Predicted happy but actually unhappy ❌ (Type I error)
-   **FN**: Predicted unhappy but actually happy ❌ (Type II error)

Whether **FP** or **FN** matter more depends on the problem at hand. For example:

-   in spam detection, **FP** matter more (the consequences of misclassifying a legitimate email as spam can be significant as it may result in important communications being filtered out and so the cost of a false positive is higher than that of a false negative in this case)
-   in disease detection (e.g cancer), **FN** matter more as a false negative means a patient that has his/her disease being misdiagnosed as being healthy will have his/her diagnosis deleayed and the appropriate treatment delayed (which will likely result in the deterioration of the patient's condition in the meanwhile or even in death!)

**What does the confusion matrix tell us about our model in this case?**

-   **18 countries** were correctly identified as *not happy* — the model recognized these correctly.
-   **13 countries** were correctly identified as *happy* — those are also correct predictions.
-   **2 countries** were *actually not happy*, but the model predicted them as *happy* — those are mistakes in the “happy” direction.
-   **7 countries** were *actually happy*, but the model predicted them as *not happy* — those are mistakes in the “unhappy” direction.

The model seems to:

-   do **well at recognizing “not happy” countries** (18 correct vs only 2 wrong).
-   be **less confident or accurate with “happy” countries**, since it misses several (7 wrong vs 13 correct).

In other words, it’s a bit **conservative** — it tends to predict “not happy” more often and misses some that are actually happy.

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

## 3.2 Classification Metrics

Now let's move on to some of the main classification metrics.

```{r classification_metrics}
# Calculate key metrics
metrics_class <- yardstick::metric_set(
  yardstick::accuracy,
  yardstick::precision,
  yardstick::recall,
  yardstick::f_meas
)
results_class <- predictions %>%
  metrics_class(truth = is_happy, estimate = .pred_class)

results_class
```

Let's understand each metric in the context of **predicting happy countries**:

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

### 3.2.1 Accuracy

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

**What it means:** Overall, what % of points did we classify correctly?

```{r explain_accuracy}
acc <- results_class %>% filter(.metric == "accuracy") %>% pull(.estimate)
cat(sprintf("Accuracy: %.1f%%\n", acc * 100))
cat(sprintf("We correctly classified %.1f%% of countries\n", acc * 100))
```

**⚠️ Accuracy can be misleading!**

It's only a good metric when the classes within your dataset are balanced. Can you guess what the problem might be in case of imbalance?

In our case, the data is pretty much balanced so accuracy is still acceptable.

> Accuracy of 77.5% means that out of all the countries in your dataset, the model correctly predicted whether each one was “happy” or “not happy” about 77.5% of the time.

So:

-   Roughly 3 out of every 4 countries were classified correctly.

-   The remaining 22.5% were misclassified — either a happy country labeled as not happy, or vice versa.

In plain terms:

> The model’s overall performance is reasonably good — it usually gets the classification right, but there’s still a meaningful number of mistakes to investigate.

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

### 3.2.2 Precision

$$\text{Precision} = \frac{TP}{TP + FP}$$

Precision is the proportion of true positives among all the predicted positive cases.

**What it means in our case:** Of all countries we predicted as "happy", what % were actually happy?

It is a good metric when the cost of false positives is high, e.g spam detection (the consequences of misclassifying a legitimate email as spam can be significant as it may result in important communications being filtered out and so the cost of a false positive is higher than that of a false negative in this case).

```{r explain_precision}
prec <- results_class %>% filter(.metric == "precision") %>% pull(.estimate)
cat(sprintf("Precision: %.1f%%\n", prec * 100))
cat(sprintf("When we say 'happy', we're right %.1f%% of the time\n", prec * 100))
```

**When precision matters in our case:**

-   Resource allocation (limited aid budgets - don't waste on countries that don't need it)
-   Targeted interventions (want to be sure you're helping the right countries)

**Trade-off:** High precision → fewer false positives BUT might miss some truly happy countries

**How can we interpret precision in our case?**

> Precision of 72.0% means that when the model predicts a country is “happy,” it’s correct 72% of the time.

So, among all the countries labeled as happy by the model:

-   72% truly are happy,

-   28% actually aren’t — those are the false positives.

In short:

> The model is fairly good at identifying happy countries, but it sometimes calls a country “happy” when it isn’t.

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

### 3.2.3 Recall (Sensitivity)

$$\text{Recall} = \frac{TP}{TP + FN}$$

Recall (or sensitivity) is the proportion of true positives among all actual positive cases.

**What it means in our case:** Of all countries that ARE happy, what % did we correctly identify?

```{r explain_recall}
rec <- results_class %>% filter(.metric == "recall") %>% pull(.estimate)
cat(sprintf("Recall: %.1f%%\n", rec * 100))
cat(sprintf("We found %.1f%% of all happy countries\n", rec * 100))
```

**When recall matters:**

-   Finding success stories to learn from
-   Comprehensive identification (don't miss any happy countries)
-   Complete coverage is critical

**Trade-off:** High recall → catch more cases BUT more false alarms

**How do we interpret recall in our context?**

> Recall of 90.0% means the model successfully identified 90% of all truly happy countries in the data.

So, out of every 10 countries that really are happy:

-   The model correctly predicted 9 of them as happy,

-   But missed 1, predicting it as not happy instead.

In plain language:

> The model is very good at finding happy countries, with only a few slipping through as false negatives.

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

### 3.2.4 F1 Score

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

**What it means:** Harmonic mean of precision and recall - balances both

```{r explain_f1}
f1 <- results_class %>% filter(.metric == "f_meas") %>% pull(.estimate)
cat(sprintf("F1 Score: %.3f\n", f1))
```

**When to use:** You care about both false positives AND false negatives

**How can we interpret this metric in our context?**

> F1 Score of 0.800 (or 80%) represents the balance between precision and recall for predicting “happy” countries.

So in plain terms:

-   The model does fairly well at both correctly identifying happy countries (high recall) and being right when it says a country is happy (reasonable precision).

-   The F1 score combines those two into a single number — showing overall effectiveness in identifying happiness without being too optimistic or too cautious.

In short:

> The model maintains a good balance between finding happy countries and being correct when it predicts happiness.

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

## 3.3 The Precision-Recall Tradeoff

```{r precision_recall_tradeoff, fig.width=10, fig.height=5}
# Calculate metrics at different thresholds
thresholds <- seq(0.1, 0.9, 0.05)

threshold_metrics <- map_df(thresholds, function(thresh) {
  preds_new <- predictions %>%
    mutate(.pred_class_new = factor(
      if_else(.pred_Yes >= thresh, "Yes", "No"),
      levels = c("No", "Yes")
    ))
  
  tibble(
    threshold = thresh,
    precision = precision_vec(preds_new$is_happy, preds_new$.pred_class_new),
    recall = recall_vec(preds_new$is_happy, preds_new$.pred_class_new),
    f1 = f_meas_vec(preds_new$is_happy, preds_new$.pred_class_new)
  )
})

threshold_metrics %>%
  pivot_longer(cols = c(precision, recall, f1), names_to = "metric") %>%
  ggplot(aes(x = threshold, y = value, color = metric)) +
  geom_line(size = 1.5) +
  geom_vline(xintercept = 0.5, linetype = "dashed") +
  labs(
    title = "Precision-Recall Tradeoff",
    subtitle = "Different thresholds optimize different metrics",
    x = "Classification Threshold",
    y = "Metric Value",
    color = "Metric"
  ) +
  scale_color_manual(values = c("precision" = "blue", "recall" = "red", "f1" = "purple"))
```

This graph illustrates the tradeoff between precision, recall, and F1 score as the classification threshold changes when converting predicted probabilities into class labels.

-   The x-axis is the **classification threshold**: the cutoff above which predicted probabilities are classified as positive (e.g., class = 1).
-   The y-axis is the **value of the metric** (precision, recall, or F1 score).
-   The blue line shows **precision**, which measures the accuracy of positive predictions (higher means fewer false positives).
-   The red line shows **recall**, which measures the ability to identify all positive cases (higher means fewer false negatives).
-   The purple line shows the **F1 score**, the harmonic mean of precision and recall, balancing both metrics.

**Key interpretation points:**

-   At low thresholds (\~0.1), recall is moderate (\~0.85) while precision is higher but slightly declining.
-   Around the threshold of 0.5 (marked by the dashed vertical line), recall holds steady at \~0.90 but precision dips sharply (\~0.70), causing the F1 score to drop as well.
-   Increasing the threshold beyond 0.6 increases recall sharply to 0.95 but precision stays near 0.7, slightly improving F1.
-   The curve shows how changing the threshold optimizes different metrics: higher thresholds improve recall but harm precision, while lower thresholds improve precision but reduce recall.

The takeaway is that the choice of classification threshold depends on which metric you prioritize. For example: - If finding most positives is critical (high recall), accept more false positives (lower precision) with a higher threshold. - If avoiding false positives is vital (high precision), choose a lower threshold, accepting lower recall.

The F1 score helps balance these competing goals, but depending on your use case, you may select a threshold that favors precision or recall more strongly.

This figure guides threshold selection for classification tasks by showing how performance metrics trade off as threshold varies.

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

# Part 4: Feature Selection 🎯

## 4.1 Why Feature Selection?

With our happiness data, we have \~20 predictors. Do we need all of them?

**Problems with too many features:**

1.  **Overfitting**: Model memorizes noise
2.  **Interpretability**: Hard to explain with 20 variables
3.  **Computational cost**: Slower training and prediction
4.  **Multicollinearity**: Correlated predictors (remember `rule_of_law` ↔ `gni_ppp`?)

**Solution:** Select only the most important features!

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

## 4.2 Method 1: Information Gain

**Idea:** Measure how much each feature reduces uncertainty about happiness

```{r information_gain, fig.width=10, fig.height=6}
# Calculate information gain
info_gain <- information_gain(
  formula = is_happy ~ .,
  data = happy_train_processed %>% select(-country, -iso3c, -region, -happiness_score),
  type = "infogain"
)

# Visualize
info_gain %>%
  arrange(desc(importance)) %>%
  slice(1:15) %>%  # Top 15 features
  ggplot(aes(x = reorder(attributes, importance), y = importance)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(
    title = "Feature Importance via Information Gain",
    subtitle = "Higher values = more predictive of happiness",
    x = "Feature",
    y = "Information Gain"
  )

# Show top 10
info_gain %>% 
  arrange(desc(importance)) %>% 
  head(10)
```

**Interpretation:**

-   `life_expectancy` and `physicians_per_1000` provide most information
-   Makes intuitive sense
-   Healthier countries tend to be happier!
-   Can set threshold (e.g., keep features with importance \> 0.05)

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

## 4.3 Method 2: Recursive Feature Elimination (RFE)

**Idea:**

1.  Build model with all features
2.  Remove least important feature
3.  Rebuild model
4.  Repeat until desired number of features

```{r rfe, warning=FALSE, message=FALSE}
# Set up RFE
rfe_control <- rfeControl(
  functions = lrFuncs,  # Logistic regression
  method = "cv",
  number = 5
)

# Prepare data for RFE
rfe_data <- happy_train_processed %>% 
  select(-country, -iso3c, -region, -happiness_score)

# Run RFE
set.seed(789)
rfe_results <- rfe(
  x = rfe_data %>% select(-is_happy),
  y = rfe_data$is_happy,
  sizes = c(3, 5, 7, 10),
  rfeControl = rfe_control
)

# View results
print(rfe_results)

# Best features
cat("\nOptimal number of features:", rfe_results$optsize, "\n")
cat("Selected features:", paste(rfe_results$optVariables, collapse = ", "), "\n")
```

**Advantages:**

-   Considers feature interactions
-   Finds optimal subset size automatically
-   Based on actual model performance (uses cross-validation - we'll see what this means next week!)

**Disadvantages:**

-   Computationally expensive (trains many models)
-   Can be slow with many features

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

## 4.4 Method 3: Sequential Forward Selection (SFS)

**Idea:** Start with nothing, add features one by one

```{r sfs, warning=FALSE, message=FALSE}
# Manual SFS implementation
available_features <- rfe_data %>% select(-is_happy) %>% names()
selected_features <- character(0)
best_accuracy <- 0

cat("Sequential Forward Selection:\n")
cat("============================\n\n")

for(i in 1:min(8, length(available_features))) {  # Try up to 8 features
  best_feature_this_round <- NULL
  best_acc_this_round <- 0
  
  for(feature in setdiff(available_features, selected_features)) {
    # Try adding this feature
    test_features <- c(selected_features, feature)
    formula_str <- paste("is_happy ~", paste(test_features, collapse = " + "))
    
    # Build model
    temp_model <- logistic_reg() %>%
      set_engine("glm") %>%
      set_mode("classification") %>%
      fit(as.formula(formula_str), data = rfe_data)
    
    # Evaluate on same training data (for simplicity)
    temp_pred <- temp_model %>% augment(rfe_data)
    temp_acc <- accuracy_vec(temp_pred$is_happy, temp_pred$.pred_class)
    
    if(temp_acc > best_acc_this_round) {
      best_acc_this_round <- temp_acc
      best_feature_this_round <- feature
    }
  }
  
  # Add best feature if it improves performance
  if(best_acc_this_round > best_accuracy) {
    selected_features <- c(selected_features, best_feature_this_round)
    best_accuracy <- best_acc_this_round
    cat("Step", i, "- Added:", best_feature_this_round, 
        "| Accuracy:", round(best_accuracy, 3), "\n")
  } else {
    cat("\nStopping - no improvement\n")
    break
  }
}

cat("\nFinal selected features:", paste(selected_features, collapse = ", "))
```

**Advantages:**

-   Intuitive and interpretable
-   Shows incremental value of each feature
-   Relatively fast

**Disadvantages:**

-   Greedy algorithm (might miss better combinations)
-   Doesn't consider removing features once added

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

## 4.5 Comparing Feature Selection Methods

```{r compare_methods}
# Build models with different feature sets

# Info gain top 5
ig_features <- info_gain %>% 
  arrange(desc(importance)) %>% 
  slice(1:5) %>% 
  pull(attributes)

model_ig <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification") %>%
  fit(as.formula(paste("is_happy ~", paste(ig_features, collapse = " + "))),
      data = rfe_data)

# RFE selected
model_rfe <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification") %>%
  fit(as.formula(paste("is_happy ~", paste(rfe_results$optVariables, collapse = " + "))),
      data = rfe_data)

# SFS selected
model_sfs <- logistic_reg() %>%
  set_engine("glm") %>%
  set_mode("classification") %>%
  fit(as.formula(paste("is_happy ~", paste(selected_features, collapse = " + "))),
      data = rfe_data)

# Full model for comparison
model_full <- logistic_model

# Evaluate all on test set
evaluate_classif <- function(model, test_data, name, n_features) {
  preds <- model %>% augment(test_data)
  tibble(
    Method = name,
    N_Features = n_features,
    Accuracy = accuracy_vec(preds$is_happy, preds$.pred_class),
    Precision = precision_vec(preds$is_happy, preds$.pred_class),
    Recall = recall_vec(preds$is_happy, preds$.pred_class),
    F1 = f_meas_vec(preds$is_happy, preds$.pred_class)
  )
}

comparison <- bind_rows(
  evaluate_classif(model_ig, happy_test_processed, "Information Gain (Top 5)", length(ig_features)),
  evaluate_classif(model_rfe, happy_test_processed, "RFE", length(rfe_results$optVariables)),
  evaluate_classif(model_sfs, happy_test_processed, "Sequential FS", length(selected_features)),
  evaluate_classif(model_full, happy_test_processed, "Full Model", 
                  ncol(rfe_data) - 1)
)

comparison %>%
  mutate(across(where(is.numeric) & !N_Features, ~round(., 3))) %>%
  knitr::kable()
```

**Key takeaways:**

-   **Simpler models often generalize better** (fewer features → less overfitting)
-   **Different methods select different features** (no single "right" answer)
-   **Consider both performance AND interpretability**
-   **Diminishing returns**: Adding more features beyond a point doesn't help much

**Which method to choose?**

Look at the performance metrics per feature selection model!

Also, think about you want a more interpretable model (i.e with fewer features) over a more complex one.

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

# 🎬 Wrap-Up & Key Takeaways

## What We Learned Today

**Part 1: Completing Regression Workflow**

-   ✅ Proper handling of missing data (imputation AFTER train/test split)
-   ✅ Recipe/prep/bake workflow prevents data leakage
-   ✅ Regularization (Ridge, Lasso, Elastic Net) helps with multicollinearity
-   ✅ Compare models on test set, not just training metrics

**Part 2: Classification Fundamentals**

-   ✅ Same data, different question type (continuous → binary)
-   ✅ Linear regression fails for categories (predictions escape \[0,1\])
-   ✅ Logistic regression uses sigmoid function to constrain probabilities
-   ✅ Interpret coefficients as odds ratios (multiplicative effects)

**Part 3: Classification Evaluation**

-   ✅ Confusion matrix shows four outcomes (TP, TN, FP, FN)
-   ✅ Accuracy can be misleading (especially with imbalanced data)
-   ✅ Precision: Of predicted positives, how many correct?
-   ✅ Recall: Of actual positives, how many found?
-   ✅ F1: Balances precision and recall (harmonic mean)
-   ✅ Precision-recall tradeoff: adjust threshold based on context

**Part 4: Feature Selection**

-   ✅ Three methods: Information Gain, RFE, Sequential FS
-   ✅ Simpler models often generalize better
-   ✅ Different methods may select different features
-   ✅ Balance performance with interpretability

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

## Decision Guide: Which Metric to Use?

| Scenario | Primary Metric | Why |
|------------------------|------------------------|------------------------|
| Balanced classes | **Accuracy** | Simple and interpretable |
| Imbalanced classes | **F1 or Precision/Recall** | Accuracy misleading |
| Cost of false positives high | **Precision** | Avoid false alarms |
| Cost of false negatives high | **Recall** | Don't miss true cases |
| Need balance | **F1** | Harmonic mean of both |

**For our happiness example:**

-   Accuracy works well (balanced 50/50 split)
-   But if we were predicting "crisis countries" (rare), recall would be more important!

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

## Next Week: More classification (new models and evaluation techniques/metrics)!

We'll cover:

-   **ROC curves and AUC** (comprehensive performance visualization)
-   **Cross-validation** (systematic hyperparameter tuning)
-   **Multi-class classification** (\>2 categories)
-   **Handling severe class imbalance** (SMOTE, class weights)

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

## 🏡 Practice Problems

1.  **Threshold Experimentation**:

    -   Try thresholds of 0.3, 0.5, 0.7 on the happiness data
    -   Calculate precision, recall, F1 for each
    -   Which would you choose and why?

2.  **Feature Engineering**:

    -   Create interaction terms (e.g., `life_expectancy * gni_ppp`)
    -   Create polynomial features (e.g., `gni_ppp^2`)
    -   Do these improve model performance?

3.  **Different Cutoffs**:

    -   Instead of median, try creating "is_happy" using:
        -   Top quartile (top 25% happiest)
        -   Bottom quartile (bottom 25%)
    -   How does class imbalance affect your metrics?

4.  **Compare All Methods**:

    -   Apply Ridge, Lasso, Elastic Net to the classification problem
    -   Compare with simple logistic regression (check both training and test set metrics!)
    -   Visualize how predictions differ

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

## Additional Resources

**Books:**

-   Hastie et al. (2021) - [*An Introduction to Statistical Learning*](https://www.statlearning.com/), 2nd Edition, Chapter 4

**Online:**

-   [Tidymodels documentation](https://www.tidymodels.org/)
-   [Feature Selection in R](https://topepo.github.io/caret/feature-selection-overview.html)
-   [Handling Class Imbalance](https://www.tidymodels.org/learn/models/class-imbalance/)

**Papers:**

-   Isabelle Guyon and André Elisseeff. 2003. *An introduction to variable and feature selection*. J. Mach. Learn. Res. 3, null (3/1/2003), 1157–1182.<https://dl.acm.org/doi/10.5555/944919.944968>

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

## The Complete Supervised Learning Journey

**Where we've been:**

1.  **Week 1-2**: Introduction to ML, linear regression fundamentals
2.  **Week 3** (today):
    -   Completed full regression pipeline with real data
    -   Introduced classification and logistic regression
    -   Learned proper evaluation metrics
    -   Applied feature selection methods

**Where we're going:**

-   **Week 4**: Cross-validation, hyperparameter tuning, ROC curves, KNN
-   **Week 5**: Tree-based models (Decision Trees, Random Forests, XGBoost/LightGBM) + SVM
-   **Week 6+**: Unsupervised learning and text mining

**You now know the fundamentals of supervised learning!** 🎉

**See you in next week!** 🚀
