Nieklasyczne metody statystyki
Regularyzacja
Małgorzata Dudanowicz
2024-12-22
To laboratorium na temat Regresji grzbietowej (Ridge Regression - RR)
i Lasso w R pochodzi ze stron 251-255 książki “Introduction to
Statistical Learning with Applications in R” autorstwa Garetha Jamesa,
Danieli Witten, Trevora Hastie i Roberta Tibshirani. Zostało ono
ponownie zaimplementowane jesienią 2016 roku w formacie
tidyverse przez Amelię McNamarę i R. Jordana Crousera w
Smith College.
W tym tygodniu omówimy dwie alternatywne formy regresji liniowej zwane regresją grzbietową i regresją LASSO. Te dwie metody są przykładami metod regularyzacji lub zmniejszania, w których zachęca się do tego, aby parametry modelu były małe.
Regresja Grzbietowa i Lasso
Wykorzystamy pakiet glmnet w celu przeprowadzenia
regresji ridge i lasso. Główną funkcją w tym pakiecie jest
glmnet(), która może być użyta do dopasowania modeli
regresji grzbietowej, modeli lasso i innych.
Funkcja ta ma nieco inną składnię niż inne funkcje dopasowujące modele, z którymi zetknęliśmy się do tej pory. W szczególności, musimy przekazać macierz \(x\) jak również wektor \(y\) i nie używamy składni \(y \sim x\).
Zanim przejdziemy dalej, upewnijmy się najpierw, że brakujące wartości zostały zostały usunięte z danych, jak opisano w poprzednim laboratorium.
Hitters = na.omit(Hitters)W raporcie tym przeprowadzimy regresję grzbietową i lasso, aby
przewidzieć Salary na danych Hitters.
Skonfigurujmy nasze dane:
x = model.matrix(Salary~., Hitters)[,-1] # przycinam pierwszą kolumnę
# zostawiam predyktory
y = Hitters %>%
select(Salary) %>%
unlist() %>%
as.numeric()Funkcja model.matrix() jest szczególnie przydatna do
tworzenia \(x\); nie tylko nie tylko
tworzy macierz odpowiadającą 19 predyktorom, ale również automatycznie
przekształca wszelkie zmienne jakościowe w zmienne dummy.
Ta ostatnia właściwość jest ważna, ponieważ glmnet()
może przyjmować tylko numeryczne, ilościowe dane wejściowe.
Bias vs Variance
Wybór modelu w problemach uczenia nadzorowanego wiąże się z realizacją dwóch sprzecznych celów:
1.) Model powienien być dobrze dopasowany do danych uczących, aby uchwycić zależność pomiędzy danymi.
2.) Model powinien dobrze przybliżać nieznane dane (zapewniać mały błąd generalizacji).
Modele złożone dobrze dopasowują się do danych wyjściowych, ale charakteryzują się dużą zmiennością wartości wyjściowych. Ryzykiem jest nadmierne dopasowanie = overfitting!
Modele prostsze są obciążone dużym błędem systematyczny (bias) i ich zastosowanie niesie ryzyko niewystarczającego dopasowania (underfitting)!
Składnikiem błędów generalizacji jest nieredukowalny błąd związany ze zmiennością danych.
Regularyzacja
Duża liczna zmiennych objaśniających (predyktorów): Metoda OLS nie daje jednoznacznego rozwiązania, gdy macierz XTX nie jest odwracalna (tzn. gdy zmienne objaśniające są liniowo zależne).
Taka sytuacja może mieć miejsce gdy zmiennych objaśniających jest tyle samo lub więcej niż obserwacji.
Duża wartość θi oznacza dużą wrażliwość funkcji regresji na drobne fluktuacje cechy!
Lepszym rozwiązaniem jest gorsze dopasowanie do danych uczących przy równoczesnym ograniczeniu parametrów świadczących o potencjalnie dużym błędzie generalizacji.
Regresja Grzbietowa
Wprowadzenie
Regresja grzbietowa (ang. Ridge regression) to technika regresji liniowej, która wprowadza regularyzację \(L_2\) do estymacji współczynników modelu. Regularyzacja \(L_2\) polega na dodaniu do funkcji celu kary proporcjonalnej do kwadratu wartości współczynników regresji.
Podstawową ideą regresji grzbietowej jest minimalizacja funkcji celu, która składa się z dwóch składników: błędu dopasowania (sumy kwadratów różnic pomiędzy rzeczywistymi wartościami odpowiedzi a przewidywanymi wartościami modelu) i kary regularyzacyjnej \(L_2\).
Wzór funkcji celu dla regresji grzbietowej można przedstawić jako: Minimize: RSS + \(\lambda \|\beta\|_2^2\), gdzie:
RSS to suma kwadratów różnic pomiędzy rzeczywistymi wartościami odpowiedzi a przewidywanymi wartościami modelu (błąd dopasowania),
\(\lambda\) (lambda) to parametr regularyzacji, który kontroluje siłę regularyzacji,
\(\|\beta\|_2^2\) to norma \(L_2\) współczynników regresji podniesiona do kwadratu.
Dodanie kary regularyzacyjnej \(L_2\) powoduje, że współczynniki regresji są skupione wokół zera, ale nie dokładnie równe zeru (chyba że \(\lambda\)=0).
Regresja grzbietowa zmniejsza wartości współczynników, ale nie powoduje, że stają się one równe zero. Im większa wartość \(\lambda\), tym bardziej są “sciskane” współczynniki regresji.
Regresja grzbietowa jest szczególnie przydatna, gdy mamy do czynienia z modelem, w którym występuje nadmierna wielowymiarowość lub wysokie korelacje między zmiennymi niezależnymi.
Poprzez zmniejszanie wartości współczynników, regresja grzbietowa może pomóc w redukcji wpływu mało istotnych cech, poprawić stabilność modelu i zmniejszyć ryzyko przeuczenia (overfitting).
Jednym ze sposobów kontroli złożoności modelu jest penalizacja jego wielkości. Na przykład, w problemie regresji liniowej:
\[ \min_{\beta \in \mathbb{R}^p} \sum_{i=1}^n (y_i - x_i^\top \beta)^2, \]
możemy kontrolować wielkość współczynników \(\beta\). Oczywiście wielkość \(\beta\) można zdefiniować na różne sposoby, np. norma-2: \(\|\beta\|_2\), norma-1: \(\|\beta\|_1\) czy norma-nieskończoność: \(\|\beta\|_{\infty}\). Regresja grzbietowa wiąże się z karą dwóch norm:
\[ \min_{\beta \in \mathbb{R}^p} \sum_{i=1}^n (y_i - x_i^\top \beta)^2 + \lambda \|\beta\|_2^2 \]
gdzie \(\lambda\) jest parametrem kontrolującym poziom regularyzacji. Zauważ, że \(X\) to macierz \(n\) na \(p\) wymiarów z wierszami: \(x_i^\top\), oraz \(Y\) to \(n\) na 1 wektor \(y_i\). Załóżmy, że \(X^\top X + \lambda I\) jest odwracalna, mamy dokładne rozwiązanie problemu regresji grzbietowej:
\[ \hat \beta_{ridge} = (X^\top X + \lambda I)^{-1}X^\top Y. \]
Przypomnijmy, że rozwiązaniem zwykłej regresji najmniejszych kwadratów jest (zakładając odwracalność macierzy \(X^\top X\)):
\[ \hat \beta_{ols} = (X^\top X)^{-1}X^\top Y. \]
Dwa fakty: kiedy \(\lambda \to 0\), \(\hat \beta_{ridge} \to \hat \beta_{ols}\); kiedy \(\lambda \to \infty\), \(\hat \beta_{ridge} \to 0\).
W szczególnych przypadkach \(X\) jest ortogonalna (tzn. kolumny \(X\) są ortogonalne), mamy:
\[ \hat \beta_{ridge} = \frac{\hat \beta_{ols}}{1 + \lambda}. \]
Widzimy więc, że estymator grzbietowy ma dodatkowo \(1/(1 + \lambda)\) tzw. “shrinkage factor”. W związku z tym na estymatorze grzbietowym występuje obciążliwość (bias).
Przykład
Funkcja glmnet() posiada argument alfa, który określa,
jaki typ modelu jest dopasowywany.
Jeśli alfa = 0 to dopasowywany jest model regresji
grzbietowej, a jeśli alfa = 1 to dopasowywany jest model
lasso.
Najpierw dopasowujemy model regresji grzbietowej:
grid = 10^seq(10, -2, length = 100)
ridge_mod = glmnet(x, y, alpha = 0, lambda = grid)Domyślnie funkcja glmnet() wykonuje regresję grzbietową
dla automatycznie wybranego wybranego zakresu wartości \(\lambda\). Jednakże, tutaj wybraliśmy
implementację funkcję w zakresie wartości od \(\lambda = 10^{10}\) do \(\lambda = 10^{-2}\), zasadniczo pokrywając
pełen zakres scenariuszy od modelu zerowego zawierającego tylko
przechwyt, do dopasowania najmniejszego kwadratu.
Jak widać, możemy również obliczyć dopasowanie modelu dla konkretnej wartości \(\lambda\), która nie jest jedną z oryginalnych wartości siatki.
Zauważ, że domyślnie funkcja glmnet() standaryzuje
zmienne tak, by były w tej samej skali. Aby wyłączyć to domyślne
ustawienie, użyj argumentu standardize = FALSE.
Z każdą wartością \(\lambda\)
związany jest wektor współczynników regresji grzbietowej, przechowywany
w macierzy, do której można uzyskać dostęp przez coef(). W
tym przypadku jest to macierz \(20 \times
100\), z 20 wierszami (po jednym dla każdego predyktora, plus
intercept) i 100 kolumnami (po jednej dla każdej wartości \(\lambda\)).
dim(coef(ridge_mod))## [1] 20 100plot(ridge_mod) # wykres współczynników
Spodziewamy się, że oszacowania współczynników będą znacznie mniejsze, w sensie normy \(l_2\), gdy używana jest duża wartość \(\lambda\), w porównaniu z małą wartością \(\lambda\).
Oto współczynniki, gdy \(\lambda = 11498\), wraz z ich normą \(l_2\):
ridge_mod$lambda[50] # Wyświetl 50-tą wartość lambdy## [1] 11497.57coef(ridge_mod)[,50] # Wyświetl współczynniki związane z 50-tą wartością lambdy## (Intercept) AtBat Hits HmRun Runs
## 407.356050200 0.036957182 0.138180344 0.524629976 0.230701523
## RBI Walks Years CAtBat CHits
## 0.239841459 0.289618741 1.107702929 0.003131815 0.011653637
## CHmRun CRuns CRBI CWalks LeagueN
## 0.087545670 0.023379882 0.024138320 0.025015421 0.085028114
## DivisionW PutOuts Assists Errors NewLeagueN
## -6.215440973 0.016482577 0.002612988 -0.020502690 0.301433531sqrt(sum(coef(ridge_mod)[-1,50]^2)) # Oblicz normę l2## [1] 6.360612Dla kontrastu, oto współczynniki, gdy \(\lambda = 705\), wraz z ich \(l_2\) normą. Zwróć uwagę na znacznie większą normę \(l_2\) współczynników związanych z tą mniejszą wartością \(\lambda\).
ridge_mod$lambda[60] # Wyświetl 60-tą wartość lambdy## [1] 705.4802coef(ridge_mod)[,60] # Wyświetl współczynniki powiązane z 60-tą wartość lambdy## (Intercept) AtBat Hits HmRun Runs RBI
## 54.32519950 0.11211115 0.65622409 1.17980910 0.93769713 0.84718546
## Walks Years CAtBat CHits CHmRun CRuns
## 1.31987948 2.59640425 0.01083413 0.04674557 0.33777318 0.09355528
## CRBI CWalks LeagueN DivisionW PutOuts Assists
## 0.09780402 0.07189612 13.68370191 -54.65877750 0.11852289 0.01606037
## Errors NewLeagueN
## -0.70358655 8.61181213sqrt(sum(coef(ridge_mod)[-1,60]^2)) # Oblicz normę l2## [1] 57.11001Funkcję predict() możemy wykorzystać do wielu celów. Na
przykład, możemy uzyskać współczynniki regresji grzbietowej dla nowej
wartości \(\lambda\), powiedzmy 50:
predict(ridge_mod, s = 50, type = "coefficients")[1:20,]## (Intercept) AtBat Hits HmRun Runs
## 4.876610e+01 -3.580999e-01 1.969359e+00 -1.278248e+00 1.145892e+00
## RBI Walks Years CAtBat CHits
## 8.038292e-01 2.716186e+00 -6.218319e+00 5.447837e-03 1.064895e-01
## CHmRun CRuns CRBI CWalks LeagueN
## 6.244860e-01 2.214985e-01 2.186914e-01 -1.500245e-01 4.592589e+01
## DivisionW PutOuts Assists Errors NewLeagueN
## -1.182011e+02 2.502322e-01 1.215665e-01 -3.278600e+00 -9.496680e+00Podzielimy teraz próbki na zbiór treningowy i testowy w celu oszacować błąd testu regresji grzbietowej i lasso.
set.seed(1)
train = Hitters %>%
sample_frac(0.5)
test = Hitters %>%
setdiff(train)
x_train = model.matrix(Salary~., train)[,-1]
x_test = model.matrix(Salary~., test)[,-1]
y_train = train %>%
select(Salary) %>%
unlist() %>%
as.numeric()
y_test = test %>%
select(Salary) %>%
unlist() %>%
as.numeric()Następnie dopasowujemy model regresji grzbietowej na zbiorze
treningowym i oceniamy jego MSE na zbiorze testowym, używając \(\lambda = 4\). Zwróć uwagę na użycie
funkcji predict(). Ponownie: tym razem otrzymujemy
przewidywania dla zbioru testowego, zastępując
type="coefficients" argumentem newx.
ridge_mod = glmnet(x_train, y_train, alpha=0, lambda = grid, thresh = 1e-12)
ridge_pred = predict(ridge_mod, s = 4, newx = x_test)
mean((ridge_pred - y_test)^2)## [1] 139858.6Testowe MSE wynosi 139858. Zauważ, że gdybyśmy zamiast tego dopasowali po prostu model tylko z wyrazem wolnym, przewidywalibyśmy każdą obserwację testową używając średniej z obserwacji zbioru treningowego. W takim przypadku moglibyśmy obliczyć MSE zestawu testowego w ten sposób:
mean((mean(y_train) - y_test)^2)## [1] 224692.1Moglibyśmy również uzyskać ten sam wynik, dopasowując model regresji
grzbietowej z bardzo dużą wartością \(\lambda\). Zauważ, że 1e10
oznacza \(10^{10}\).
ridge_pred = predict(ridge_mod, s = 1e10, newx = x_test)
mean((ridge_pred - y_test)^2)## [1] 224692.1Tak więc dopasowanie modelu regresji grzbietowej z \(\lambda = 4\) prowadzi do znacznie niższego testu MSE niż dopasowanie modelu z samym przechwytem.
Sprawdzimy teraz, czy jest jakaś korzyść z wykonania regresji grzbietowej z \(\lambda = 4\) zamiast po prostu wykonać regresję najmniejszych kwadratów.
Przypomnijmy, że najmniejsza kwadratura to po prostu regresja grzbietowa z \(\lambda = 0\).
* Uwaga: Aby glmnet() dawał dokładne
(exact) współczynniki najmniejszego kwadratu, gdy \(\lambda = 0\), używamy argumentu
exact=T przy wywołaniu funkcji predict(). W
przeciwnym razie, funkcja predict() będzie interpolować nad
siatką wartości \(\lambda\) użytą w
dopasowaniu modelu glmnet(), dając przybliżone wyniki.
Nawet gdy użyjemy exact = T, pozostaje niewielka
rozbieżność na trzecim miejscu po przecinku między wynikami
glmnet(), gdy \(\lambda =
0\) i wyjściem z lm(); jest to spowodowane
numerycznym przybliżeniem ze strony glmnet().
ridge_pred = predict(ridge_mod, s = 0, newx = x_test)
mean((ridge_pred - y_test)^2)## [1] 174060lm(Salary~., data = train)##
## Call:
## lm(formula = Salary ~ ., data = train)
##
## Coefficients:
## (Intercept) AtBat Hits HmRun Runs RBI
## 2.398e+02 -1.639e-03 -2.179e+00 6.337e+00 7.139e-01 8.735e-01
## Walks Years CAtBat CHits CHmRun CRuns
## 3.594e+00 -1.309e+01 -7.136e-01 3.316e+00 3.407e+00 -5.671e-01
## CRBI CWalks LeagueN DivisionW PutOuts Assists
## -7.525e-01 2.347e-01 1.322e+02 -1.346e+02 2.099e-01 6.229e-01
## Errors NewLeagueN
## -4.616e+00 -8.330e+01predict(ridge_mod, s = 0, type="coefficients")[1:20,]## (Intercept) AtBat Hits HmRun Runs
## 239.89368111 -0.01946204 -2.07305757 6.44254692 0.64610179
## RBI Walks Years CAtBat CHits
## 0.82179888 3.62448842 -13.28142313 -0.70314292 3.26064805
## CHmRun CRuns CRBI CWalks LeagueN
## 3.33170237 -0.54000590 -0.72015101 0.22582579 131.41324242
## DivisionW PutOuts Assists Errors NewLeagueN
## -134.76073238 0.20949301 0.61942855 -4.58545824 -82.35090554Wygląda na to, że rzeczywiście poprawiamy się w stosunku do zwykłego najmniejszego kwadratu!
Uwaga: ogólnie, jeśli chcemy dopasować (niespenalizowany) model
najmniejszych kwadratów, to powinniśmy użyć funkcji lm(),
ponieważ ta funkcja dostarcza bardziej użytecznych wyjścia, takie jak
błędy standardowe i wartości \(p\) dla
współczynników.
Zamiast arbitralnie wybierać \(\lambda =
4\), lepiej byłoby użyć walidacji krzyżowej do wyboru parametru
dostrojenia \(\lambda\). Możemy to
zrobić używając wbudowanej funkcji walidacji krzyżowej,
cv.glmnet(). Domyślnie funkcja ta wykonuje 10-krotną
walidację krzyżową, choć można to zmienić używając argumentu argumentu
folds. Zauważ, że najpierw ustawiamy losowe ziarno, aby
nasze wyniki były powtarzalne, ponieważ wybór krotności walidacji
krzyżowej jest losowy.
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] 326.1406Widzimy zatem, że wartość \(\lambda\), która powoduje najmniejszy błąd walidacji krzyżowej to 326. Możemy również wykreślić MSE jako funkcję \(\lambda\):
plot(cv.out) # Narysuj wykres treningowego MSE jako funkcję lambda
Jaki jest testowy MSE związany z tą wartością \(\lambda\)?
ridge_pred = predict(ridge_mod, s = bestlam, newx = x_test) # Użyj najlepszej lambdy do przewidywania danych testowych
mean((ridge_pred - y_test)^2) # Oblicz testowe MSE## [1] 140056.2Stanowi to dalszą poprawę w stosunku do testowego MSE, które uzyskaliśmy używając \(\lambda = 4\). Ostatecznie, ponownie wyznaczamy nasz model regresji grzbietowej na pełnym zestawie danych, używając wartości \(\lambda\) wybranej w walidacji krzyżowej, i sprawdzamy oszacowania współczynników.
out = glmnet(x, y, alpha = 0) # Dopasuj model regresji grzbietowej do pełnego zbioru danych
predict(out, type = "coefficients", s = bestlam)[1:20,] # Wyświetlanie współczynników przy użyciu lambda wybranego przez CV## (Intercept) AtBat Hits HmRun Runs RBI
## 15.44834992 0.07716945 0.85906253 0.60120338 1.06366687 0.87936073
## Walks Years CAtBat CHits CHmRun CRuns
## 1.62437580 1.35296285 0.01134998 0.05746377 0.40678422 0.11455696
## CRBI CWalks LeagueN DivisionW PutOuts Assists
## 0.12115916 0.05299953 22.08942756 -79.03490992 0.16618830 0.02941513
## Errors NewLeagueN
## -1.36075645 9.12528397Zgodnie z oczekiwaniami, żaden ze współczynników nie jest dokładnie zerowy - regresja grzbietowa nie dokonuje selekcji zmiennych!
Regresja Lasso
Wprowadzenie
Zamiast regularyzacji \(L_2\), LASSO używa penalizacji \(L_1\), to znaczy:
\[ \min_{\beta \in \mathbb{R}^p} \sum_{i=1}^n (y_i - x_i^\top \beta)^2 + \lambda \|\beta\|_1. \]
Ze względu na charakter normy \(L_1\), LASSO ma tendencję do dawania bardziej rzadkich rozwiązań niż regresja grzbietowa. Jest to typowo użyteczne w ustawieniach wielowymiarowych, gdy prawdziwy model jest w rzeczywistości niskowymiarowym osadzeniem.
Model regresji lasso został pierwotnie opracowany w 1989 roku. Jest to alternatywa dla klasycznego oszacowania metodą najmniejszych kwadratów, która unika wielu problemów z nadmiernym dopasowaniem (overfittingiem), gdy mamy dużą liczbę niezależnych zmiennych.
Regresja Lasso (Least Absolute Shrinkage and Selection Operator) to technika regresji liniowej stosowana do oszacowania współczynników modelu, która wprowadza regularyzację \(L_1\). Regularyzacja L1 polega na dodaniu do funkcji celu kary proporcjonalnej do wartości bezwzględnej współczynników regresji.
Regresja Lasso ma zdolność do jednoczesnego wykonania selekcji cech i regularyzacji, co oznacza, że może pomóc w identyfikacji najbardziej istotnych cech modelu, a także zmniejszyć wpływ mniej istotnych cech.
Podstawowym celem regresji Lasso jest minimalizacja funkcji celu, która składa się z dwóch składników: błędu dopasowania (sumy kwadratów różnic pomiędzy rzeczywistymi wartościami odpowiedzi a przewidywanymi wartościami modelu) i kary regularyzacyjnej \(L_1\).
Wzór funkcji celu dla regresji Lasso może być przedstawiony jako: Minimize: RSS + \(\lambda \|\beta\|_1\), gdzie:
RSS to suma kwadratów różnic pomiędzy rzeczywistymi wartościami odpowiedzi a przewidywanymi wartościami modelu (błąd dopasowania),
\(\lambda\) (lambda) to parametr regularyzacji, który kontroluje siłę regularyzacji, a \(\|\beta\|_1\) to norma \(L_1\) współczynników regresji.
Dodanie kary regularyzacyjnej \(L_1\) powoduje, że niektóre współczynniki regresji stają się równe zero, co prowadzi do selekcji cech. Im większa wartość \(\lambda\), tym większa jest tendencja do redukcji współczynników do zera, prowadząc do bardziej rzadkiego modelu z mniejszą liczbą cech.
Regresja Lasso jest przydatna w przypadkach, gdy mamy do czynienia z wieloma cechami, z których niektóre mogą być nieistotne. Może pomóc w identyfikacji istotnych cech, redukcji nadmiaru danych i zwiększeniu interpretowalności modelu.
Przykład
Zobaczyliśmy, że regresja grzbietowa z mądrym wyborem \(\lambda\) może przewyższać metodę najmniejszych kwadratów, jak również model zerowy na zbiorze danych Hitters.
Teraz zobaczmy, czy lasso może dać albo dokładniejszy, albo bardziej interpretowalny model niż regresja grzbietowa.
W celu dopasowania modelu lasso, po raz kolejny używamy funkcji
glmnet(), jednak tym razem używamy argumentu
alpha=1. Poza tą zmianą postępujemy tak samo jak w
przypadku dopasowywania modelu regresji grzbietowej:
lasso_mod = glmnet(x_train,
y_train,
alpha = 1,
lambda = grid) # Dopasuj model lasso do danych treningowych
plot(lasso_mod) # Wykreśl współczynniki
Zauważmy, że na wykresie współczynników, w zależności od wyboru dostrojenia parametru, niektóre ze współczynników są dokładnie równe zeru. Teraz przeprowadzimy walidację krzyżową i obliczymy związany z nią błąd testu:
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
mean((lasso_pred - y_test)^2) # Oblicz MSE w próbie testowej## [1] 143273Jest to znacznie niższe MSE zbioru testowego niż modelu zerowego i modelu najmniejszych kwadratów, i bardzo podobny do MSE testu regresji grzbietowej z \(\lambda\) wybranej przez walidację krzyżową.
Jednakże lasso ma istotną przewagę nad regresją grzbietową w tym, że wynikowe oszacowania współczynników są rzadkie. Tutaj widzimy, że 12 z 19 oszacowań współczynników jest dokładnie zerowych:
out = glmnet(x, y, alpha = 1, lambda = grid) # Dopasuj model lasso do pełnego zbioru danych
lasso_coef = predict(out, type = "coefficients", s = bestlam)[1:20,] # Wyświetlanie współczynników przy użyciu lambda wybranego przez CV
lasso_coef## (Intercept) AtBat Hits HmRun Runs
## 1.27429897 -0.05490834 2.18012455 0.00000000 0.00000000
## RBI Walks Years CAtBat CHits
## 0.00000000 2.29189433 -0.33767315 0.00000000 0.00000000
## CHmRun CRuns CRBI CWalks LeagueN
## 0.02822467 0.21627609 0.41713051 0.00000000 20.28190194
## DivisionW PutOuts Assists Errors NewLeagueN
## -116.16524424 0.23751978 0.00000000 -0.85604181 0.00000000Wybierając tylko predyktory o niezerowych współczynnikach widzimy, że model lasso z \(\lambda\) wybranym przez walidację krzyżową zawiera tylko siedem zmiennych:
lasso_coef[lasso_coef != 0] # Wyświetlanie tylko niezerowych współczynników## (Intercept) AtBat Hits Walks Years
## 1.27429897 -0.05490834 2.18012455 2.29189433 -0.33767315
## CHmRun CRuns CRBI LeagueN DivisionW
## 0.02822467 0.21627609 0.41713051 20.28190194 -116.16524424
## PutOuts Errors
## 0.23751978 -0.85604181Twoja 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!
library(ggplot2)
data(diamonds)
set.seed(123)
head(diamonds)## # A tibble: 6 × 10
## carat cut color clarity depth table price x y z
## <dbl> <ord> <ord> <ord> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.23 Ideal E SI2 61.5 55 326 3.95 3.98 2.43
## 2 0.21 Premium E SI1 59.8 61 326 3.89 3.84 2.31
## 3 0.23 Good E VS1 56.9 65 327 4.05 4.07 2.31
## 4 0.29 Premium I VS2 62.4 58 334 4.2 4.23 2.63
## 5 0.31 Good J SI2 63.3 58 335 4.34 4.35 2.75
## 6 0.24 Very Good J VVS2 62.8 57 336 3.94 3.96 2.48Zbiór danych zawiera informacje o 53 940 diamentach, ich cechach fizycznych, jakości oraz cenach. Zawiera takie zmienne jak:
carat - waga diamentu w karatach cut - jakość szlifu diamentu (Fair, Good, Very Good, Premium, Ideal) color - kolor diamentu, oceniany w skali od D (najlepszy, bezbarwny) do J (najgorszy, lekko żółtawy) clarity - przejrzystość diamentu, opisująca liczbę i rodzaj inkluzji (niedoskonałości) depth - głębokość diamentu wyrażona jako procent średnicy table - szerokość górnej płaszczyzny diamentu (zwanej table) wyrażona jako procent szerokości price - cena diamentu w dolarach amerykańskich x - długość diamentu w milimetrach y - szerokość diamentu w milimetrach z - głębokość diamentu w milimetrach
Przegląd danych
Na początku zróbmy szybki przegląd naszego zbioru danych:
ggplot(diamonds, aes(x = price)) +
geom_histogram(bins = 50, color = "black", fill = "blue") +
labs(title = "Rozkład cen diamentów", x = "Cena ($)", y = "Liczba diamentów")
Z histogramu widać, że większość diamentów ma ceny skupione w niższych przedziałach (0-5000 USD), co wskazuje na prawostronnie skośny rozkład. Ceny droższych diamentów są znacznie mniej liczne.Najwięcej diamentów znajduje się w przedziale cenowym około 500-1500 USD. To sugeruje, że większość diamentów w zbiorze to produkty o stosunkowo niskiej wartości.
ggplot(diamonds, aes(x = carat, y = price)) +
geom_point(alpha = 0.5) +
labs(title = "Cena a masa diamentu", x = "Masa (carat)", y = "Cena ($)")
Ogólnie obserwujemy, że wraz ze wzrostem masy diamentu (carat) cena (price) rośnie, co wskazuje na silną dodatnią zależność między tymi zmiennymi. Większe diamenty są droższe, co jest zgodne z intuicją.
ggplot(diamonds, aes(x = cut, y = price)) +
geom_boxplot() +
labs(title = "Cena w zależności od jakości szlifu", x = "Jakość szlifu", y = "Cena ($)")
corr_data <- diamonds %>%
select_if(is.numeric) %>%
cor()
corrplot(corr_data, method = "ellipse", type = "lower", tl.cex = 0.8)
Najsilniejszą korelację z ceną ma masa (carat).
ggplot(diamonds, aes(x = clarity, fill = clarity)) +
geom_bar() +
labs(title = "Rozkład przejrzystości diamentów", x = "Przejrzystość", y = "Liczba diamentów")
Niektóre kategorie przejrzystości (SI1, VS2) są znacznie bardziej reprezentowane niż inne.
Wartości odstające
ggplot(diamonds, aes(y = price)) +
geom_boxplot(fill = "blue", color = "black", outlier.color = "red", outlier.size = 2) +
labs(title = "Rozkład cen diamentów (z wartościami odstającymi)", y = "Cena ($)", x = "") +
theme_minimal()
Q1 <- quantile(diamonds$price, 0.25)
Q3 <- quantile(diamonds$price, 0.75)
IQR <- Q3 - Q1
lower_bound <- Q1 - 1.5 * IQR
upper_bound <- Q3 + 1.5 * IQR
diamonds <- subset(diamonds, price >= lower_bound & price <= upper_bound)
ggplot(diamonds, aes(y = price)) +
geom_boxplot(fill = "green", color = "black") +
labs(title = "Rozkład cen diamentów (bez wartości odstających)", y = "Cena ($)", x = "") +
theme_minimal()
W analizie zdecydowano się na usunięcie wartości odstających z uwagi na ich potencjalnie negatywny wpływ na jakość modelu predykcyjnego oraz miary błędu.Wartości odstające mogą znacząco zwiększać wartość średniego błędu kwadratowego (MSE), co może sprawiać, że model stanie się mniej użyteczny w przewidywaniu dla typowych obserwacji.Mogą także negatywnie wpływać na współczynniki regresji, gdyż model jest zmuszony do nadmiernego dopasowania, aby uwzględnić te nietypowe dane.
Kodowanie zmiennych jakościowych i podział na zbiory
# Kodowanie zmiennych jakościowych
x <- model.matrix(price ~ carat + depth + color + clarity + cut + table + x + y + z, data = diamonds)[, -1]
# Zmienna objaśniana
y <- diamonds$price
x_scaled <- scale(x)
# Podział danych na zbiory treningowe i testowe
trainIndex <- createDataPartition(y, p = 0.7, list = FALSE)
x_train <- x[trainIndex, ]
x_test <- x[-trainIndex, ]
y_train <- y[trainIndex]
y_test <- y[-trainIndex]W zbiorze diamonds zmienne takie jak cut, color, i clarity są zmiennymi jakościowymi. Modele statystyczne, takie jak regresja grzbietowa i Lasso, wymagają danych numerycznych. Funkcja model.matrix automatycznie przekształca zmienne jakościowe na zmienne typu dummy (zero-jedynkowe), dzięki czemu można je uwzględnić w modelu.Funkcja createDataPartition zapewnia, że dane są podzielone w sposób losowy, ale uwzględnia proporcje w przypadku problemów klasyfikacyjnych (stratified sampling).
Regresja OLS
ols_model <- lm(price ~ carat + depth + + color + clarity + cut + table + x + y, data = diamonds[trainIndex, ])
summary(ols_model)##
## Call:
## lm(formula = price ~ carat + depth + +color + clarity + cut +
## table + x + y, data = diamonds[trainIndex, ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -11621.0 -450.4 -126.7 296.6 5469.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1756.742 343.853 5.109 3.26e-07 ***
## carat 8767.981 60.054 146.002 < 2e-16 ***
## depth -32.410 3.625 -8.940 < 2e-16 ***
## color.L -1429.951 15.241 -93.823 < 2e-16 ***
## color.Q -486.987 13.861 -35.134 < 2e-16 ***
## color.C -132.360 12.922 -10.243 < 2e-16 ***
## color^4 46.805 11.859 3.947 7.93e-05 ***
## color^5 -48.905 11.159 -4.382 1.18e-05 ***
## color^6 -29.567 10.088 -2.931 0.00338 **
## clarity.L 3191.578 25.838 123.522 < 2e-16 ***
## clarity.Q -1602.252 23.920 -66.983 < 2e-16 ***
## clarity.C 695.128 20.493 33.920 < 2e-16 ***
## clarity^4 -323.818 16.431 -19.707 < 2e-16 ***
## clarity^5 216.747 13.503 16.051 < 2e-16 ***
## clarity^6 37.983 11.793 3.221 0.00128 **
## clarity^7 82.010 10.447 7.850 4.26e-15 ***
## cut.L 486.616 19.391 25.095 < 2e-16 ***
## cut.Q -256.109 15.473 -16.552 < 2e-16 ***
## cut.C 116.556 13.387 8.707 < 2e-16 ***
## cut^4 -23.609 10.714 -2.204 0.02755 *
## table -18.508 2.531 -7.313 2.66e-13 ***
## x -478.285 34.582 -13.830 < 2e-16 ***
## y 58.850 27.327 2.154 0.03128 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 792 on 35259 degrees of freedom
## Multiple R-squared: 0.9182, Adjusted R-squared: 0.9182
## F-statistic: 1.799e+04 on 22 and 35259 DF, p-value: < 2.2e-16Większość zmiennych jest statystycznie istotna, co oznacza, że mają realny wpływ na cenę diamentów. Masa (carat) ma najsilniejszy wpływ na cenę diamentu. Każdy dodatkowy karat zwiększa cenę średnio o 11236 USD. Cechy jakościowe, takie jak kolor (color) i przejrzystość (clarity), również znacząco wpływają na cenę.Bardzo wysoka wartość R² (91.97%) wskazuje, że model dobrze wyjaśnia zmienność cen.
# Predykcje na zbiorze testowym
ols_pred <- predict(ols_model, newdata = diamonds[-trainIndex, ])
ols_mse <- mean((diamonds$price[-trainIndex] - ols_pred)^2)
# Obliczenie R²
ols_r2 <- 1 - (sum((diamonds$price[-trainIndex] - ols_pred)^2) /
sum((diamonds$price[-trainIndex] - mean(diamonds$price[-trainIndex]))^2))
# Wyświetlenie wyników
cat("OLS MSE:", ols_mse, "\nOLS R^2:", ols_r2, "\n")## OLS MSE: 627997.5
## OLS R^2: 0.9173359Model OLS działa dobrze pod względem wyjaśniania zmienności danych (wysokie R²), ale wysoka wartość MSE wskazuje na problemy z dokładnością prognoz, szczególnie dla ekstremalnych wartości ceny.
Regresja ridge
X_train_matrix <- as.matrix(x_train)
X_test_matrix <- as.matrix(x_test)
ridge_model <- cv.glmnet(x_train, y_train, alpha = 0)
plot(ridge_model)
# Optymalna lambda
best_lambda_ridge <- ridge_model$lambda.min
best_lambda_ridge## [1] 253.4909Model użył 10-krotnej walidacji krzyżowej do oszacowania błędu średniokwadratowego (MSE) dla różnych wartości parametru lambda.Na wykresie widać, że dla bardzo małych wartości lambda model jest mniej regularyzowany i osiąga niższy błąd MSE, ale może być bardziej podatny na przeuczenie.Dla bardzo dużych wartości lambda model jest silnie regularyzowany, co prowadzi do uproszczenia modelu, ale także do zwiększenia błędu (niedopasowania).Wybrana wartość lambda.min minimalizuje błąd walidacji krzyżowej, zapewniając najlepszy kompromis między dokładnością modelu a jego prostotą.
# Predykcje na zbiorze testowym
ridge_pred <- predict(ridge_model, s = best_lambda_ridge, newx = x_test)
# Obliczenie MSE
ridge_mse <- mean((y_test - ridge_pred)^2)
# Obliczenie R²
ridge_r2 <- 1 - (sum((y_test - ridge_pred)^2) / sum((y_test - mean(y_test))^2))
# Wyświetlenie wyników
cat("Ridge MSE:", ridge_mse, "\nRidge R^2:", ridge_r2, "\n")## Ridge MSE: 777390.9
## Ridge R^2: 0.897671Regresja Lasso
lasso_model <- cv.glmnet(x_train, y_train, alpha = 1)
plot(lasso_model)
Regresja Lasso jest metodą regularyzacyjną, która, w przeciwieństwie do regresji grzbietowej, może całkowicie wyeliminować nieistotne predyktory poprzez “zerowanie” ich współczynników.
best_lambda_lasso <- lasso_model$lambda.min
best_lambda_lasso## [1] 1.232625lasso_pred <- predict(lasso_model, s = best_lambda_lasso, newx = x_test)
# Obliczenie MSE dla Lasso
lasso_mse <- mean((y_test - lasso_pred)^2)
# Obliczenie R² dla Lasso
lasso_r2 <- 1 - (sum((y_test - lasso_pred)^2) / sum((y_test - mean(y_test))^2))
# Wyświetlenie wyników
cat("Lasso MSE:", lasso_mse, "\nLasso R^2:", lasso_r2, "\n")## Lasso MSE: 628548.7
## Lasso R^2: 0.9172633Odpowiedzi na pytania:
- Który zbiór danych wybrałeś?
Wybrany zbiór to diamonds z pakietu ggplot2, zawierający dane o 53 940 diamentach, ich cechach fizycznych i cenach.
- Jaka była Twoja zmienna zależna (tzn. co próbowałeś modelować)?
Zmienną zależną jest price (cena diamentów w dolarach amerykańskich).
- Czy oczekiwałeś, że regresja grzbietowa będzie lepsza od lasso, czy odwrotnie? Jak wypada w stosunku do OLS? Pokaż odpowiednie raporty, miary dopasowania i krótko je omów (porównaj).
Regresja OLS: Model OLS ma wysokie R² (91.97%), co oznacza, że bardzo dobrze wyjaśnia zmienność ceny. Najsilniejszym predyktorem jest masa diamentu (carat). Wysoka wartość MSE sugeruje, że OLS może mieć trudności z prognozowaniem cen dla rzadkich, ekstremalnych wartości.
Regresja Ridge: Ridge radzi sobie dobrze przy minimalizowaniu problemów współliniowości, ale w tym przypadku, gdzie dane są dobrze dopasowane przez OLS, Ridge wypadło gorzej. MSE Ridge jest wyższe, ponieważ regularyzacja “zmiękcza” wpływ najistotniejszych predyktorów, takich jak carat.
Regresja Lasso: Lasso dodatkowo eliminuje nieistotne predyktory poprzez zerowanie ich współczynników. Jednak w Twoich danych większość predyktorów wydaje się istotna, co mogło wpłynąć na relatywnie gorsze wyniki Lasso w porównaniu do OLS.
Najważniejsze wnioski:
results <- data.frame(
Model = c("OLS", "Ridge", "Lasso"),
MSE = c(ols_mse, ridge_mse, lasso_mse),
R2 = c(ols_r2, ridge_r2, lasso_r2)
)
kable(results, caption = "Podsumowanie wyników modeli") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = F, position = "center") %>%
column_spec(2, width = "10em") %>%
column_spec(3, width = "10em")| Model | MSE | R2 |
|---|---|---|
| OLS | 627997.5 | 0.9173359 |
| Ridge | 777390.9 | 0.8976710 |
| Lasso | 628548.7 | 0.9172633 |
OLS wypadło najlepiej, ponieważ:
W zbiorze danych diamonds zmienna carat (masa diamentu w karatach) ma wyjątkowo silną, dodatnią korelację z ceną (price). Korelacja ta wynosi powyżej 0.9, co oznacza, że carat wyjaśnia większość zmienności ceny. W przypadku OLS model nie jest ograniczany przez regularyzację, więc zmienne takie jak carat mogą w pełni pokazać swój wpływ. W Ridge i Lasso regularyzacja ogranicza współczynniki regresji, w tym także współczynnika dla carat. W rezultacie modele te “osłabiają” wpływ najważniejszego predyktora, co prowadzi do zwiększenia błędu (MSE).
Ridge i Lasso wykazują swoje zalety w danych, gdzie predyktory są silnie współliniowe (tj. istnieją silne korelacje między zmiennymi objaśniającymi).W przypadku danych diamonds współliniowość między predyktorami, takimi jak carat, depth, table, i x, jest umiarkowana i nie powoduje istotnych problemów dla OLS.
W zbiorze diamonds wszystkie predyktory (takie jak carat, clarity, color, cut, x, y, z) mają intuicyjny wpływ na cenę diamentu. Nie ma “zbędnych” predyktorów, które należy usunąć lub wyzerować.
W Ridge i Lasso regularyzacja wprowadziła jedynie niepotrzebne “osłabienie” wpływu predyktorów, które są rzeczywiście ważne.
Dane diamonds są stosunkowo “czyste” i dobrze zorganizowane. Po usunięciu wartości odstających nie zawierają dużego poziomu szumu.