---
title: "Predicting National Happiness: Case Study"
subtitle: "Supervised Learning and Linear Regression in R"
author: "Dr Ghita Berrada"
format: html
self-contained: true
execute:
  echo: true
  warning: false
  error: false
---

## Goal of the Analysis

In this case study, we aim to **illustrate supervised learning concepts** using **linear regression** and its variants.

We will predict **national happiness scores** using socioeconomic and environmental indicators from the World Bank and World Happiness Report.

Key objectives:

1.  Introduce **data preprocessing** and exploration.\
2.  Explain **missingness mechanisms** and how to impute missing values.\
3.  Show **training/test split** and discuss **data leakage**.\
4.  Fit **linear regression models** (base R and tidymodels).\
5.  Evaluate models with **residual plots** and metrics (R², RMSE, MAE).\
6.  Apply **regularization**: Ridge, Lasso, Elastic Net.\
7.  Discuss **interpretation and ethics** of predictions.

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

## 1. Load Libraries and Data

```{r}
# Load libraries
library(tidyverse)
library(tidymodels)
library(vroom)
library(naniar)
library(visdat)
library(glmnet)
library(GGally) #make sure you update your ggplot2 to version 4.0.0 for this to work!
library(janitor)
library(countrycode)
library(missForest)
```

```{r}

setwd('..')
# Load cleaned World Bank data
wb_data <- read_csv("data/wb_macro_2020.csv")

# Load WHR 2020 happiness score
whr_data <- read_csv("data/whr_2020.csv")
whr_data <- whr_data |> clean_names()
```

Let's have a look at our data:

```{r}
glimpse(whr_data)
```

The only variable of interest to us here is `ladder_score` (which we'll rename to `happiness_score`): it's the happiness score per country and it is our outcome variable.

```{r}
whr_data <- whr_data |> rename('happiness_score'='ladder_score')
```

We won't be using the other variables in the dataset (except `country_name` to do the merge with the World Bank dataset): the values for the variables are mostly simply values that explain the contribution of said variable to the happiness score so they are not independent of the outcome variable. If we incorporated them, we'd risk data leakage.

::: callout-note
Data leakage occurs when information from outside the training dataset (we'll come back to this notion in a moment) is unintentionally used to create the model, giving it access to data that it wouldn’t have in a real-world scenario. This leads to overly optimistic performance estimates during training and testing, but the model fails when applied to new, unseen data.

In simpler terms: it’s when your model “cheats” by seeing answers it shouldn’t see.

**Key Points**

-   Source: Leakage usually happens when features in your training data contain information that is directly or indirectly correlated with the target variable in a way that wouldn’t be available during actual prediction.

-   Consequence: The model appears very accurate in testing but performs poorly in production.

*Example:*

-   Predicting whether a patient has a disease using lab tests, and accidentally including a feature like “medication prescribed” that only exists after diagnosis.

-   Using future data points to predict past events in time series forecasting.

**Types of Data Leakage**

-   Target leakage – When a feature includes information from the target variable.

-   Train-test contamination – When information from the test set leaks into the training set (e.g., scaling or feature selection done on the entire dataset before splitting).
:::

```{r}
glimpse(wb_data)
```

Let's check how many countries each of these datasets contain:

```{r}
whr_data |> select(country_name) |> unique() |> tally()
```

```{r}
wb_data |> select(country) |> unique() |> tally()
```

The World Bank dataset has many more countries than the World Happiness Report one! When we merge both those sets, we'll have to make sure to remove the countries that are not present in `whr_data`. There is another problem to tackle though before we merge the datasets: some countries are actually present in both datasets but are named differently (e.g "Egypt" in `whr_data` and "Egypt, Arab Rep." in `wb_data`). If we merged the datasets without any precaution, those countries would be dropped unnecessarily.

We have to standardize the country names before merging datasets:

```{r}
whr_final <- whr_data %>%
  mutate(
    country_name = trimws(country_name),  # remove any stray whitespace
    
    # 1️⃣ Manual mapping FIRST
    country_clean = case_when(
      country_name == "Czech Republic" ~ "Czechia",
      country_name == "Slovakia" ~ "Slovak Republic",
      country_name == "South Korea" ~ "Korea, Rep.",
      country_name == "Russia" ~ "Russian Federation",
      country_name == "Kyrgyzstan" ~ "Kyrgyz Republic",
      country_name == "Hong Kong S.A.R. of China" ~ "Hong Kong SAR, China",
      country_name == "Vietnam" ~ "Viet Nam",
      country_name == "Ivory Coast" ~ "Cote d'Ivoire",
      country_name == "Congo (Brazzaville)" ~ "Congo, Rep.",
      country_name == "Congo (Kinshasa)" ~ "Congo, Dem. Rep.",
      country_name == "Macedonia" ~ "North Macedonia",
      country_name == "Turkey" ~ "Turkiye",
      country_name == "Venezuela" ~ "Venezuela, RB",
      country_name == "Laos" ~ "Lao PDR",
      country_name == "Gambia" ~ "Gambia, The",
      country_name == "Iran" ~ "Iran, Islamic Rep.",
      country_name == "Palestinian Territories" ~ "West Bank and Gaza",
      country_name == "Swaziland" ~ "Eswatini",
      country_name == "Egypt" ~ "Egypt, Arab Rep.",
      country_name == "Yemen" ~ "Yemen, Rep.",
      
      # These are the missing ones
      country_name %in% c("Bosnia & Herzegovina", "Bosnia and Herzegovina") ~ "Bosnia and Herzegovina",
      country_name %in% c("Trinidad & Tobago", "Trinidad and Tobago") ~ "Trinidad and Tobago",
      country_name %in% c("Myanmar (Burma)", "Burma") ~ "Myanmar",
      country_name == "Kosovo" ~ "Kosovo",
      
      TRUE ~ NA_character_
    ),
    
    # 2️⃣ Fill missing manual mappings with standardized country names
    country_clean = if_else(
      is.na(country_clean),
      countrycode(country_name, "country.name", "country.name", warn = FALSE),
      country_clean
    )
  )
```

Now, let's merge our datasets:

```{r}
# Merge datasets
survey_results <- wb_data %>%
  inner_join(
    whr_final %>% select(country_clean, happiness_score),
    by = c("country" = "country_clean")
  )
```

**Why merge?**

We need the **outcome variable (happiness_score)** from WHR, while predictors come from World Bank. The merge aligns countries, ensuring that each row has predictors and the outcome.

> Effect: Countries in WB but missing in WHR will have `NA` for happiness_score, and vice versa.

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

## 2. Variable Overview

| Variable | Description | Source |
|----|----|----|
| `country` | Country name | World Bank/World Happiness Report |
| `iso3c` | ISO3 code | World Bank |
| `region` | Geographic region | World Bank |
| `income` | Income category (Low/Lower-Middle/Upper-Middle/High) | World Bank |
| `social_safety_net` | Adequacy of social safety net programs (% of total welfare of beneficiary households) | [World Bank](https://data.worldbank.org/indicator/per_sa_allsa.adq_pop_tot) |
| `cpi_gender_eq` | CPIA gender equality rating (1=low, 6=high) | [World Bank](https://data.worldbank.org/indicator/IQ.CPA.GNDR.XQ) |
| `rule_of_law` | Rule of Law estimate | [World Bank](https://data.worldbank.org/indicator/RL.EST) |
| `renewable_energy_pct` | Renewable energy consumption (% of total final energy) | [World Bank](https://data.worldbank.org/indicator/EG.FEC.RNEW.ZS) |
| `ghg_per_capita` | Per capita GHG emissions (Mt CO₂e / population) | World Bank (see below for details) |
| `electricity_access_pct` | Access to electricity (% of population) | [World Bank](https://data.worldbank.org/indicator/EG.ELC.ACCS.ZS) |
| `drinking_water_pct` | People using at least basic drinking water services (%) | [World Bank](https://data.worldbank.org/indicator/SH.H2O.BASW.ZS) |
| `education_upper_secondary_pct` | Population 25+ with at least upper secondary education (%) | [World Bank](https://data.worldbank.org/indicator/SE.SEC.CUAT.UP.ZS) |
| `migrant_stock_pct` | International migrant stock (% of population) | [World Bank](https://data.worldbank.org/indicator/SM.POP.TOTL.ZS) |
| `internet_users_pct` | Individuals using the Internet (%) | [World Bank](https://data.worldbank.org/indicator/IT.NET.USER.ZS) |
| `urban_population` | Number of urban residents | [World Bank](https://data360.worldbank.org/en/indicator/WB_WDI_SP_URB_TOTL_IN_ZS) |
| `suicide_rate` | Suicide mortality rate (per 100,000 population) | [World Bank](https://data360.worldbank.org/en/indicator/WB_WDI_SH_STA_SUIC_P5) |
| `youth_neet_pct` | Share of youth not in education, employment, or training (%) | [World Bank](https://data.worldbank.org/indicator/SL.UEM.NEET.ZS) |
| `population_density` | Population density (people per sq.km) | [World Bank](https://data.worldbank.org/indicator/EN.POP.DNST) |
| `gni_ppp` | GNI, PPP (current international \$) | [World Bank](https://data360.worldbank.org/en/indicator/WB_WDI_NY_GNP_MKTP_PP_CD) |
| `unemployment_pct` | Unemployment rate (%) | [World Bank](https://data360.worldbank.org/en/indicator/WB_WDI_SL_UEM_TOTL_ZS) |
| `physicians_per_1000` | Number of physicians per 1,000 people | [World Bank](https://data.worldbank.org/indicator/SH.MED.PHYS.ZSt) |
| `nurses_midwives_per_1000` | Number of nurses and midwives per 1,000 people | [World Bank](https://data.worldbank.org/indicator/SH.MED.NUMW.P3) |
| `oop_health_pct` | Out-of-pocket expenditure (% of health expenditure) | [World Bank](https://data.worldbank.org/indicator/SH.XPD.OOPC.CH.ZS) |
| `community_health_workers_per_1000` | Community health workers per 1,000 people | [World Bank](https://data.worldbank.org/indicator/SH.MED.CMHW.P3) |
| `gini_index` | Gini index of income inequality | [World Bank](https://data.worldbank.org/indicator/SI.POV.GINI) |
| `adult_literacy_pct` | Adult literacy rate (%) | [World Bank](https://data.worldbank.org/indicator/SE.ADT.LITR.ZS) |
| `life_expectancy` | Life expectancy at birth (years) | [World Bank](https://data360.worldbank.org/en/indicator/WB_WDI_SP_DYN_LE00_IN) |
| `maternal_death_risk_pct` | Lifetime risk of maternal death (%) | [World Bank](https://data.worldbank.org/indicator/SH.MMR.RISK.ZS) |
| `happiness_score` | Outcome: WHR happiness score | [World Happiness Report 2020](https://www.kaggle.com/datasets/londeen/world-happiness-report-2020) |

**Note:** `ghg_per_capita` was calculated as `ghg_total_mt / population_total` (i.e [Total greenhouse gas emissions including LULUCF (Mt CO2e)](https://data.worldbank.org/indicator/EN.GHG.ALL.LU.MT.CE.AR5) and [Population, total](https://data.worldbank.org/indicator/SP.POP.TOTL)). After this calculation, `ghg_total_mt` and `population_total` can be excluded from modeling to avoid redundancy.

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

## 3. Handling Missing Data

### Why Impute?

Missing data can bias results or reduce statistical power. **We must handle missingness** carefully.

### Missingness Mechanisms

1.  **MCAR (Missing Completely at Random)** Missingness does not depend on any variable. *Example:* A survey file is lost randomly for some countries.

2.  **MAR (Missing at Random)** Missingness depends on observed variables. *Example:* Low-income countries are less likely to report healthcare workforce data (`physicians_per_1000`) because weaker reporting systems exist — but this depends on `gni_ppp`, an observed variable.

3.  **MNAR (Missing Not at Random)** Missingness depends on the missing value itself. *Example:* Countries with extremely high unemployment may avoid reporting `unemployment_pct`.

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

### Imputation Using Random Forest

**Intuition:** Random Forest imputation predicts missing values using other variables as predictors. It iteratively fits trees on observed data and fills in missing entries.

Advantages:

-   Handles **non-linear relationships**.
-   Automatically considers interactions.
-   Preserves **distribution and variance** of original data better than mean/median imputation.

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

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

# Or a tidy summary
survey_results %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  tidyr::pivot_longer(everything(), names_to = "variable", values_to = "missing_count") %>%
  arrange(desc(missing_count))
```

```{r}
# Visualize missing data
gg_miss_var(survey_results)
```

```{r}
# Visualize patterns of missingness
vis_miss(survey_results)
```

```{r}
# Correlation of missingness
gg_miss_upset(survey_results)
```

For now, we'll drop the variables that have too many missing values i.e

```{r}
survey_final <- survey_results |> select(-community_health_workers_per_1000,-social_safety_net,-adult_literacy_pct,-cpi_gender_eq)
```

We won't impute straight away: we'll only do so after having split our data into training/test sets (data leakage!)

> Always perform imputation **after splitting into training/test sets** to avoid data leakage.

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

## 4. Exploratory Data Analysis (EDA)

```{r}
# Visualize missingness
vis_miss(survey_final)
```

```{r}
# Summary statistics
skimr::skim(survey_final)
```

```{r}
# Correlations
numeric_vars <- survey_final %>% select(c(-happiness_score,-year))
numeric_vars %>% 
  select(where(is.numeric)) %>%  # Keep only numeric columns
  GGally::ggpairs()
```

**Observations:**

Some very high correlations (multicollinearity):

-   `life_expectancy ↔ maternal_death_risk_pct`: `r = -0.762` (very strong negative)
-   `electricity_access_pct ↔ drinking_water_pct`: `r = 0.903` (redundant!)
-   `internet_users_pct ↔ gni_ppp`: `r = 0.870`
-   `nurses_midwives_per_1000 ↔ gni_ppp`: `r = 0.801`
-   `physicians_per_1000 ↔ nurses_midwives_per_1000`: `r = 0.870`
-   `internet_users_pct ↔ electricity_access_pct`: `r = 0.867`
-   `education_upper_secondary_pct ↔ internet_users_pct`: `r = 0.857`
-   `life_expectancy ↔ multiple variables`: 0.769, 0.790, 0.805, 0.869

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

## 5. Training/Test Split

```{r}
set.seed(123)
data_split <- initial_split(survey_final, prop=0.8, strata=region)
train_data <- training(data_split)
test_data <- testing(data_split)
```

`r strata=region` ensures that countries from each region are proportionally represented in both the training and test sets — which is often crucial when we suspect regional patterns influence the model.

**Why split?**

-   **Training set:** Fit the model.
-   **Test set:** Evaluate performance on unseen data.
-   Prevents **overfitting**.

**Data leakage:** Ensure no variable in training contains future information or derived from the outcome.

We can now safely impute the missing values (on the training/test set separately):

```{r}
# 1️⃣ Define a recipe for the survey data
survey_recipe <- recipe(happiness_score ~ ., data = train_data) %>%
  
  # 2️⃣ Remove columns that shouldn't be used for modeling
  #    (identifiers, country codes, region names, etc.)
  update_role(country, iso3c, region, new_role = "id") %>%
  
  # 3️⃣ Handle *new* factor levels in future/test data
  #    step_novel() ensures that unseen categories (like "Not classified")
  #    are safely converted to "novel" instead of causing an error
  step_novel(all_nominal_predictors(), new_level = "novel") %>%
  
  # 4️⃣ Handle missing factor values
  #    step_unknown() replaces NA values in factors with "unknown"
  step_unknown(all_nominal_predictors(), new_level = "unknown") %>%
  
  # 5️⃣ Collapse rare levels in training data into "other"
  #    This helps reduce noise from very infrequent categories
  step_other(all_nominal_predictors(), threshold = 0.01) %>%

  step_dummy(all_nominal_predictors()) %>%        # <<< convert factors to dummy variables
  
  # 6️⃣ Impute missing numeric values
  #    step_impute_bag() fills in missing numbers by finding patterns in other variables
  #    (no need to know the math now — just know it gives a smarter guess than the mean)
  step_impute_bag(all_numeric_predictors())

# 7️⃣ Prep the recipe (learn everything from training data)
survey_recipe_prep <- prep(survey_recipe, training = train_data)

# 8️⃣ Bake the recipe (apply to train and test data)
train_imputed <- bake(survey_recipe_prep, new_data = NULL)
test_imputed  <- bake(survey_recipe_prep, new_data = test_data)
```

### Simple intuition for `step_impute_bag`

> `step_impute_bag` **fills in missing numbers by looking at patterns in the other variables**. It’s smarter than just using the mean because it considers how columns relate to each other. You don’t need to understand the details yet (we'll see trees in a few weeks!) — just know it guesses missing values in a reasonable way.

::: callout-note
Here's an explanation of the `recipe`/`prep`/`bake` workflow in `tidymodels`. You'll see more of this in the coming weeks.

## **1️⃣ recipe() – define a preprocessing blueprint**

A **recipe** is like a **plan for preprocessing your data**. It doesn’t touch the data yet; it just specifies what you want to do.

Example:

``` r
library(recipes)

survey_recipe <- recipe(happiness_score ~ ., data = train_data) %>%
  update_role(country, iso3c, region, new_role = "id") %>% # don't use these in modeling
  step_impute_rf(all_numeric_predictors())               # impute numeric variables with Random Forest              
```

**Key points:**

-   `recipe()` defines the **target** (`happiness_score`) and **predictors**.

-   `step_*()` functions define transformations like:

    -   `step_impute_rf()`: Random Forest imputation
    -   `step_normalize()`: scale/center numeric data
    -   `step_dummy()`: convert categorical variables into dummy columns

-   Nothing is applied yet; it’s just a **blueprint**.

## **2️⃣ prep() – learn the transformations from training data**

`prep()` actually **trains or fits the recipe** on the dataset you provide.

``` r
survey_recipe_prep <- prep(survey_recipe, training = train_data)
```

What happens here:

-   For `step_impute_rf()`, it **builds a Random Forest** for each numeric column with missing values **using only training data**.
-   Other steps like scaling, encoding, etc., also **learn parameters** from the training set.

**Important:** This is where tidymodels prevents leakage. The test set is never used here.

## **3️⃣ bake() – apply the learned transformations**

`bake()` applies the **trained recipe** to a dataset.

``` r
# Apply to training data
train_imputed <- bake(survey_recipe_prep, new_data = NULL)

# Apply to test data
test_imputed <- bake(survey_recipe_prep, new_data = test_data)
```

-   `new_data = NULL` → bake on **training data**.
-   `new_data = test_data` → bake on **test data** using **same rules learned from training**.
-   This ensures the **same columns, same scaling, same imputations** are applied consistently.

### ✅ **Why this workflow is powerful**

1.  Keeps **training/test separation** intact — prevents leakage.
2.  Standardizes preprocessing for modeling pipelines.
3.  Supports complex sequences of steps (imputation, scaling, encoding, feature engineering) in **one object**.
4.  Integrates directly with `workflow()` and models in tidymodels.

💡 **Analogy:**

-   `recipe()` = your cooking **recipe** (ingredients & instructions)
-   `prep()` = **gathering ingredients and prepping them** (learning quantities from training data)
-   `bake()` = actually **cooking your meal** (applying the recipe to either training or test data)
:::

![](../figures/mermaid-drawing.png)

::: callout-note
-   `recipe()` → define the steps you want to do to your data.
-   `prep()` → calculate everything needed from the training data.
-   `bake()` → apply the transformations to any dataset (train or test).
:::

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

## 6. Linear Regression (Base R)

We're done with the preparation. We can now try and fit a linear regression model to our data. Let's first do it the way you're used to (the base R way). We'll use all variables except (`country`, `iso3c`, `year` and `region`) to do the prediction.

```{r}
# Fit model
lm_model <- lm(happiness_score ~  rule_of_law +
                 renewable_energy_pct + ghg_per_capita + electricity_access_pct +
                 drinking_water_pct + education_upper_secondary_pct +
                 migrant_stock_pct + internet_users_pct + urban_population +
                 suicide_rate + youth_neet_pct + population_density +
                 gni_ppp + unemployment_pct + physicians_per_1000 +
                 nurses_midwives_per_1000 + oop_health_pct + life_expectancy + maternal_death_risk_pct,
               data=train_imputed)

summary(lm_model)
```

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

### Residual Plots

Now that we have our model, it's time to evaluate it!🔎

```{r}
par(mfrow=c(2,2))
plot(lm_model)
```

**Explanation and interpretation of the residuals plots**

*1. Residuals vs Fitted*

**Purpose:** To check **linearity** and **homoscedasticity** (constant variance).

**Principle:** If the model fits well, residuals should:

-   Be randomly scattered around zero (no pattern)
-   Show roughly constant spread across fitted values

**Interpretation (from our plot):**

-   Shows relatively random scatter around zero, which is good
-   Observations 87, 51, and 72 are potential outliers
-   No strong systematic patterns, suggesting linearity assumption is reasonably met

*2. Normal Q–Q Plot*

**Purpose:** To assess whether residuals are **normally distributed**.

**Principle:**

-   Points should lie on the diagonal line if residuals are normal.
-   Deviations (especially in tails) indicate skewness or outliers.

**Interpretation (from our plot):**

-   Points generally follow the diagonal line well in the middle
-   Some deviation in the tails (especially observations 87, 51, 72)
-   Suggests mild departures from normality, but not severe

*3. Scale–Location Plot (Spread–Location Plot)*

**Purpose:** To check **homoscedasticity** again, but standardized for easier interpretation.

**Principle:**

-   Plots √\|standardized residuals\| vs fitted values.
-   The red line should be roughly horizontal for constant variance.

**Interpretation (from our plot):**

-   Relatively flat trend line with some scatter
-   No strong funnel pattern, indicating homoscedasticity (constant variance) is reasonably satisfied
-   Again, the same observations show larger standardized residuals

*4. Residuals vs Leverage (with Cook’s Distance)*

**Purpose:** To identify **influential observations** that may unduly affect the model.

**Principle:**

-   **Leverage** measures how far an observation’s predictors are from the mean of the predictors.
-   **Cook’s distance** shows combined influence on the fitted model.
-   Points beyond Cook’s contours (typically \> 0.5 or \> 1) are influential.

**Interpretation (from our plot):**

-   Most points have low leverage
-   Observations 91, 49, and 87 are noted
-   No points fall outside Cook's distance contours (0.5 and 1), meaning no highly influential outliers that are drastically affecting the model

### **Residual Diagnostics Summary**

The regression diagnostics suggest that the model fits the data reasonably well but with some areas for improvement. The **Residuals vs Fitted** plot shows a slight curvature, indicating that the relationship between predictors and the response variable may not be perfectly linear. The **Normal Q–Q** plot reveals minor deviations from normality at the tails, implying possible outliers or a slightly skewed error distribution. The **Scale–Location** plot displays a mild increase in residual spread as fitted values grow, suggesting **heteroscedasticity** (non-constant variance). Finally, the **Residuals vs Leverage** plot highlights a few observations (notably cases 51, 60, 87, 89, and 104) with relatively high leverage or Cook’s distance, meaning these data points have a larger-than-average influence on the regression results.

*Model Improvements*

1.  **Address nonlinearity:**

    -   We could consider adding **polynomial terms** (e.g., quadratic or cubic) or using **log/sqrt transformations** of predictors or the outcome.
    -   Alternatively, we could try a **nonlinear or generalized additive model (GAM)** if appropriate.

2.  **Check influential observations:**

    -   We could investigate the flagged points (51, 60, 87, etc.) to determine whether they are genuine or data-entry errors.
    -   We could consider refitting the model without them to assess their influence on coefficients.

3.  **Handle heteroscedasticity:**

    -   We could try **weighted least squares** or transform the response (e.g., log(Y)) to stabilize variance.

4.  **Normality of residuals:**

    -   If residuals are slightly non-normal but the sample size is large, this may not be a major issue.
    -   For small samples (like ours in fact), robust regression techniques could improve reliability.

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

## 7. Linear Regression (Tidymodels)

We can rewrite the linear regression model using `tidymodels`. You'll see more of tidymodels in the coming weeks.

```{r}
lm_spec <- linear_reg() %>%
  set_engine("lm") %>%
  set_mode("regression")

lm_fit <- lm_spec %>%
  fit(happiness_score ~ rule_of_law +
        renewable_energy_pct + ghg_per_capita + electricity_access_pct +
        drinking_water_pct + education_upper_secondary_pct +
        migrant_stock_pct + internet_users_pct + urban_population +
        suicide_rate + youth_neet_pct + population_density +
        gni_ppp + unemployment_pct + physicians_per_1000 +
        nurses_midwives_per_1000 + oop_health_pct + life_expectancy + maternal_death_risk_pct,
      data=train_imputed)

lm_fit
```

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

## 8. Model Evaluation

```{r}
# Predictions
preds <- predict(lm_fit, test_imputed) %>% bind_cols(test_imputed)

# Metrics

metrics <- metric_set(rmse, rsq, mae)
metrics(preds, truth=happiness_score, estimate=.pred)

```

**Regression Performance Metrics Interpretation**

| Metric | Meaning | Your Value | Interpretation |
|----|----|----|----|
| **RMSE (Root Mean Square Error)** | Measures the average magnitude of error (in the same units as the response variable). It penalizes large errors more heavily than small ones. | **0.746** | On average, predictions deviate from the observed values by about **0.75 units**. Lower RMSE indicates better fit. |
| **MAE (Mean Absolute Error)** | The mean of the absolute differences between predicted and observed values. It is more robust to outliers than RMSE. | **0.595** | The model’s predictions are off by roughly **0.6 units on average** — this confirms generally reasonable prediction errors. |
| **R² (Coefficient of Determination)** | Represents the proportion of variance in the dependent variable explained by the model. Ranges from 0 to 1. | **0.598** | The model explains about **60% of the variation** in the response — indicating an moderate fit. |

Given that our outcome ranges from 2.6 to 7.8 (a span of 5.2 units), an RMSE of 0.746 represents about 14% of the total range - this is reasonable but leaves room for improvement.

These were the test set metrics: let's compare them with the training set metrics:

```{r}
# Predictions
preds_train <- predict(lm_fit, train_imputed) %>% bind_cols(train_imputed)

# Metrics

metrics_train <- metric_set(rmse, rsq, mae)
metrics_train(preds_train, truth=happiness_score, estimate=.pred)

```

Both RMSE and MAE are lower for the training set (0.493 and 0.387 respectively). $R^2$ is higher (0.797). Not surprisingly the model performs better on the training set. However, the gap in performance between training/test set hints at possible model overfitting.

**Overall Interpretation**

The model is reasonably sound with moderate predictive power.

The assumptions appear largely satisfied, though we have a few observations worth investigating further (particularly 87, 51, 72). The lack of points beyond Cook's distance suggests these aren't severely distorting your results, but examining what makes these cases unusual could provide insights for model improvement.

**Avenues for improvement**

Here are several strategies to improve our model, roughly ordered from easiest to most complex:

1.  Investigate the Outliers:

We could start by examining observations 87, 51, and 72:

-   Are these data entry errors or measurement issues?
-   Do they represent a distinct subgroup or special cases?
-   If they're legitimate but unusual cases, we could consider whether they should be included or if a separate analysis is warranted
-   Sometimes understanding *why* these cases are unusual reveals missing predictors

2.  Feature Engineering

With 60% variance explained, there's substantial room to add information:

-   **Polynomial terms**: we could try squared or cubic terms for continuous predictors if we suspect non-linear relationships
-   **Interactions**: we could test whether effects of one variable depend on another (e.g., GNI × life expectancy, educational attainment x internet use, GNI x rule of law)
-   **Transformations**: Consider log, sqrt, or other transformations of skewed predictors (in particular if some variables have a vastly different scale than others)
-   **Domain knowledge**: What factors does theory suggest should matter but aren't in the model?

3.  Add More Predictors

-   We could look at what's currently unmeasured or unused in our dataset
-   The 40% unexplained variance has to come from somewhere

4.  Check for Non-linearity

Even though residuals look reasonable:

-   we could plot each predictor against residuals individually to spot non-linear patterns
-   we could consider using splines or generalized additive models (GAMs) for flexible non-linear fitting

5.  Try Different Model Types

If linear regression assumptions are limiting:

-   **Regularization** (Ridge, Lasso, Elastic Net): Helps if you have many predictors or multicollinearity
-   Non-linear models such as random forest (we'll see these in W05)

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

## 9. Regularization

Let's make sure `country` and other categorical variables are not part of our modeling:

```{r}
# After baking, remove ID columns
train_modeling <- train_imputed %>% select(-country, -iso3c, -region)
test_modeling <- test_imputed %>% select(-country, -iso3c, -region)
```

### Ridge Regression (L2)

```{r}
# Example: Ridge
ridge_spec <- linear_reg(penalty=0.1, mixture=0) %>%
  set_engine("glmnet") %>%
  set_mode("regression")

ridge_fit <- ridge_spec %>% fit(happiness_score ~ ., data=train_modeling)

### Lasso Regression (L1)

# Lasso Regression (mixture = 1)
lasso_spec <- linear_reg(penalty = 0.1, mixture = 1) %>%
  set_engine("glmnet") %>%
  set_mode("regression")

lasso_fit <- lasso_spec %>% fit(happiness_score ~ ., data=train_modeling)
```

### Elastic Net

```{r}
# ElasticNet Regression (mixture = 0.5)
elastic_spec <- linear_reg(penalty = 0.1, mixture = 0.5) %>%
  set_engine("glmnet") %>%
  set_mode("regression")

elastic_fit <- elastic_spec %>% fit(happiness_score ~ ., data=train_modeling)
```

Now that we've defined our regularized models, let's look at their performance:

```{r}
# -------------------------------
# 3️⃣ Predictions and Evaluation
# -------------------------------
ridge_preds <- predict(ridge_fit, new_data = test_modeling) %>%
  bind_cols(test_data %>% select(happiness_score))
lasso_preds <- predict(lasso_fit, new_data = test_modeling) %>%
  bind_cols(test_data %>% select(happiness_score))
elastic_preds <- predict(elastic_fit, new_data = test_modeling) %>%
  bind_cols(test_data %>% select(happiness_score))

# Function to compute metrics
get_metrics <- function(pred_df) {
  metrics(pred_df, truth = happiness_score, estimate = .pred)
}

ridge_metrics   <- get_metrics(ridge_preds)
lasso_metrics   <- get_metrics(lasso_preds)
elastic_metrics <- get_metrics(elastic_preds)

# Combine results
model_metrics <- bind_rows(
  ridge = ridge_metrics,
  lasso = lasso_metrics,
  elastic = elastic_metrics,
  .id = "model"
)

print(model_metrics)
```

**Interpretation:**

-   Ridge shrinks coefficients, reducing variance.
-   Lasso shrinks and can zero out coefficients → variable selection.
-   Elastic Net mixes both, good for correlated predictors.

**Evaluation**

Here's where we're at:

| Model           | RMSE     | R²       | MAE      |
|-----------------|----------|----------|----------|
| **Original LM** | 0.746    | 0.598    | 0.595    |
| **Ridge**       | 0.691 ✅ | 0.650 ✅ | 0.568 ✅ |
| **Lasso**       | 0.722 ✅ | 0.621 ✅ | 0.573 ✅ |
| **Elastic Net** | 0.685 ✅ | 0.657 ✅ | 0.543 ✅ |

*Key Wins:*

-   **Elastic Net is our best model**: R² improved from 59.8% → **65.7%** (explaining 6% more variance!)
-   **RMSE dropped** from 0.746 → **0.685** (8% reduction in prediction error)
-   **MAE dropped** from 0.595 → **0.543** (9% improvement)

The **Elastic Net** combines the best of both worlds (Ridge + Lasso), which is why it's outperforming.

Let's look at the regularized models coefficients.

```{r}

# ========================================
# Extract Coefficients from All Models
# ========================================

ridge_coefs <- tidy(ridge_fit) %>%
  filter(term != "(Intercept)") %>%
  mutate(model = "Ridge")

lasso_coefs <- tidy(lasso_fit) %>%
  filter(term != "(Intercept)") %>%
  mutate(model = "Lasso")

elastic_coefs <- tidy(elastic_fit) %>%
  filter(term != "(Intercept)") %>%
  mutate(model = "Elastic Net")

all_coefs <- bind_rows(ridge_coefs, lasso_coefs, elastic_coefs)

```

```{r}

# ========================================
# 1️⃣ Side-by-Side Dot Plot
# ========================================

ggplot(all_coefs, aes(x = estimate, y = reorder(term, estimate), color = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_point(size = 2.5, alpha = 0.7) +
  facet_wrap(~model, ncol = 3) +
  labs(
    title = "Coefficient Comparison Across Regularization Methods",
    x = "Coefficient Estimate",
    y = "Predictor"
  ) +
  theme_minimal() +
  theme(
    legend.position = "none",
    axis.text.y = element_text(size = 8)
  )
```

```{r}

# ========================================
# 2️⃣ Top 15 Features Comparison
# ========================================

# Identify top 15 by average absolute coefficient
top_vars <- all_coefs %>%
  group_by(term) %>%
  summarise(avg_abs = mean(abs(estimate))) %>%
  arrange(desc(avg_abs)) %>%
  head(15) %>%
  pull(term)

top_coefs <- all_coefs %>%
  filter(term %in% top_vars)

ggplot(top_coefs, aes(x = estimate, y = reorder(term, estimate), color = model)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_point(size = 3, alpha = 0.8, position = position_dodge(width = 0.5)) +
  labs(
    title = "Top 15 Most Important Predictors",
    subtitle = "Comparing coefficients across Ridge, Lasso, and Elastic Net",
    x = "Coefficient Estimate",
    y = "Predictor",
    color = "Model"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

```

### 🔍 Key insights from the coefficient plots

**1. The big winners (Strongest positive effects on happiness):**

-   **`ghg_per_capita`** (\~10,000) - Surprisingly massive coefficient! This is unusual and might indicate:

    -   Highly correlated with wealth/development
    -   Possible scaling issue (this variable wasn't normalized so that might be the problem)
    -   Could be capturing industrial development

-   **`rule_of_law`** (\~7,500) - Strong governance = much happier populations

-   **`maternal_death_risk_pct`** (\~5,000) - Lower maternal mortality = higher happiness

-   **`nurses_midwives_per_1000`** & **`gini_index`** - Healthcare access and income equality matter

**2. The big detractors (negative effects):**

-   **`income_Lower.middle.income`** & **`income_Low.income`** - Being in lower income categories strongly reduces happiness (as expected)
-   **`unemployment_pct`** & **`youth_neet_pct`** - Joblessness hurts happiness

**3. Model Behavior Differences:**

**Lasso vs Ridge vs Elastic Net:**

-   All three models largely **agree** on the top predictors - good sign of stability!
-   **Lasso** (red) tends to have slightly smaller coefficients, consistent with its aggressive shrinkage
-   **Ridge** (blue) keeps coefficients similar to Elastic Net
-   **Elastic Net** (green) sits nicely between them

The fact that all three models identify the same top features suggests these are **robust, reliable predictors**.

### ⚠️ Red Flags to Address:

**The `ghg_per_capita` Coefficient is Suspiciously Large**

A coefficient of 10,000 suggests:

1.  **Scaling issue** - The magnitude is way larger than others. Hints at the fact the variable wasn't normalized
2.  **Multicollinearity** - Likely highly correlated with another variable (`gni_ppp` (wealth)?)

Going to our earlier recipe, we didn't use `step_normalize()` in there! This is **critical** for regularized regression. Without it, variables with different scales get penalized unfairly.

**Here's what we should add:**

``` r
survey_recipe <- recipe(happiness_score ~ ., data = train_data) %>%
  # ... our existing steps ...
  step_impute_bag(all_numeric_predictors()) %>%
  step_normalize(all_numeric_predictors())  # ← ADD THIS!
```

Then we'd refit our models. The coefficients will be on the same scale and much easier to interpret.

Once we fix the scaling, we'll be able to confidently say:

-   **Governance** (`rule_of_law`) is crucial for happiness
-   **Healthcare access** (nurses, maternal health) matters significantly\
-   **Economic status** (income categories) has huge impact
-   **Youth employment** is important for societal happiness

These are meaningful, interpretable results that align well with happiness research!

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

## 10. Ethical Considerations

-   Predicting national happiness may **oversimplify complex sociocultural factors**.
-   Always interpret models **cautiously**, and be aware of **cultural biases** in survey responses.

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

## Next Steps

You can try

-   **feature selection** (have a look at sequential feature selection, recursive feature elimination or feature selection through information gain - you might want to look at the `colino` library for some of these),
-   **interaction terms**,
-   or **log-transform skewed variables**. And normalization too!
