Nieklasyczne metody statystyki

Regularyzacja

library(ISLR)
library(dplyr)
library(glmnet)
data(College)

College<-na.omit(College)

x = model.matrix(Grad.Rate~., College)[,-1] 
                                        
y = College %>%
  select(Grad.Rate) %>%
  unlist() %>%
  as.numeric()

Analizowaną zmienną w tym raporcie jest zmienna Grad.Rate ze zbioru danych College z pakietu ISLR

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

Wartość alpha została w tym przypadku ustalona na 0 co oznacza wybór metody regresji grzbietowej

dim(coef(ridge_mod))

## [1]  18 100

plot(ridge_mod) 

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)    PrivateYes          Apps        Accept        Enroll 
##  6.524147e+01  1.924300e-02  9.686567e-07  7.024066e-07 -6.112685e-07 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  7.153334e-04  6.143038e-04 -4.144832e-07 -4.310351e-06  3.619870e-06 
##    Room.Board         Books      Personal           PhD      Terminal 
##  9.873884e-06  1.194325e-07 -1.014868e-05  4.754458e-04  5.000512e-04 
##     S.F.Ratio   perc.alumni        Expend 
## -1.970232e-03  1.009748e-03  1.899832e-06

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

## [1] 0.01940515

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)    PrivateYes          Apps        Accept        Enroll 
##  6.222586e+01  2.855277e-01  1.464045e-05  1.093961e-05 -8.384396e-06 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  1.046539e-02  9.031997e-03 -6.110823e-06 -6.531529e-05  5.334163e-05 
##    Room.Board         Books      Personal           PhD      Terminal 
##  1.453208e-04 -8.045918e-06 -1.524621e-04  6.855198e-03  7.141166e-03 
##     S.F.Ratio   perc.alumni        Expend 
## -2.809567e-02  1.502890e-02  2.710929e-05

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

## [1] 0.2878028

Szacujemy współczynniki dla lambdy wynoszącej 50

predict(ridge_mod, s = 50, type = "coefficients")[1:12,]

##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
##  4.671830e+01  1.903498e+00  1.234870e-04  1.146584e-04  3.405181e-06 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  6.176559e-02  5.629912e-02 -3.664026e-05 -5.169350e-04  3.385363e-04 
##    Room.Board         Books 
##  9.118980e-04 -5.822815e-04

Dzielimy dane na zbiór treningowy i testowy

set.seed(1)

train = College %>%
  sample_frac(0.7)

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()

Dla lambdy równej 4 otrzymujemy błąd MSE wynaszący 175,5626

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] 175.5626

Błąd ten dla modelu z samym wyrazem wolnym wyniósłby 288,5538

mean((mean(y_train) - y_test)^2)

## [1] 288.5538

Ten sam wynik uzyskujemy również przy modelu regresji grzbietowej z bardzo wysoką wartością lambda

ridge_pred = predict(ridge_mod, s = 1e10, newx = x_test)
mean((ridge_pred - y_test)^2)

## [1] 288.5538

Dla modelu klasyczną metodą najmniejszych kwadratów wartość błędu MSE wynosi 175,1703, a więc jest on mniejszy niż dla poprzednich modeli z regresją grzbietową

ridge_pred = predict(ridge_mod, s = 0, newx = x_test)
mean((ridge_pred - y_test)^2)

## [1] 175.1703

lm(Grad.Rate~., data = train)

## 
## Call:
## lm(formula = Grad.Rate ~ ., data = train)
## 
## Coefficients:
## (Intercept)   PrivateYes         Apps       Accept       Enroll    Top10perc  
##  32.3042366    4.1341442    0.0010918   -0.0008182    0.0033567    0.1151044  
##   Top25perc  F.Undergrad  P.Undergrad     Outstate   Room.Board        Books  
##   0.0692132   -0.0003865   -0.0016357    0.0008178    0.0023210   -0.0037467  
##    Personal          PhD     Terminal    S.F.Ratio  perc.alumni       Expend  
##  -0.0024720    0.0633880   -0.0368708    0.2048758    0.3086942   -0.0003419

predict(ridge_mod, s = 0, type="coefficients")[1:12,]

##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
## 32.3089038088  4.1316615324  0.0010763654 -0.0007834882  0.0032723816 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  0.1158221837  0.0690767008 -0.0003767153 -0.0016347369  0.0008157950 
##    Room.Board         Books 
##  0.0023202998 -0.0037430587

Szukamy wartości dla lambdy która da nam najniższą wartość błędu MSE

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] 3.680579

Najniższy błąd MSE będzie miała lambda o wartości 3,680579

plot(cv.out) # Narysuj wykres treningowego MSE jako funkcję lambda

Jak widać na wykresie wartość błędu MSE rośnie wraz ze wzrostem wartości Lambdy w okolicach wartości logarytmu z Lambdy wynoszącego około 2

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] 175.5099

out = glmnet(x, y, alpha = 0) # Dopasuj model regresji grzbietowej do pełnego zbioru danych
predict(out, type = "coefficients", s = bestlam)[1:12,] # Wyświetlanie współczynników przy użyciu lambda wybranego przez CV

##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
## 35.5843548058  3.5253980915  0.0004716526  0.0003717433  0.0001752357 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  0.0928485480  0.1050442817 -0.0001180038 -0.0012969894  0.0007123055 
##    Room.Board         Books 
##  0.0017566454 -0.0022095533

Dla najlepszej Lambdy wartość MSE nadal jest większa niż przy metodzie najmniejszych kwadratów

Następnie podobne kroki przeprowadzamy dla metody Lasso (Alpha = 1)

#lasso

lasso_mod = glmnet(x_train, 
                   y_train, 
                   alpha = 1, 
                   lambda = grid) # Dopasuj model lasso do danych treningowych

plot(lasso_mod) 

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] 174.7583

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:12,] # Wyświetlanie współczynników przy użyciu lambda wybranego przez CV
lasso_coef

##   (Intercept)    PrivateYes          Apps        Accept        Enroll 
##  3.465761e+01  3.067055e+00  7.727418e-04  0.000000e+00  0.000000e+00 
##     Top10perc     Top25perc   F.Undergrad   P.Undergrad      Outstate 
##  6.028698e-02  1.263351e-01 -1.948574e-05 -1.501766e-03  9.362671e-04 
##    Room.Board         Books 
##  1.734008e-03 -1.649290e-03

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

##   (Intercept)    PrivateYes          Apps     Top10perc     Top25perc 
##  3.465761e+01  3.067055e+00  7.727418e-04  6.028698e-02  1.263351e-01 
##   F.Undergrad   P.Undergrad      Outstate    Room.Board         Books 
## -1.948574e-05 -1.501766e-03  9.362671e-04  1.734008e-03 -1.649290e-03

Dla metody LASSO błąd MSE wyniósł 174,7583 a więc jest on mniejszy zarówno od błędu przy regresji grzbietowej jak i przy najmniejszych kwadratach.