Nieklasyczne metody statystyki

Regularyzacja

Twoja kolej!

Teraz nadszedł czas na przetestowanie tych metod (regresja grzbietowa i lasso) oraz metod oceny (zestaw walidacyjny, walidacja krzyżowa) na innych zbiorach danych. Możesz pracować z zespołem nad tą częścią laboratorium.

Możesz użyć dowolnego zbioru danych zawartego w ISLR lub wybrać jeden z pakietów danych na Kaggle/Data World itp. (zmienna zależna musi być ciągła).

Pobierz zbiór danych i spróbuj określić optymalny zestaw parametrów, które należy użyć do jego modelowania!

#install.packages('ISLR')
#install.packages('dplyr')
#install.packages('glmnet')
#install.packages('caret')
#install.packages('tidyverse')
library(ISLR)

## Warning: pakiet 'ISLR' został zbudowany w wersji R 4.3.2

library(dplyr)

## Warning: pakiet 'dplyr' został zbudowany w wersji R 4.3.2

## 
## Dołączanie pakietu: 'dplyr'

## Następujące obiekty zostały zakryte z 'package:stats':
## 
##     filter, lag

## Następujące obiekty zostały zakryte z 'package:base':
## 
##     intersect, setdiff, setequal, union

library(glmnet)

## Warning: pakiet 'glmnet' został zbudowany w wersji R 4.3.2

## Ładowanie wymaganego pakietu: Matrix

## Loaded glmnet 4.1-8

library(caret)    

## Warning: pakiet 'caret' został zbudowany w wersji R 4.3.3

## Ładowanie wymaganego pakietu: ggplot2

## Warning: pakiet 'ggplot2' został zbudowany w wersji R 4.3.3

## Ładowanie wymaganego pakietu: lattice

library(tidyverse)

## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.3.3

## Warning: pakiet 'stringr' został zbudowany w wersji R 4.3.2

## Warning: pakiet 'lubridate' został zbudowany w wersji R 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.3     ✔ tibble    3.2.1
## ✔ purrr     1.0.2     ✔ tidyr     1.3.0
## ✔ readr     2.1.4

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::expand() masks Matrix::expand()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ purrr::lift()   masks caret::lift()
## ✖ tidyr::pack()   masks Matrix::pack()
## ✖ tidyr::unpack() masks Matrix::unpack()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

head(Wage)

##        year age           maritl     race       education             region
## 231655 2006  18 1. Never Married 1. White    1. < HS Grad 2. Middle Atlantic
## 86582  2004  24 1. Never Married 1. White 4. College Grad 2. Middle Atlantic
## 161300 2003  45       2. Married 1. White 3. Some College 2. Middle Atlantic
## 155159 2003  43       2. Married 3. Asian 4. College Grad 2. Middle Atlantic
## 11443  2005  50      4. Divorced 1. White      2. HS Grad 2. Middle Atlantic
## 376662 2008  54       2. Married 1. White 4. College Grad 2. Middle Atlantic
##              jobclass         health health_ins  logwage      wage
## 231655  1. Industrial      1. <=Good      2. No 4.318063  75.04315
## 86582  2. Information 2. >=Very Good      2. No 4.255273  70.47602
## 161300  1. Industrial      1. <=Good     1. Yes 4.875061 130.98218
## 155159 2. Information 2. >=Very Good     1. Yes 5.041393 154.68529
## 11443  2. Information      1. <=Good     1. Yes 4.318063  75.04315
## 376662 2. Information 2. >=Very Good     1. Yes 4.845098 127.11574

data("Wage")
Wage = na.omit(Wage)

Przygotowanie danych

Zbiór danych, który analizujemy, pochodzi z pakietu ISLR i dotyczy wynagrodzeń pracowników. Zbiór ten zawiera różne cechy demograficzne i zawodowe pracowników, które mogą wpływać na wysokość wynagrodzenia.

Zmienne w analizie

• Zmienne zależna (y):

  • wage – Wynagrodzenie pracownika (wartość ciągła, główna zmienna, którą chcemy modelować).

• Zmienne niezależne (X):
  Zmiennymi predykcyjnymi są:

  1. age – Wiek pracownika.

  2. education – Poziom wykształcenia (kategoria nominalna, np. 1. < HS Grad, 2. HS Grad, 3. Some College, 4. College Grad, 5. Advanced Degree).

  3. jobclass – Klasa wykonywanej pracy (kategoria binarna: Industrial lub Information).

  4. health – Samoocena stanu zdrowia pracownika (<=Good lub >=Very Good).

  5. health_ins – Informacja o posiadaniu ubezpieczenia zdrowotnego (Yes lub No).

# Zmienna zależna
y <- Wage$wage
# Zmienne niezależne
X <- model.matrix(wage ~ age + education + jobclass + health + health_ins, data = Wage)[,-1]

Podział na zbiór treningowy i testowy

set.seed(123) 
train_index <- createDataPartition(y, p = 0.8, list = FALSE)
X_train <- X[train_index, ]
X_test <- X[-train_index, ]
y_train <- y[train_index]
y_test <- y[-train_index]

Dzielimy dane na zbiór treningowy i zbiór testowy, aby móc ocenić jakość modelu na danych, których wcześniej “nie widział”. Korzystamy z funkcji createDataPartition, aby losowo wybrać 80% obserwacji jako zbiór treningowy, a pozostałe 20% jako zbiór testowy. Następnie dzielimy macierz zmiennych niezależnych (X) oraz wektor zmiennej zależnej (y) na odpowiednie zestawy: treningowe (X_train, y_train) i testowe (X_test, y_test). Dzięki temu uczymy model na danych treningowych, a jego jakość weryfikujemy na danych testowych, co pozwala na obiektywną ocenę jego skuteczności.

Regresja OLS

W kolejnym kroku dopasowaliśmy model regresji liniowej OLS (ang. Ordinary Least Squares) na zbiorze treningowym. Celem było oszacowanie wpływu zmiennych niezależnych na zmienną zależną oraz ocena jakości predykcji modelu na zbiorze testowym. Nie ma uniwersalnej reguły, która pozwoliłaby nam stwierdzić na tym etapie który z modeli sprawdzi się lepiej ridge lub lasso. Spodziewamy się natomiast, że obie metody regularizacji poradzą sobie lepiej z ewentualną współliniowością i lepiej uogólnią wyniki w porównaniu z klasycznym modelem OLS. W praktyce:

# Dopasowanie modelu OLS
ols_model <- lm(wage ~ age + education + jobclass + health + health_ins, data = Wage[train_index, ])
summary(ols_model)

## 
## Call:
## lm(formula = wage ~ age + education + jobclass + health + health_ins, 
##     data = Wage[train_index, ])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -99.490 -19.656  -3.572  13.878 214.065 
## 
## Coefficients:
##                              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  67.19440    3.91866  17.147  < 2e-16 ***
## age                           0.48904    0.06278   7.789 9.95e-15 ***
## education2. HS Grad           7.94836    2.67757   2.968  0.00302 ** 
## education3. Some College     17.53744    2.83965   6.176 7.71e-10 ***
## education4. College Grad     31.22250    2.86542  10.896  < 2e-16 ***
## education5. Advanced Degree  54.09232    3.14541  17.197  < 2e-16 ***
## jobclass2. Information        4.14379    1.49159   2.778  0.00551 ** 
## health2. >=Very Good          8.75607    1.59218   5.499 4.21e-08 ***
## health_ins2. No             -18.39689    1.58634 -11.597  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 34.51 on 2393 degrees of freedom
## Multiple R-squared:  0.3105, Adjusted R-squared:  0.3082 
## F-statistic: 134.7 on 8 and 2393 DF,  p-value: < 2.2e-16

Wszystkie uwzględnione zmienne okazały się statystycznie istotne (p < 0,05), co oznacza, że wiek, poziom wykształcenia, rodzaj wykonywanej pracy, stan zdrowia oraz posiadanie ubezpieczenia zdrowotnego mają wpływ na wysokość płacy. Szczególnie duże znaczenie ma wykształcenie — Advanced Degree wiąże się z aż o około 54 USD/godz. wyższymi zarobkami w porównaniu do osób z wykształceniem niższym niż szkoła średnia. Równie istotny, choć negatywny wpływ, ma brak ubezpieczenia zdrowotnego, który obniża płacę o około 18 USD/godz. Z kolei każdy rok więcej w metryce odpowiada wzrostowi wynagrodzenia o ok. 0,49 USD/godz., a lepszy stan zdrowia (>= Very Good) przekłada się na dodatkowe +8,76 USD/godz. Pracownicy sektora „Information” zarabiają też średnio o 4,14 USD/godz. więcej niż osoby zatrudnione w „Industrial”.

Wskaźnik R2≈0,31 oznacza, że model wyjaśnia ok. 31% zmienności płacy, co w badaniach społeczno-ekonomicznych można uznać za umiarkowane dopasowanie (biorąc pod uwagę, że często istnieje wiele pominiętych czynników wpływających na zarobki). MSE = 1288,62 pokazuje skalę przeciętnego kwadratowego błędu predykcji, a Residual Standard Error = 34,51 wskazuje, o ile średnio pomyłka w prognozie odbiega od obserwowanej wartości wynagrodzenia. Choć model nie tłumaczy całej złożoności czynników kształtujących płace, jasno wskazuje, że wyższe wykształcenie, posiadanie ubezpieczenia zdrowotnego, lepszy stan zdrowia, praca w sektorze informacyjnym i wiek to kluczowe determinanty wyższych zarobków.

Regresja Ridge

# Przygotowanie danych dla glmnet
X_train_matrix <- as.matrix(X_train)
X_test_matrix <- as.matrix(X_test)

# Dopasowanie modelu Ridge (α = 0)
ridge_model <- cv.glmnet(X_train_matrix, y_train, alpha = 0)
plot(ridge_model)

ridge_lambda <- ridge_model$lambda.min
cat("Optimal Ridge Lambda:", ridge_lambda, "\n")

## Optimal Ridge Lambda: 1.580485

# Predykcje i ocena
ridge_pred <- predict(ridge_model, s = ridge_lambda, newx = X_test_matrix)
ridge_mse <- mean((y_test - ridge_pred)^2)
ridge_r2 <- 1 - (sum((y_test - ridge_pred)^2) / sum((y_test - mean(y_test))^2))

cat("Ridge MSE:", ridge_mse, "\nRidge R^2:", ridge_r2, "\n")

## Ridge MSE: 1295.187 
## Ridge R^2: 0.2880225

Z kolei, przechodząc do regresji grzbietowej (ridge), zastosowano walidację krzyżową w celu ustalenia optymalnej wartości parametru regularyzacji λ. Z przedstawionego wykresu widać, jak średni błąd kwadratowy (MSE) zmienia się w zależności od log⁡(λ): w okolicach log⁡(λ)≈0.46 (czyli λ≈1,58bserwujemy jego minimum. Ostatecznie, przy tym „najlepszym” λ Ridge MSE wynosi 1295,19, a R^2 to 0,288.

W porównaniu z modelem OLS, w którym MSE było na poziomie ok. 1288,62, a R2≈0,292R^2 , różnice są niewielkie. Ridge delikatnie ustępuje OLS pod względem błędu, ale w zamian zyskujemy większą stabilność i kontrolę nad współliniowością (zwłaszcza przy większej liczbie zmiennych). W tej sytuacji przewaga ridge nie jest mocno widoczna w samych miarach dopasowania, choć w praktyce może poprawiać odporność modelu na potencjalne zakłócenia czy wahania w danych.

Regresja Lasso

# Dopasowanie modelu Lasso (α = 1)
lasso_model <- cv.glmnet(X_train_matrix, y_train, alpha = 1)
plot(lasso_model)

# Optymalne lambda
lasso_lambda <- lasso_model$lambda.min
cat("Optimal Lasso Lambda:", lasso_lambda, "\n")

## Optimal Lasso Lambda: 0.02575796

# Predykcje i ocena
lasso_pred <- predict(lasso_model, s = lasso_lambda, newx = X_test_matrix)

lasso_mse <- mean((y_test - lasso_pred)^2)
lasso_r2 <- 1 - (sum((y_test - lasso_pred)^2) / sum((y_test - mean(y_test))^2))

cat("Lasso MSE:", lasso_mse, "\nLasso R^2:", lasso_r2, "\n")

## Lasso MSE: 1288.678 
## Lasso R^2: 0.2916008

Kontynuując analizę, w przypadku regresji Lasso zastosowano podobną procedurę walidacji krzyżowej, która wyłoniła optymalne λ≈0,0258 Dla tej wartości kary Lasso MSE wynosi 1288,68, a R2≈0,292. Wyniki są więc niemal identyczne jak w modelu OLS (gdzie MSE i R2R^2R2 miały odpowiednio wartości ~1288,62 i 0,2916). Lasso nie przyniosło w tym konkretnym zestawie danych istotnego zysku w predykcji, ale w przeciwieństwie do OLS potrafi zredukować znaczenie mniej istotnych zmiennych, „zerując” niektóre współczynniki. W efekcie może dawać prostszy, bardziej zinterpretowalny model kosztem zbliżonego lub lekko wyższego błędu. Jednak w naszej sytuacji różnica w stosunku do OLS jest minimalna, co pokazuje, że żaden z dodanych mechanizmów regularizacji (ridge czy lasso) nie zapewnia wyraźnej przewagi nad klasyczną regresją liniową w tym zestawie danych.

Ważność predyktorów (Lasso)

# Ważność zmiennych dla Lasso
lasso_coefficients <- coef(lasso_model, s = lasso_lambda)
lasso_coefficients_df <- data.frame(
  Variable = rownames(lasso_coefficients),
  Coefficient = as.vector(lasso_coefficients)
)
lasso_coefficients_df <- lasso_coefficients_df %>% filter(Coefficient != 0)

print(lasso_coefficients_df)

##                      Variable Coefficient
## 1                 (Intercept)  67.8501543
## 2                         age   0.4869979
## 3         education2. HS Grad   7.3559273
## 4    education3. Some College  16.9337841
## 5    education4. College Grad  30.6286741
## 6 education5. Advanced Degree  53.4917935
## 7      jobclass2. Information   4.1429969
## 8        health2. >=Very Good   8.7235583
## 9             health_ins2. No -18.4037668

W podsumowaniu porównań między OLS, ridge i lasso okazało się dość niespodziewanie, że model OLS osiągnął minimalnie lepsze wskaźniki dopasowania (np. niższe MSE) niż metody z regularizacją. Mimo to, ostatecznie zdecydowano się wybrać model Lasso. Głównym powodem jest jego potencjał do stabilnej selekcji zmiennych i ograniczania ewentualnego przeuczenia, chociaż w tym konkretnym zestawie danych żaden współczynnik nie został wyzerowany (czyli Lasso wybrało identyczny zestaw predyktorów, co OLS). Z punktu widzenia interpretacji najważniejsze wnioski pozostają spójne z wynikami modelu klasycznego: najwyższe dodatnie przełożenie na płacę zapewnia wyższe wykształcenie (szczególnie Advanced Degree), a największy negatywny wpływ ma brak ubezpieczenia zdrowotnego. Choć rezultaty były nieco zaskakujące (spodziewano się, że regularizacja poprawi wyniki względem OLS), Lasso nadal oferuje korzystne cechy (np. odporność na współliniowość i możliwość eliminacji mniej istotnych zmiennych w innych scenariuszach), dlatego ostatecznie wybrano je do dalszych analiz.