Nieklasyczne metody statystyki

Regularyzacja

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 dobdrze 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 100

plot(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.57

coef(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.301433531

sqrt(sum(coef(ridge_mod)[-1,50]^2)) # Oblicz normę l2

## [1] 6.360612

Dla 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.4802

coef(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.61181213

sqrt(sum(coef(ridge_mod)[-1,60]^2)) # Oblicz normę l2

## [1] 57.11001

Funkcję 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+00

Podzielimy 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.6

Testowe 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.1

Moglibyś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.1

Tak 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] 174060

lm(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+01

predict(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.35090554

Wyglą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.1406

Widzimy 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.2

Stanowi 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.12528397

Zgodnie 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] 143273

Jest 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.00000000

Wybierają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.85604181

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! Aby zaliczyć to laboratorium, zamieść odpowiedzi na następujące pytania:

Zbiór danych

• Który zbiór danych wybrałeś?

Wybrany datase: Statystyki dla dużej liczby amerykańskich szkół wyższych z wydania US News and World Report z 1995 roku.

# Wybrany datase:
# Statystyki dla dużej liczby amerykańskich szkół wyższych z wydania US News and World Report z 1995 roku.

data("College")
head(College)

##                              Private Apps Accept Enroll Top10perc Top25perc
## Abilene Christian University     Yes 1660   1232    721        23        52
## Adelphi University               Yes 2186   1924    512        16        29
## Adrian College                   Yes 1428   1097    336        22        50
## Agnes Scott College              Yes  417    349    137        60        89
## Alaska Pacific University        Yes  193    146     55        16        44
## Albertson College                Yes  587    479    158        38        62
##                              F.Undergrad P.Undergrad Outstate Room.Board Books
## Abilene Christian University        2885         537     7440       3300   450
## Adelphi University                  2683        1227    12280       6450   750
## Adrian College                      1036          99    11250       3750   400
## Agnes Scott College                  510          63    12960       5450   450
## Alaska Pacific University            249         869     7560       4120   800
## Albertson College                    678          41    13500       3335   500
##                              Personal PhD Terminal S.F.Ratio perc.alumni Expend
## Abilene Christian University     2200  70       78      18.1          12   7041
## Adelphi University               1500  29       30      12.2          16  10527
## Adrian College                   1165  53       66      12.9          30   8735
## Agnes Scott College               875  92       97       7.7          37  19016
## Alaska Pacific University        1500  76       72      11.9           2  10922
## Albertson College                 675  67       73       9.4          11   9727
##                              Grad.Rate
## Abilene Christian University        60
## Adelphi University                  56
## Adrian College                      54
## Agnes Scott College                 59
## Alaska Pacific University           15
## Albertson College                   55

Zmienna zależna

• Jaka była Twoja zmienna zależna (tzn. co próbowałeś modelować)?

Wybrana zmienna objaśniana - Wskaźnik ukończenia studiów

y = College %>%
  select(Grad.Rate) %>%
  unlist() %>%
  as.numeric()

x = model.matrix(Grad.Rate~., College)[,-1]

Regresja grzbietowa, Lasso, OLS

• 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).

W pierwszej kolejnosci należy podzielić próbki na zbiór treningowy i testowy w celu oszacować błąd testu regresji grzbietowej i lasso oraz OLS.

# Zastosowany został stosunek próby treningowej do testowej 80/20
set.seed(1)

train = College %>%
  sample_frac(0.8) 

test = College %>%
  setdiff(train)

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

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

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

1. Regresja grzbietowa

W pierwszej kolejności przetestowana zostanie regresja grzbietowa (alfa=0)

grid = 10^seq(10, -2, length = 100)
ridge_mod = glmnet(x, y, alpha = 0, lambda = grid)
dim(coef(ridge_mod))

## [1]  18 100

plot(ridge_mod)

Następnie identyfikujemy wartość lambdy, która najbardziej minimalizuje błąd MSE

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

## [1] 3.349577

plot(cv.out)

W kolejnym kroku wyznaczamy prognozę na bazie danych treningowych. Oszacowany został również błąd MSE oraz wartości współczynników w modelu.

ridge_pred = predict(ridge_mod, s = bestlam, newx = x_test) # Zastosowanie najlepiszej wartości lambdy
MSE_reg_g = mean((ridge_pred - y_test)^2) # MSE
print(paste("Wartość błędu MSE dla prognozy za pomocą regresji grzbietowej:",  MSE_reg_g))

## [1] "Wartość błędu MSE dla prognozy za pomocą regresji grzbietowej: 143.303113635089"

out = glmnet(x, y, alpha = 0) 
Reg_g_coef = predict(out, type = "coefficients", s = bestlam)[1:18,] # Współczynniki
Reg_g_coef 

##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
## 35.4259636579  3.5360447063  0.0004891764  0.0003725834  0.0001803018 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  0.0929335051  0.1061424028 -0.0001242029 -0.0013135925  0.0007248959 
##    Room.Board         Books      Personal           PhD      Terminal 
##  0.0017740362 -0.0022337145 -0.0018475835  0.0491413814 -0.0171970338 
##     S.F.Ratio   perc.alumni        Expend 
##  0.0425557858  0.2401109479 -0.0001708884

2. Model regresji Lasso

lasso_mod = glmnet(x_train, y_train, alpha = 1,lambda = grid) # Model lasso
plot(lasso_mod)  # Wykres współczynników

## Warning in regularize.values(x, y, ties, missing(ties), na.rm = na.rm):
## collapsing to unique 'x' values

Identyfikacja wartości lambdy, która najbardziej minimalizuje błąd MSE w modelu Lasso.

set.seed(1)
cv.out = cv.glmnet(x_train, y_train, alpha = 1) # Dopasowanie modelu Lasso
bestlam = cv.out$lambda.min # Najlepsza lambda
plot(cv.out) 

Ponownie wyznaczamy wartości prognozowane oraz liczymy błąd MSE. Następnie sprawdzamy współczynniki modelu.

lasso_pred = predict(lasso_mod, s = bestlam, newx = x_test) # Predykcja z najlepszą lambdą
MSE_reg_l = mean((lasso_pred - y_test)^2) # MSE dla Lasso
print(paste("Wartość błędu MSE dla prognozy za pomocą regresji Lasso:",  MSE_reg_l))

## [1] "Wartość błędu MSE dla prognozy za pomocą regresji Lasso: 150.170305331426"

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:18,] # Wszystkie współczynniki
lasso_coef

##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
##  3.601812e+01  1.696896e+00  5.109904e-04  0.000000e+00  0.000000e+00 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  1.668854e-02  1.437067e-01  0.000000e+00 -1.248776e-03  8.922760e-04 
##    Room.Board         Books      Personal           PhD      Terminal 
##  1.452065e-03  0.000000e+00 -1.539697e-03  0.000000e+00  0.000000e+00 
##     S.F.Ratio   perc.alumni        Expend 
##  0.000000e+00  2.523755e-01 -5.372494e-06

Tym razem 7 współczynników zostało wyzerowane

lasso_coef_2 = lasso_coef[lasso_coef != 0] # Wyświetlanie tylko niezerowych współczynników
lasso_coef_2

##   (Intercept)    PrivateYes          Apps     Top10perc     Top25perc 
##  3.601812e+01  1.696896e+00  5.109904e-04  1.668854e-02  1.437067e-01 
##   P.Undergrad      Outstate    Room.Board      Personal   perc.alumni 
## -1.248776e-03  8.922760e-04  1.452065e-03 -1.539697e-03  2.523755e-01 
##        Expend 
## -5.372494e-06

Finalnie pozostaje 11 istotnych współczynników w modelu

3. Klasyczny model OLS

x_train <- as.data.frame(x_train)  # Konwersja macierzy na ramkę danych
y_train <- as.data.frame(y_train) 
train_data <- cbind(y_train, x_train)

ols_mod <- lm(y_train ~ ., data = train_data)  # OLS: wszystkie zmienne objaśniające
summary(ols_mod)  # Wyświetlenie szczegółów modelu

## 
## Call:
## lm(formula = y_train ~ ., data = train_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -51.834  -6.962  -0.523   7.188  54.344 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 30.3348289  5.5530915   5.463 6.86e-08 ***
## PrivateYes   3.9113363  1.9446868   2.011 0.044738 *  
## Apps         0.0010006  0.0005299   1.888 0.059489 .  
## Accept      -0.0004933  0.0010687  -0.462 0.644525    
## Enroll       0.0030146  0.0027903   1.080 0.280416    
## Top10perc    0.0804254  0.0812505   0.990 0.322646    
## Top25perc    0.0817411  0.0613172   1.333 0.183006    
## F.Undergrad -0.0004387  0.0004752  -0.923 0.356260    
## P.Undergrad -0.0016990  0.0004353  -3.903 0.000106 ***
## Outstate     0.0008562  0.0002601   3.292 0.001053 ** 
## Room.Board   0.0022091  0.0006675   3.309 0.000991 ***
## Books       -0.0023413  0.0032875  -0.712 0.476637    
## Personal    -0.0015716  0.0008594  -1.829 0.067947 .  
## PhD          0.0678093  0.0642605   1.055 0.291744    
## Terminal    -0.0252143  0.0708814  -0.356 0.722170    
## S.F.Ratio    0.1375400  0.1806895   0.761 0.446838    
## perc.alumni  0.3113216  0.0558734   5.572 3.80e-08 ***
## Expend      -0.0003238  0.0001660  -1.950 0.051627 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.97 on 604 degrees of freedom
## Multiple R-squared:  0.4598, Adjusted R-squared:  0.4446 
## F-statistic: 30.24 on 17 and 604 DF,  p-value: < 2.2e-16

summary(ols_mod)$coefficients

##                  Estimate   Std. Error    t value     Pr(>|t|)
## (Intercept) 30.3348289190 5.5530914628  5.4626921 6.862151e-08
## PrivateYes   3.9113363224 1.9446868419  2.0112937 4.473820e-02
## Apps         0.0010005716 0.0005299303  1.8881194 5.948895e-02
## Accept      -0.0004933097 0.0010686657 -0.4616127 6.445253e-01
## Enroll       0.0030145501 0.0027903341  1.0803545 2.804155e-01
## Top10perc    0.0804254087 0.0812504950  0.9898452 3.226462e-01
## Top25perc    0.0817410923 0.0613171903  1.3330861 1.830061e-01
## F.Undergrad -0.0004386809 0.0004751626 -0.9232225 3.562599e-01
## P.Undergrad -0.0016989840 0.0004353336 -3.9027175 1.058208e-04
## Outstate     0.0008561899 0.0002600847  3.2919656 1.052934e-03
## Room.Board   0.0022090732 0.0006675431  3.3092591 9.910922e-04
## Books       -0.0023412591 0.0032875211 -0.7121655 4.766373e-01
## Personal    -0.0015715846 0.0008594283 -1.8286396 6.794660e-02
## PhD          0.0678093163 0.0642605140  1.0552252 2.917442e-01
## Terminal    -0.0252143361 0.0708813729 -0.3557258 7.221701e-01
## S.F.Ratio    0.1375400010 0.1806895325  0.7611952 4.468376e-01
## perc.alumni  0.3113216154 0.0558734003  5.5719110 3.798553e-08
## Expend      -0.0003237719 0.0001660288 -1.9500945 5.162711e-02

W kolejnym kroku sekwencyjnie usuwamy nieistotne zmienne objaśniające

# Usuwanie niestotnych zmiennych
stepwise_mod <- step(ols_mod, direction = "backward")

## Start:  AIC=3205.7
## y_train ~ PrivateYes + Apps + Accept + Enroll + Top10perc + Top25perc + 
##     F.Undergrad + P.Undergrad + Outstate + Room.Board + Books + 
##     Personal + PhD + Terminal + S.F.Ratio + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## - Terminal     1      21.3 101634 3203.8
## - Accept       1      35.8 101649 3203.9
## - Books        1      85.3 101698 3204.2
## - S.F.Ratio    1      97.5 101710 3204.3
## - F.Undergrad  1     143.4 101756 3204.6
## - Top10perc    1     164.8 101778 3204.7
## - PhD          1     187.3 101800 3204.8
## - Enroll       1     196.4 101809 3204.9
## - Top25perc    1     299.0 101912 3205.5
## <none>                     101613 3205.7
## - Personal     1     562.6 102175 3207.1
## - Apps         1     599.8 102213 3207.4
## - Expend       1     639.8 102253 3207.6
## - PrivateYes   1     680.6 102293 3207.9
## - Outstate     1    1823.2 103436 3214.8
## - Room.Board   1    1842.4 103455 3214.9
## - P.Undergrad  1    2562.4 104175 3219.2
## - perc.alumni  1    5223.0 106836 3234.9
## 
## Step:  AIC=3203.83
## y_train ~ PrivateYes + Apps + Accept + Enroll + Top10perc + Top25perc + 
##     F.Undergrad + P.Undergrad + Outstate + Room.Board + Books + 
##     Personal + PhD + S.F.Ratio + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## - Accept       1      35.0 101669 3202.0
## - S.F.Ratio    1      98.4 101733 3202.4
## - Books        1     100.1 101734 3202.4
## - F.Undergrad  1     152.3 101787 3202.8
## - Top10perc    1     183.9 101818 3203.0
## - Enroll       1     200.1 101834 3203.1
## - PhD          1     225.5 101860 3203.2
## - Top25perc    1     283.0 101917 3203.6
## <none>                     101634 3203.8
## - Personal     1     553.4 102188 3205.2
## - Apps         1     600.8 102235 3205.5
## - Expend       1     657.8 102292 3205.8
## - PrivateYes   1     710.5 102345 3206.2
## - Outstate     1    1802.2 103436 3212.8
## - Room.Board   1    1821.1 103455 3212.9
## - P.Undergrad  1    2557.1 104191 3217.3
## - perc.alumni  1    5201.9 106836 3232.9
## 
## Step:  AIC=3202.05
## y_train ~ PrivateYes + Apps + Enroll + Top10perc + Top25perc + 
##     F.Undergrad + P.Undergrad + Outstate + Room.Board + Books + 
##     Personal + PhD + S.F.Ratio + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## - Books        1      95.7 101765 3200.6
## - S.F.Ratio    1     103.4 101773 3200.7
## - F.Undergrad  1     144.1 101813 3200.9
## - Enroll       1     165.8 101835 3201.1
## - PhD          1     209.9 101879 3201.3
## - Top25perc    1     259.8 101929 3201.6
## - Top10perc    1     282.8 101952 3201.8
## <none>                     101669 3202.0
## - Personal     1     542.0 102211 3203.4
## - Expend       1     629.6 102299 3203.9
## - PrivateYes   1     686.0 102355 3204.2
## - Apps         1    1314.0 102983 3208.0
## - Outstate     1    1772.2 103441 3210.8
## - Room.Board   1    1829.8 103499 3211.1
## - P.Undergrad  1    2522.5 104192 3215.3
## - perc.alumni  1    5312.2 106981 3231.7
## 
## Step:  AIC=3200.63
## y_train ~ PrivateYes + Apps + Enroll + Top10perc + Top25perc + 
##     F.Undergrad + P.Undergrad + Outstate + Room.Board + Personal + 
##     PhD + S.F.Ratio + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## - S.F.Ratio    1      99.2 101864 3199.2
## - F.Undergrad  1     141.3 101906 3199.5
## - Enroll       1     160.6 101926 3199.6
## - Top25perc    1     245.6 102011 3200.1
## - PhD          1     250.3 102015 3200.2
## - Top10perc    1     275.9 102041 3200.3
## <none>                     101765 3200.6
## - Personal     1     630.4 102395 3202.5
## - Expend       1     663.7 102429 3202.7
## - PrivateYes   1     685.8 102451 3202.8
## - Apps         1    1315.4 103080 3206.6
## - Room.Board   1    1754.1 103519 3209.3
## - Outstate     1    1778.5 103543 3209.4
## - P.Undergrad  1    2538.9 104304 3214.0
## - perc.alumni  1    5438.3 107203 3231.0
## 
## Step:  AIC=3199.24
## y_train ~ PrivateYes + Apps + Enroll + Top10perc + Top25perc + 
##     F.Undergrad + P.Undergrad + Outstate + Room.Board + Personal + 
##     PhD + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## - F.Undergrad  1     136.8 102001 3198.1
## - Enroll       1     163.1 102027 3198.2
## - Top25perc    1     238.7 102103 3198.7
## - PhD          1     262.6 102127 3198.8
## - Top10perc    1     275.1 102139 3198.9
## <none>                     101864 3199.2
## - PrivateYes   1     623.5 102488 3201.0
## - Personal     1     683.3 102547 3201.4
## - Expend       1    1012.9 102877 3203.4
## - Apps         1    1347.7 103212 3205.4
## - Outstate     1    1732.8 103597 3207.7
## - Room.Board   1    1740.8 103605 3207.8
## - P.Undergrad  1    2534.2 104398 3212.5
## - perc.alumni  1    5358.4 107223 3229.1
## 
## Step:  AIC=3198.07
## y_train ~ PrivateYes + Apps + Enroll + Top10perc + Top25perc + 
##     P.Undergrad + Outstate + Room.Board + Personal + PhD + perc.alumni + 
##     Expend
## 
##               Df Sum of Sq    RSS    AIC
## - Enroll       1      29.3 102030 3196.3
## - Top25perc    1     232.9 102234 3197.5
## - PhD          1     253.8 102255 3197.6
## - Top10perc    1     268.8 102270 3197.7
## <none>                     102001 3198.1
## - Personal     1     747.8 102749 3200.6
## - PrivateYes   1     762.5 102763 3200.7
## - Expend       1     988.4 102989 3202.1
## - Apps         1    1292.2 103293 3203.9
## - Room.Board   1    1721.7 103723 3206.5
## - Outstate     1    1776.4 103777 3206.8
## - P.Undergrad  1    2972.5 104973 3213.9
## - perc.alumni  1    5502.9 107504 3228.8
## 
## Step:  AIC=3196.25
## y_train ~ PrivateYes + Apps + Top10perc + Top25perc + P.Undergrad + 
##     Outstate + Room.Board + Personal + PhD + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## - Top25perc    1     241.7 102272 3195.7
## - Top10perc    1     263.9 102294 3195.9
## - PhD          1     268.1 102298 3195.9
## <none>                     102030 3196.3
## - Personal     1     724.9 102755 3198.7
## - PrivateYes   1     733.5 102764 3198.7
## - Expend       1    1026.2 103056 3200.5
## - Room.Board   1    1694.4 103725 3204.5
## - Outstate     1    1792.6 103823 3205.1
## - P.Undergrad  1    3004.2 105034 3212.3
## - Apps         1    4005.6 106036 3218.2
## - perc.alumni  1    5529.0 107559 3227.1
## 
## Step:  AIC=3195.72
## y_train ~ PrivateYes + Apps + Top10perc + P.Undergrad + Outstate + 
##     Room.Board + Personal + PhD + perc.alumni + Expend
## 
##               Df Sum of Sq    RSS    AIC
## <none>                     102272 3195.7
## - PhD          1     370.6 102643 3196.0
## - PrivateYes   1     719.4 102991 3198.1
## - Personal     1     736.3 103008 3198.2
## - Expend       1    1308.2 103580 3201.6
## - Room.Board   1    1725.5 103997 3204.1
## - Outstate     1    1820.9 104093 3204.7
## - Top10perc    1    2118.2 104390 3206.5
## - P.Undergrad  1    2981.6 105254 3211.6
## - Apps         1    4216.9 106489 3218.9
## - perc.alumni  1    5735.5 108007 3227.7

# Predykcja na danych testowych
x_test <- as.data.frame(x_test)  # Konwersja macierzy testowej na ramkę danych
ols_pred <- predict(stepwise_mod, newdata = x_test)
summary(stepwise_mod)$coefficients

##                  Estimate   Std. Error   t value     Pr(>|t|)
## (Intercept) 33.3222467245 3.4747006252  9.589962 2.201419e-20
## PrivateYes   3.7740985934 1.8204574735  2.073159 3.857600e-02
## Apps         0.0008993941 0.0001791899  5.019223 6.814264e-07
## Top10perc    0.1647836641 0.0463218480  3.557364 4.036245e-04
## P.Undergrad -0.0017216581 0.0004079227 -4.220550 2.807359e-05
## Outstate     0.0008331982 0.0002526161  3.298278 1.029290e-03
## Room.Board   0.0020796857 0.0006477270  3.210744 1.393386e-03
## Personal    -0.0017470704 0.0008329927 -2.097342 3.637366e-02
## PhD          0.0636885579 0.0428038725  1.487916 1.372890e-01
## perc.alumni  0.3207703658 0.0547984130  5.853643 7.851020e-09
## Expend      -0.0004122359 0.0001474575 -2.795626 5.342737e-03

OLS_coef=stepwise_mod$coefficients

Jak widać, finalnie również pozostało 11 istotnych zmiennych objaśniających. Następnie liczymy błąd MSE.

MSE_ols= mean((ols_pred - y_test)^2) # MSE dla OLS
print(paste("Wartość błędu MSE dla prognozy za pomocą OLS:",  MSE_ols))

## [1] "Wartość błędu MSE dla prognozy za pomocą OLS: 147.770118442805"

plot(stepwise_mod)

Wyniki

Na końcu należy porównać uzyskane wyniki za pomocą wyznaczonych miar błędu MSE.

data.frame(Model=c("Regresja grzbietowa", "Regresja Lasso", "Klasyczny OLS"),MSE=c(MSE_reg_g,MSE_reg_l,MSE_ols))

##                 Model      MSE
## 1 Regresja grzbietowa 143.3031
## 2      Regresja Lasso 150.1703
## 3       Klasyczny OLS 147.7701

Na podstawie wskzanych wartości można powiedzieć, że najlepszy rezultat prognozy uzyskano za pomocą regresji grzbietowej. Należy jednak dodać, że w tym przypadku zostały uwzględnione wszystkie predyktory, co może skutkować overfittingiem. Na drugim miejscu znalazła się klasyczna metoda najmniejszych kwadratów, natomiast wynik nie odbiega znacznie od regresji Lasso. Dodatkowo w metodzie OLS nie zostały zweryfikowane wszelkie testy związane z klasycznymi założeniami KNMK, które prawdopodbnie nie zostały całkowicie spełnione.

Predyktory

• Które predyktory okazały się ważne w ostatecznym modelu (modelach)?

list(Reg_g_coef, lasso_coef_2, OLS_coef)

## [[1]]
##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
## 35.4259636579  3.5360447063  0.0004891764  0.0003725834  0.0001803018 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  0.0929335051  0.1061424028 -0.0001242029 -0.0013135925  0.0007248959 
##    Room.Board         Books      Personal           PhD      Terminal 
##  0.0017740362 -0.0022337145 -0.0018475835  0.0491413814 -0.0171970338 
##     S.F.Ratio   perc.alumni        Expend 
##  0.0425557858  0.2401109479 -0.0001708884 
## 
## [[2]]
##   (Intercept)    PrivateYes          Apps     Top10perc     Top25perc 
##  3.601812e+01  1.696896e+00  5.109904e-04  1.668854e-02  1.437067e-01 
##   P.Undergrad      Outstate    Room.Board      Personal   perc.alumni 
## -1.248776e-03  8.922760e-04  1.452065e-03 -1.539697e-03  2.523755e-01 
##        Expend 
## -5.372494e-06 
## 
## [[3]]
##   (Intercept)    PrivateYes          Apps     Top10perc   P.Undergrad 
## 33.3222467245  3.7740985934  0.0008993941  0.1647836641 -0.0017216581 
##      Outstate    Room.Board      Personal           PhD   perc.alumni 
##  0.0008331982  0.0020796857 -0.0017470704  0.0636885579  0.3207703658 
##        Expend 
## -0.0004122359

W modeli regresji grzbietowej nie została odrzucona żadna zmienna. W przypadku modelu regresji Lasso oraz OLS, zostało po 11 istotnych zmiennych. Oba modele różnią się wyłącznie 2 zmiennymi (Top25perc oraz PhD).