Nieklasyczne metody statystyki

Regularyzacja

Twoja kolej!

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

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

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

Credit=na.omit(Credit)
x = model.matrix(Balance~., Credit)[,-1] # przycinam pierwszą kolumnę
                                         # zostawiam predyktory
y = Credit %>%
  select(Balance) %>%
  unlist() %>%
  as.numeric()

modellm <- lm( y~x)
summary(modellm)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -166.48  -77.62  -14.37   56.21  316.52 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -487.07424   36.73407 -13.259  < 2e-16 ***
## xID                    0.04105    0.04343   0.945   0.3452    
## xIncome               -7.80740    0.23431 -33.321  < 2e-16 ***
## xLimit                 0.19052    0.03279   5.811  1.3e-08 ***
## xRating                1.14249    0.49100   2.327   0.0205 *  
## xCards                17.83639    4.34324   4.107  4.9e-05 ***
## xAge                  -0.62955    0.29449  -2.138   0.0332 *  
## xEducation            -1.09831    1.59817  -0.687   0.4924    
## xGenderFemale         -9.54615    9.98431  -0.956   0.3396    
## xStudentYes          426.16715   16.73077  25.472  < 2e-16 ***
## xMarriedYes           -8.78055   10.36758  -0.847   0.3976    
## xEthnicityAsian       16.85752   14.12112   1.194   0.2333    
## xEthnicityCaucasian    9.29289   12.24194   0.759   0.4483    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 98.8 on 387 degrees of freedom
## Multiple R-squared:  0.9552, Adjusted R-squared:  0.9538 
## F-statistic: 687.7 on 12 and 387 DF,  p-value: < 2.2e-16

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

dim(coef(ridge_mod))

## [1]  13 100

plot(ridge_mod) # wykres współczynników

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)                 ID             Income              Limit 
##       4.437406e+02       7.098518e-04       2.072944e-01       6.255481e-03 
##             Rating              Cards                Age          Education 
##       9.352902e-02       1.094867e+00      -7.433736e-03      -3.648999e-02 
##       GenderFemale         StudentYes         MarriedYes     EthnicityAsian 
##       7.279655e-01       1.522470e+01      -2.740740e-01      -3.232180e-01 
## EthnicityCaucasian 
##      -8.954612e-02

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

## [1] 15.28924

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

##    (Intercept)             ID         Income          Limit         Rating 
##  -387.01506763     0.02788993    -4.71033241     0.11026631     1.60846689 
##          Cards            Age      Education   GenderFemale     StudentYes 
##    16.10495492    -1.00880354    -0.40547730    -3.16706653   372.95067139 
##     MarriedYes EthnicityAsian 
##   -12.41062933    12.27125777

set.seed(1)

train = Credit %>%
  sample_frac(0.5)

test = Credit %>%
  setdiff(train)

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

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

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

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

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

## [1] 194030.9

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

## [1] 10660.7

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

#Widzimy zatem, że wartość λ która powoduje najmniejszy błąd walidacji krzyżowej to 42

plot(cv.out) # Narysuj wykres treningowego MSE jako funkcję 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] 16116.16

out = glmnet(x, y, alpha = 0)
predict(out, type = "coefficients", s = bestlam)[1:13,]

##        (Intercept)                 ID             Income              Limit 
##      -402.69174822         0.02968823        -5.08362998         0.11363315 
##             Rating              Cards                Age          Education 
##         1.65055380        15.94593669        -0.97611320        -0.46045505 
##       GenderFemale         StudentYes         MarriedYes     EthnicityAsian 
##        -3.91626143       380.23436720       -12.29347706        13.01622062 
## EthnicityCaucasian 
##         8.53524735

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

plot(lasso_mod)    # Wykreśl współczynniki

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

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

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

##        (Intercept)                 ID             Income              Limit 
##      -490.32404851         0.03266285        -7.67489924         0.17202850 
##             Rating              Cards                Age          Education 
##         1.38723835        16.02401187        -0.59472497        -0.72481881 
##       GenderFemale         StudentYes         MarriedYes     EthnicityAsian 
##        -7.47592875       421.49160417        -7.03306237        11.52382785 
## EthnicityCaucasian 
##         4.72253992

Aby zaliczyć to laboratorium, zamieść odpowiedzi na następujące pytania:

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

• Jaka była Twoja zmienna zależna (tzn. co próbowałeś modelować)? - Średnie saldo karty kredytowej w $

• 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). Oczkekiwaliśmy, że lepszą będzie regresja grzbietowa, ze względu na możliwość występowania korelacji w predyktorach. W modelu OLS R^2 wyniosło 0,955 ( R-squared: 0.9552), co oznacza, że model jest bardzo dobrze dopasowany. MSE dla regresji grzbietowej wyniosi 16116, a dla regresji LASSO 10506, co stanowi mniejszą wartość. W przypadku regresji grzbietowej żaden ze wsppółczynnikow nie jest dokładnie zerowy, regresja ta nie dokonuje selekcji zmiennych. W przypadku LASSO można zauważyć, że występują współczynniki, które wynoszą dokładnie 0. Za model optymalny uznano zatem model regresji LASSO, ze względu na najniższą wartość MSE.

• Które predyktory okazały się ważne w ostatecznym modelu (modelach)? Wszystkie predyktory okazały się być ważne, jednak najważniejszym było posiadanie statusu studenta.