Raport 5 - Regularyzacja

Główne pytania do Raportu 5:

Raport 5

  1. Do analizy wybrałam zbiór Credit z pakietu ISLR.

  2. Zmienna objaśniana w modelach to dochód (Income).

  3. Hipoteza początkowa: Oczekuje, że lepszą okaże się być któraś z metod nieklasycznych tzn. metoda lasso, regresja grzbietowa lub RSS (Regularyzacja Elastic Net) w porównaniu do regresji liniowej.

biblioteki

library(ISLR)
library(glmnet)
library(dplyr)
library(tidyr)
library(rsample)
library(tibble)
library(ggplot2)
library(rlang)
library(recipes)
library(Matrix)

Spis treści

  1. opis danych

  2. Model liniowy

  3. Model regresji grzbietowej

  4. Model metodą lasso

  5. Model metodą RSS

  6. Porównanie modeli

  7. Wnioski

Opis danych

Do modelowania zostanie wykorzystany publicznie dostępny zbiór danych - credit, zawierający informacje o zadłużeniu na karcie kredytowej dla 400 klientów.

Zmienne modelu

Zmienna objaśniana

income (in thousands of dollars) - dochód w tysiącach dolarów

Zmienne objaśniające

  • age - wiek w latach

  • cards (number of credit cards) - liczba kart kredytowych

  • education (years of education) - liczba lat edukacji

  • Limit (credit limit) - limit kredytowy Rating (credit rating)

  • Gender - Male/Female

  • Student - TAK/NIE

  • Ethnicity (3 poziomy) - Caucasian/African American/Asian

  • Balance - saldo na karcie

Oszacowania dla 4 metod

0. OLS

Poniżej przedstawiono oszacowany model regresji liniowej oraz krótka charakterystyka Dochodów (histogram, wykresy)

## 
## Call:
## lm(formula = Income ~ ., data = Credit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.974  -7.719  -3.341   5.191  34.294 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -52.033752   4.108162 -12.666  < 2e-16 ***
## ID                   0.004828   0.004789   1.008  0.31402    
## Limit                0.018842   0.003647   5.167 3.83e-07 ***
## Rating               0.143524   0.054043   2.656  0.00824 ** 
## Cards                1.321744   0.484734   2.727  0.00669 ** 
## Age                 -0.007460   0.032671  -0.228  0.81950    
## Education           -0.139600   0.176236  -0.792  0.42878    
## GenderFemale        -1.272507   1.100628  -1.156  0.24833    
## StudentYes          41.268951   2.171062  19.009  < 2e-16 ***
## MarriedYes          -0.519615   1.144256  -0.454  0.65001    
## EthnicityAsian       1.728813   1.557888   1.110  0.26781    
## EthnicityCaucasian   0.595004   1.350901   0.440  0.65986    
## Balance             -0.094978   0.002850 -33.321  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.9 on 387 degrees of freedom
## Multiple R-squared:  0.9073, Adjusted R-squared:  0.9044 
## F-statistic: 315.5 on 12 and 387 DF,  p-value: < 2.2e-16
## 
## Call:
## lm(formula = Income ~ Limit + Rating + Cards + Student + Balance, 
##     data = Credit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -17.440  -7.805  -3.212   5.151  34.635 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -53.373691   2.431849 -21.948  < 2e-16 ***
## Limit         0.019027   0.003601   5.285 2.09e-07 ***
## Rating        0.139953   0.053288   2.626  0.00897 ** 
## Cards         1.343453   0.480853   2.794  0.00546 ** 
## StudentYes   41.049813   2.142737  19.158  < 2e-16 ***
## Balance      -0.094740   0.002785 -34.024  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.87 on 394 degrees of freedom
## Multiple R-squared:  0.9062, Adjusted R-squared:  0.905 
## F-statistic: 760.9 on 5 and 394 DF,  p-value: < 2.2e-16
## 
##  Shapiro-Wilk normality test
## 
## data:  Credit$Income
## W = 0.80965, p-value < 2.2e-16

W powyższym modelu wykorzystano wszystkie 12 zmiennych objaśniających - 5 okazało się być istotnych statystycznie.

R-squared wynosi 90,73% - co sugeruje dobre dopasowanie zmiennych teoretycznych do zmiennych rzeczywistych.

W klasycznym podejściu, aby w najprostrzy sposów wyelimować nie istotne zmienne w modelu liniowym i stworzyć model składający się tylko zmiennych isotnych statystycznie, stostuje się usuwananie z modelu po jednej zmiennej w kolejności od najwyższego p-value (test pominiętych zmiennych). Potrafi być to czasochłonne.

W porónaniu, metoda Lasso pozwala na szybszą weryfikacje istotności statystycznej zmiennych zależnych, co jest jej zaletą.

1. metoda grzbietowa

Definiowanie zmiennych w modelu

Zdefiniowano zmienną objaśnianą (y=wynagrodzenia) i zmienne obiaśniające (wszystkie pozostałe)

grid = 10^seq(10, -2, length = 100)
ridge_mod = glmnet(x, y, alpha = 0, lambda = grid)
dim(coef(ridge_mod))
## [1]  14 100
plot(ridge_mod, main="Wynagrodzenia~.")    # wykres współczynników

ridge_mod1 = glmnet(x1, y, alpha = 0, lambda = grid)
dim(coef(ridge_mod1, main="Wynagrodzenia~."))
## [1]  12 100
plot(ridge_mod1, main="Wynagrodzenia~wszystkie (bez czy_student)") #ciekawszy wykres

Dla modelu z dwunastoma zmiennymi zależnymi dostrzegamy, że jedna zmienna istotnie odchyla się od reszty (wykres “Wynagrodzenia~.”). Jest to zmienna CzyStudent. Jest to już pierwszy sygnał, że metoda grzbietowa może nie dać oczekiwanych pozytywnych rezultatów.

Po usunięciu tylko tej zmiennej z modelu wykres współczynników względem lambdy nie daje jednoznacznego wyniku gdyż dla lambda>4 większość współczynników zmiennych (poza dwoma) zmniejsza swoją wartość w kierunku zera. Natomiast wspomniane dwie zmienne (linia ciemno i jasno niebieska na wykresie “Wynagrodzenia~wszystkie pozostałe zmienne(bez czy_student)” ) oddalają się od zera do +nieskończoności.

Oszacowanie parametrów modelu metodą grzbietową

ridge_mod$lambda[4] # Wyświetl 4tą wartość lambdy
## [1] 4328761281
coef(ridge_mod)[,4] # Wyświetl współczynniki związane z 4-tą wartością lambdy
##        (Intercept)        (Intercept)                 ID              Limit 
##       4.521888e+01       0.000000e+00       9.222134e-11       9.834909e-11 
##             Rating              Cards                Age          Education 
##       1.465871e-09      -3.818968e-09       2.913149e-09      -2.539482e-09 
##       GenderFemale         StudentYes         MarriedYes     EthnicityAsian 
##      -6.150725e-09       1.873151e-08       2.094719e-08      -1.125396e-08 
## EthnicityCaucasian            Balance 
##      -1.127823e-08       2.890255e-10
sqrt(sum(coef(ridge_mod)[-1,4]^2)) # Oblicz normę l2
## [1] 3.336324e-08

Podzielimy teraz próbki na zbiór treningwy i testowy (80:20), by oszacować błąd testu reogresji grzbietowej i lasso.

set.seed(1)

train = Credit %>%
  sample_frac(0.8)

test = Credit %>%
  setdiff(train)

x_train = model.matrix(Income~., train)[,-1]
x_test = model.matrix(Income~., test)[,-1]

y_train = train %>%
  select(Income) %>%
  unlist() %>%
  as.numeric()

y_test = test %>%
  select(Income) %>%
  unlist() %>%
  as.numeric()

Dopasowujemy model regresji grzbietowej na zbiorze treningowym i oceniamy jego MSE na zbiorze testowym, używając \(\lambda = 1\). Wybrano taką lambde w oparciu o powyższy wykres współczynników.

## [1] 167.2863

Testowe MSE dla lambda =1 wynosi 167.2863.

Porównanie współczynników regresji liniowej i regresji grzbietowej (lambda=1)

## [1] 167.2863
## 
## Call:
## lm(formula = Income ~ ., data = train)
## 
## Coefficients:
##        (Intercept)                  ID               Limit              Rating  
##         -54.714721            0.001461            0.019307            0.136947  
##              Cards                 Age           Education        GenderFemale  
##           1.592745            0.011159           -0.042706           -1.704036  
##         StudentYes          MarriedYes      EthnicityAsian  EthnicityCaucasian  
##          39.878980           -0.782886            1.873665            0.202290  
##            Balance  
##          -0.093616
##        (Intercept)                 ID              Limit             Rating 
##      -47.545585970        0.001615479        0.012173692        0.178664295 
##              Cards                Age          Education       GenderFemale 
##        0.732122387        0.075317175       -0.067209435       -1.380117333 
##         StudentYes         MarriedYes     EthnicityAsian EthnicityCaucasian 
##       29.915957471       -0.734158442        1.969257581        0.157966524

Wygląda na to, że rzeczywiście poprawiamy się w stosunku do zwykłego modelu liniowego, ponieważ współczynniki regresji grzbietwej są bliższe zeru dla prawie każdej zmiennej. Wyższe wspołczynniki tylko dla rating, Age i Ethincity.

Minimalizacja lambdy

set.seed(1)
cv.out = cv.glmnet(x_train, y_train, alpha = 0) # Dopasuj model regresji grzbietowej na danych treningowych
bestlam = cv.out$lambda.min  # Wybierz lamdę, która minimalizuje treningowy MSE 
bestlam
## [1] 2.861825

Widzimy zatem, że wartość \(\lambda\), która powoduje najmniejszy błąd dla danych treningowych walidacji krzyżowej to 2,86

Poniżej funkcja MSE i logarytmu z \(\lambda\) dla zbioru danych treningowych:

MSE dla najlepszej lambdy (2,86)(dane trningowe) wynosi 252.0325. Co pokazuje, że błąd jest wyższy niż dla lambda=1 (167.2863) wybranego w oparciu o wykres współczynników i od modelu liniowego (MSE=135).

Zatem, widzimy, że należało by jeszcze zmniejszyć lambda, do np. 0.05.

Uzyskano wymik MSE=135.5891. Zatem dla modelu przy 12 zmiennych, tylko bardzo mała wartość alfy, daje wystarcząco niski średni błąd standardowy.

## [1] 135.5891

2. metoda lasso

Wykres współczynników modelu względem lambdy przy zastosowaniu funkcji glmnet(), alpha=1

Na powyższym wykresie, że dla labda <31 współczynniki z wyjątkiem jednego są bliskie lub rónwne zero. Dla lambda>31 pozostałe współczynniki zwiększają się.

Dopasowanie modelu lasso

Przy pomocy funkcji minimalizującej błąd MSE wyznaczamy wartość najlepszego lambda.

Na tej podstawie ustalono, że najlepsza labmda=0.118251, dla najmniejszego MSE= 136.3699. Co w porównaniu do regresji grzebietowej daje mniejszy błąd MSE.

Wykres lambda względem MSE poniżej.

set.seed(1)
cv.out = cv.glmnet(x_train, y_train, alpha = 1) # Dopasuj model lasso do danych treningowych
plot(cv.out) # Narysuj wykres MSE dla próby uczącej jako funkcję lambda

bestlam = cv.out$lambda.min # Wybierz lamdę, która minimalizuje MSE w próbie uczącej
lasso_pred = predict(lasso_mod, s = bestlam, newx = x_test) # Użyj najlepszej lambdy do przewidywania danych testowych
bestlam
## [1] 0.118251
mean((lasso_pred - y_test)^2) # Oblicz MSE w próbie testowej
## [1] 136.3699

MSE dla wybranej alfy jest niskie. W porównaniu do modelu najmniejszych kwadratów, i regresji grzbietowej z \(\lambda\) wybranej przez walidację krzyżową, metoda lasso dała bardzo podobne rezultaty.

Lasso - oszacowanie parametrów

dla lambda=grid

##        (Intercept)        (Intercept)                 ID              Limit 
##      -51.116242566        0.000000000        0.003875661        0.016690854 
##             Rating              Cards                Age          Education 
##        0.165799108        1.027643462        0.000000000       -0.098168397 
##       GenderFemale         StudentYes         MarriedYes     EthnicityAsian 
##       -1.048059880       39.372118321       -0.229978110        1.054837252 
## EthnicityCaucasian            Balance 
##        0.014905780       -0.091499584

Jednakże lasso ma istotną przewagę nad regresją grzbietową w tym, że wynikowe oszacowania współczynników zapewniają także inforamcje o ich ostotności statystycznej.

Wybierając tylko predyktory o niezerowych współczynnikach widzimy, że model lasso z \(\lambda\) wybranym przez walidację krzyżową ma 11 zmiennych istotnych, na 12 zmiennych w bazowym ujęciu. Odrzucamy zmienną Age.

Dodatkowo dla lambda=bestlam=0.118251

##        (Intercept)        (Intercept)                 ID              Limit 
##      -50.378723828        0.000000000        0.003842671        0.018903719 
##             Rating              Cards                Age          Education 
##        0.133383405        1.191319854        0.000000000       -0.106405240 
##       GenderFemale         StudentYes         MarriedYes     EthnicityAsian 
##       -1.047568616       39.569480696       -0.157843713        0.992241397 
## EthnicityCaucasian            Balance 
##        0.000000000       -0.091727651

Wyświetlanie tylko niezerowych współczynników

##        (Intercept)                 ID              Limit             Rating 
##      -51.116242566        0.003875661        0.016690854        0.165799108 
##              Cards          Education       GenderFemale         StudentYes 
##        1.027643462       -0.098168397       -1.048059880       39.372118321 
##         MarriedYes     EthnicityAsian EthnicityCaucasian            Balance 
##       -0.229978110        1.054837252        0.014905780       -0.091499584

Dla porównania Model Lasso z wykorzystaniem tidyflow

## Warning: `pull_workflow_fit()` was deprecated in workflows 0.2.3.
## ℹ Please use `extract_fit_parsnip()` instead.
## ℹ The deprecated feature was likely used in the tidyflow package.
##   Please report the issue at <https://github.com/cimentadaj/tidyflow/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## ══ Tidyflow [trained] ══════════════════════════════════════════════════════════
## Data: 400 rows x 12 columns
## Split: initial_split w/ default args
## Formula: Income ~ .
## Resample: None
## Grid: None
## Model:
## Linear Regression Model Specification (regression)
## 
## Main Arguments:
##   penalty = 0.001
##   mixture = 1
## 
## Computational engine: glmnet 
## 
## ══ Results ═════════════════════════════════════════════════════════════════════
## 
## 
## Fitted model:
## 
## Call:  glmnet::glmnet(x = maybe_matrix(x), y = y, family = "gaussian",      alpha = ~1) 
## 
##    Df  %Dev  Lambda
## 1   0  0.00 27.4500
## 
## ...
## and 72 more lines.

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Formula
## Model: linear_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## Income ~ .
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## 
## Call:  stats::glm(formula = ..y ~ ., family = stats::gaussian, data = data)
## 
## Coefficients:
##        (Intercept)                  ID               Limit              Rating  
##         -52.945709            0.002519            0.017832            0.161970  
##              Cards                 Age           Education        GenderFemale  
##           1.304157            0.011654           -0.198533           -1.315769  
##         StudentYes          MarriedYes      EthnicityAsian  EthnicityCaucasian  
##          41.642923           -0.942324            2.345405            0.385694  
##            Balance  
##          -0.095553  
## 
## Degrees of Freedom: 299 Total (i.e. Null);  287 Residual
## Null Deviance:       376300 
## Residual Deviance: 33800     AIC: 2297

4. Porównanie

Porównanie współczynników z modelu regresji liniowej, grzbietowej i lasso.

coefs Linear.coefficients Ridge.coefficients lasso.coefficients
(Intercept) -50.88 -36.08 0.00
ID 0.73 0.68 0.00
Limit 44.37 23.84 0.02
Rating 21.65 23.91 0.13
Cards 1.50 -0.29 1.19
Age -0.22 1.67 0.00
Education -0.42 -0.50 -0.11
Balance -43.01 -24.18 -1.05
Gender_Female -0.81 -0.36 39.57
Student_Yes 12.76 7.38 -0.16
Married_Yes -0.33 0.03 0.99
Ethnicity_Asian -0.05 -0.45 0.00
Ethnicity_Caucasian 0.21 0.07 -0.09

Bazując o oszacowania parametrów w powyższej tabeli, parametry regresją grzbietową w 8/12 przypadkach dla zmiennych zależnych są niższe od parametrów regresji liniowej. Zatem możemy przypuszczać, że regresja istotnie może poprawić dokładność oszacowań dla badanego zjawiska.

Metoda lasso dla wielu zmiennych również dała oczekiwany efekt - współczynniki są bliżesz 0 niż w modelu liniowym. Jednak warto zauważyć, że parametr przy zmiennej Płeć jest bardzo wysoki, co jest dosyć niepokojące.

RMSE - prówananie

Regularyzacja Elastic Net - RMSE

Nie wszystko zadziałało mi, ze skryptu z prezentacji na enauczaniu, ale stworzyłam inny kod ( z pomocą chata GPT). Dzięki temu mogę porównać błędy RMSE między nieklasycznymi modelami.

Błąd RMSE dla metody elastic net jest równy 11.58164. Dodatkowo na podsatwie wykresu widzimy, że Im większa liczba regularyzacji tym błąd RMSE jest więjszy.

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard      11.7   Preprocessor1_Model1
## 2 rsq     standard       0.858 Preprocessor1_Model1
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        11.6

Model grzbietowy - RMSE

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard      10.8   Preprocessor1_Model1
## 2 rsq     standard       0.880 Preprocessor1_Model1
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        10.7

Błąd RMSE dla metody grzboetowej wynosi 10.73908.

W oparciu o wykres można zauwazyć, że błąd RMSE nie zmienia się znacznie jeśli zastosuję się między 0.01 do 0,1 regularyzacji. natomiast powyżej 0,1 wartść RMSE wzrasta wykładniczo.

Model lasso - RMSE

## # A tibble: 2 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 rmse    standard      10.8   Preprocessor1_Model1
## 2 rsq     standard       0.880 Preprocessor1_Model1
## # A tibble: 1 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard        10.7

Błąd RMSE dla metody lasso wynosi 10.73908.

## [1] 10.71915

RMSE dla modeli linowego wynosi ‘r train_rmse_lm’

Wnioski

Najniższy błąd RMSE uzyskano dla modeli metodą grzbietową, Lasso i regresji liniowej -> 10,73.

Model lasso i regresja nie wykazały znacznie niższych błędów (MSE, RMSE) niż się spodziewano. Jednak na podstawie powyższych analiz stwierdzono, że zastosowanie modeli zatówno grzebitowej jak i Lasso zminimalizuje szacowane parametry i poprawi jakość modelu. Co też potwierdza hipotezą początkową.

Wesołych Świąt i Szczęśliwego nowego Roku :)