Publicacion_lasso_ridge

11 minute read

Regresion Ridge VS Lasso

Regresión con regularización Lasso y Ridge

Cuando realizamos regresión lineal podemos en algunos casos querer ver la manera de regularizar nuestro modelo a manera de que nuestro modelo no sobreajuste o también con la finalidad de mejorar nuestra predicción en conjuntos de prueba

Otro de los puntos que también se toca es la seleccion de variables que formara nuestro modelo, existen varios metodos:

  • Stepwise Forward Regression
  • Stepwise Backward Regression
  • Backward-Forward Regresion (la mezcla de las 2 previas)
  • Exhastive Variable Selection

Para la seleccion del subconjunto de variables que mejor funciona podemos aplicar la regularizacion con penalizacion L1 también llamado Regularización Lasso, este método permite hacer la selección de variables reduciendo a cero las variables de poca importancia para el modelo y ademas funciona con mayor velocidad que los metodos previamente descritos

En cuanto al otro metodo de regularización Ridge (Penalizacion L2) no podemos categorizarlo como un metodo de selección de variables dado que aun y que hace que ciertas variables se encojan en su peso, no garantiza llevarla a ceros, por lo cual, siempre que se aplica el modelo resultante es un modelo con el mismo numero de variables de inicio, pero en el cual las variables que no aportan al modelo de regresión se encogen en su coeficiente

Las 2 formulaciones teoricas de las regresiones con penalizacion L1 y L2 (Lasso y Ridge) son las siguientes:

  • Lasso \(\sum_{i=1}^{n}(y_i-\beta_0-\sum_{j=1}^{p}\beta_jx_{ij})^2 + \lambda\sum_{j=1}^{p}|\beta_j| = RSS + \lambda\sum_{j=1}^{p}|\beta_j|\)

  • Ridge \(\sum_{i=1}^{n}(y_i-\beta_0-\sum_{j=1}^{p}\beta_jx_{ij})^2 + \lambda\sum_{j=1}^{p}\beta_j^2 = RSS + \lambda\sum_{j=1}^{p}(\beta_j)^2\)

Cuando hacer uso de un metodo de selección de variables o de regualarización

Dependiendo el objetivo del procedimiento que vayamos a realizar se podría sugerir un metodo con repecto a otro, y también dependiendo del tiempo que tengamos, dado que hay metodos muy costosos computacionalmente hablando

Si nuestro objetivo es realizar un forecasting quizas la mejor opcion sea probar varios metodos realizar validación cruzada y quedarnos con el que menor valor reporté en la funcion de perdida que elijamos (MSE,RMSE,etc…)

Si nuestro objetivo es interpretación de entrada quizas no sea buena opcción el metodo Ridge, dado que si tenemos un problema con muchas variables al finalizar nuestro modelo tendremos el mismo numero y será dificil de poder interpretar, aqui serías mas conveniente ver alguno de los metodos de selección de variables, aunque si es un problema con muchas variables no será buena opción usar el metodo exhaustivo (Revisar todos los modelos de subconjuntos posibles de variables \(2^n\) modelos posibles) quizas viendolo en tiempo y con la finalidad de interpretación se optaria por Lasso

El problema a analizar

Para revisar los 2 metodos de regularización Lasso y Ridge checare como funciona sobre un conjunto de datos, quizas no es posible con solo un conjunto de datos decir que uno le ganaa al otro pero la idea es ver como funcionan

Los datos a revisar es el conjunto de datos Carseats de la librería ISLR, como nota esta librería pertenece a un excelente libro Introduction to Statical Learning with R el cual recomiendo ampliamente

Vamos a explorar un poco el conjunto de datos que analizaremos

library(ISLR)
library(leaps)
library(glmnet)
## Loading required package: Matrix
## Loading required package: foreach
## Loaded glmnet 2.0-16
data("Carseats")
##Que variables tenemos
colnames(Carseats)
##  [1] "Sales"       "CompPrice"   "Income"      "Advertising" "Population" 
##  [6] "Price"       "ShelveLoc"   "Age"         "Education"   "Urban"      
## [11] "US"
## tipos de datos de las variables
sapply(Carseats, class)
##       Sales   CompPrice      Income Advertising  Population       Price 
##   "numeric"   "numeric"   "numeric"   "numeric"   "numeric"   "numeric" 
##   ShelveLoc         Age   Education       Urban          US 
##    "factor"   "numeric"   "numeric"    "factor"    "factor"
tipos_variables <-sapply(Carseats, class) 
variables_numericas <- names(tipos_variables[tipos_variables=="numeric"])

# un summary de las variables numericas para ver sus rangos medias etc...
summary(Carseats)
##      Sales          CompPrice       Income        Advertising    
##  Min.   : 0.000   Min.   : 77   Min.   : 21.00   Min.   : 0.000  
##  1st Qu.: 5.390   1st Qu.:115   1st Qu.: 42.75   1st Qu.: 0.000  
##  Median : 7.490   Median :125   Median : 69.00   Median : 5.000  
##  Mean   : 7.496   Mean   :125   Mean   : 68.66   Mean   : 6.635  
##  3rd Qu.: 9.320   3rd Qu.:135   3rd Qu.: 91.00   3rd Qu.:12.000  
##  Max.   :16.270   Max.   :175   Max.   :120.00   Max.   :29.000  
##    Population        Price        ShelveLoc        Age       
##  Min.   : 10.0   Min.   : 24.0   Bad   : 96   Min.   :25.00  
##  1st Qu.:139.0   1st Qu.:100.0   Good  : 85   1st Qu.:39.75  
##  Median :272.0   Median :117.0   Medium:219   Median :54.50  
##  Mean   :264.8   Mean   :115.8                Mean   :53.32  
##  3rd Qu.:398.5   3rd Qu.:131.0                3rd Qu.:66.00  
##  Max.   :509.0   Max.   :191.0                Max.   :80.00  
##    Education    Urban       US     
##  Min.   :10.0   No :118   No :142  
##  1st Qu.:12.0   Yes:282   Yes:258  
##  Median :14.0                      
##  Mean   :13.9                      
##  3rd Qu.:16.0                      
##  Max.   :18.0
#Ver si existen relaciones cercanamente lineales
pairs(Carseats[variables_numericas])

Primero usando mejor subconjunto (Exahustivo)

Dado que tenemos pocas variables, podemos darnos el lujo de realizar la seleccion del mejor subconjunto de variables, checando los \(2^p\) submodelos posibles

library(ISLR)
library(leaps)
library(glmnet)

data("Carseats")

x<-model.matrix(Sales~.,data = Carseats)[,-1]

# CReamos un nuevo data frame con las variables categoricas convertidas a dummies
Carseats_new<-data.frame(Sales=Carseats$Sales,x)
colnames(Carseats_new)
##  [1] "Sales"           "CompPrice"       "Income"         
##  [4] "Advertising"     "Population"      "Price"          
##  [7] "ShelveLocGood"   "ShelveLocMedium" "Age"            
## [10] "Education"       "UrbanYes"        "USYes"
## tomar 80% para entrenamiento y 20 para test (Los porcentajes mas comunes 80-20 o 70-30)
set.seed(123)
train<-sample(nrow(Carseats_new),nrow(Carseats_new)*.8)
#train

####  Mejor subconjunto

mejores_modelos_ex<-regsubsets(Sales~.,data =Carseats_new,subset=train,nvmax = 11,method =  "exhaustive",nbest = 1)
summ_ex<-summary(mejores_modelos_ex)
summ_ex$adjr2
##  [1] 0.2292493 0.4664869 0.6056609 0.6940840 0.7722060 0.8459569 0.8697723
##  [8] 0.8699916 0.8698297 0.8695204 0.8691000
summ_ex$cp
##  [1] 1556.412034  978.006512  639.956946  426.161481  238.427109
##  [6]   62.338356    6.397615    6.881785    8.271981   10.007667
## [11]   12.000000
summ_ex$bic
##  [1]  -72.79295 -185.75456 -277.72445 -354.21863 -443.82585 -564.26520
##  [7] -613.26419 -609.06242 -603.92648 -598.43264 -592.67229
#summ_ex$which
MSE_trainning<-summ_ex$rss/nrow(Carseats_new[train,])

best_modelo<-which.min(summ_ex$bic)
which.min(summ_ex$cp)
## [1] 7
which.max(summ_ex$adjr2)
## [1] 8
## el mejor submodelo considerando bic seria de 7 variables predictoras

library(ggplot2)
p <- ggplot(data = data.frame(n_predictores = 1:11,
                              bic_ex = summ_ex$bic[1:11]),
            aes(x = n_predictores, y = bic_ex)) +
  geom_line(color="green") +
  geom_point(color="green")+
  xlab("# Predictores")+
  ylab("Valor BIC")
p

Pero el mejor modelo segun los datos de prueba seria revisar los errores de test de cada uno Haciendo enfasis en lo solicitado se revisara el modelo con mejor error de prueba y no mejor error de entranmiento.

val.errors = rep(NA, 11)
x.test = model.matrix(Sales ~ ., data = Carseats_new[-train,])  # notice the -index!
for (i in 1:11) {
  coefi = coef(mejores_modelos_ex, id = i)
  pred = x.test[, names(coefi)] %*% coefi
  val.errors[i] = mean((Carseats_new$Sales[-train] - pred)^2)
}

val.errors
##  [1] 5.436022 4.265773 2.347872 1.861534 1.720193 1.256614 1.058922
##  [8] 1.064591 1.058951 1.046329 1.048199
n_predictores<-which.min(val.errors)
error_bestsubset<-min(val.errors)
cbind(n_predictores,error_bestsubset)
##      n_predictores error_bestsubset
## [1,]            10         1.046329

Si nos basamos en el MSE de prueba para los 11 modelos posibles, el mejor submodelo es el que tiene 10 variables aunque la diferencia entre los modelos de mas de 7 variables es minima, si queremos un modelo parsimonioso lo mejor seria quedarse con un modelo de 7 variables.

coeficientes_mejor_subconjunto<-coef(mejores_modelos_ex, id = 10)
coeficientes_mejor_subconjunto
##     (Intercept)       CompPrice          Income     Advertising 
##      5.80345180      0.09301823      0.01552543      0.12281461 
##           Price   ShelveLocGood ShelveLocMedium             Age 
##     -0.09544900      4.89251432      2.03492752     -0.04805947 
##       Education        UrbanYes           USYes 
##     -0.02754300      0.10035768     -0.08633875

Ahora probare con Ridge y Lasso

cv.out_ridge<-cv.glmnet(x[train,],Carseats_new$Sales[train],alpha=0)
cv.out_lasso<-cv.glmnet(x[train,],Carseats_new$Sales[train],alpha=1)

plot(cv.out_ridge)

plot(cv.out_lasso)

bestlam_ridge<-cv.out_ridge$lambda.min
bestlam_ridge_se<-cv.out_ridge$lambda.1se

bestlam_lasso<-cv.out_lasso$lambda.min
bestlam_lasso_se<-cv.out_lasso$lambda.1se

### Tanto para lasso como para ridge ,se puede considerar el lamba a un error estandard, esto con el fin
#### DE lograr un modelo aun mas parsimonioso, en el caso de lasso puede reducir aun mas el numero de variables
## sin sacrificar mucho en cuanto poder de prediccion

#entrenar modelo final con lambdas de 1 desviacion estandar
cv.out_lasso$lambda.1se
## [1] 0.09135486
lasso_1se<-glmnet(x[train,],Carseats_new$Sales[train],alpha=1,lambda=cv.out_lasso$lambda.1se)
coef(lasso_1se)
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                          s0
## (Intercept)      6.73364438
## CompPrice        0.07664656
## Income           0.01163099
## Advertising      0.10436310
## Population       .         
## Price           -0.08504838
## ShelveLocGood    4.36087400
## ShelveLocMedium  1.57627189
## Age             -0.04259526
## Education        .         
## UrbanYes         .         
## USYes            .
cv.out_lasso$lambda.min
## [1] 0.01878724
lasso_lmin<-glmnet(x[train,],Carseats_new$Sales[train],alpha=1,lambda=cv.out_lasso$lambda.min)
coef(lasso_lmin)
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                          s0
## (Intercept)      5.98566562
## CompPrice        0.08948958
## Income           0.01472203
## Advertising      0.11572318
## Population       .         
## Price           -0.09333222
## ShelveLocGood    4.78270630
## ShelveLocMedium  1.94021076
## Age             -0.04704128
## Education       -0.01969634
## UrbanYes         0.05853327
## USYes            .
cv.out_ridge$lambda.min
## [1] 0.1488857
ridge_lmin<-glmnet(x[train,],Carseats_new$Sales[train],alpha=0,lambda=cv.out_ridge$lambda.min)
coef(ridge_lmin)
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                            s0
## (Intercept)      6.6079168209
## CompPrice        0.0792225422
## Income           0.0142140836
## Advertising      0.1108983012
## Population      -0.0001026325
## Price           -0.0853374277
## ShelveLocGood    4.4308501551
## ShelveLocMedium  1.7279402001
## Age             -0.0459118570
## Education       -0.0254045819
## UrbanYes         0.0883677350
## USYes            0.0191305325
cv.out_ridge$lambda.1se
## [1] 0.2370681
ridge_1se<-glmnet(x[train,],Carseats_new$Sales[train],alpha=0,lambda=cv.out_ridge$lambda.1se)
coef(ridge_1se)
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                            s0
## (Intercept)      6.9649835006
## CompPrice        0.0726656417
## Income           0.0135650752
## Advertising      0.1051008796
## Population      -0.0001279344
## Price           -0.0803871937
## ShelveLocGood    4.2010786335
## ShelveLocMedium  1.5795065119
## Age             -0.0446699994
## Education       -0.0244412027
## UrbanYes         0.0828909457
## USYes            0.0653223192
erro_lasso_1se<-mean( (predict(lasso_1se,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_lasso_1se
## [1] 1.083014
erro_lasso_lmin<-mean( (predict(lasso_lmin,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_lasso_lmin
## [1] 1.041594
erro_ridge_1se<-mean( (predict(ridge_1se,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_ridge_1se
## [1] 1.156963
erro_ridge_lmin<-mean( (predict(ridge_lmin,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_ridge_lmin
## [1] 1.083586
plot(data.frame(c("lasso_1se","lasso_min","ridge_1se","ridge_min","exhaust"),c(erro_lasso_1se,erro_lasso_lmin,erro_ridge_1se,erro_ridge_lmin,error_bestsubset)),type="o",xlab="",ylab="MSE",col=c("red","green","blue","yellow","purple"))

La ultima grafica mostrada muestra los errores de prueba de los modelos, notese que el mejor en cuanto a prediccion es el modelo de Lasso tomando el lambda (que minimiza),

A simple vista dado que lasso redujo 2 variables y obtuvo mejor MSE de prueba que el exhaustivo que solo redujo una variables, pues la idea es tomar el mas parsimonioso y a la vez mejor en predicción

Aqui la comparativa de los coeficientes

coef(lasso_lmin)
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                          s0
## (Intercept)      5.98566562
## CompPrice        0.08948958
## Income           0.01472203
## Advertising      0.11572318
## Population       .         
## Price           -0.09333222
## ShelveLocGood    4.78270630
## ShelveLocMedium  1.94021076
## Age             -0.04704128
## Education       -0.01969634
## UrbanYes         0.05853327
## USYes            .
coef(mejores_modelos_ex,id = 10)
##     (Intercept)       CompPrice          Income     Advertising 
##      5.80345180      0.09301823      0.01552543      0.12281461 
##           Price   ShelveLocGood ShelveLocMedium             Age 
##     -0.09544900      4.89251432      2.03492752     -0.04805947 
##       Education        UrbanYes           USYes 
##     -0.02754300      0.10035768     -0.08633875

Otra cosa por notar es que el modelo de lasso tomando el lambda de un error estandard es comparable con el Ridge de lambda minimo y este modelo lasso_1se reduce 4 variables a cero,

Yo en mi caso eligiria el modelo Lasso que usa el lambda de un errror estandar y que elimina 4 variables del modelo

PCR

Veamos si PCR logra algo que sea mejor en MSE de prueba como para que valga la pena considerarlo

library(pls)
## 
## Attaching package: 'pls'
## The following object is masked from 'package:stats':
## 
##     loadings
set.seed(123) # se fija una semilla para los generadores de numeros aleaorios en R
pcr.fit<-pcr(Sales~.,data=Carseats_new[train,],scale=T,validation="CV") #esta funci?n realiza la regresion por componentes
#pcr.fit<-pcr(Sales~.,data=Carseats_new[train,],scale=TRUE,validation="CV") #esta funci?n realiza la regresion por componentes
summary(pcr.fit)
## Data:    X dimension: 320 11 
##  Y dimension: 320 1
## Fit method: svdpc
## Number of components considered: 11
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV           2.827    2.733    2.706    2.616    2.630    2.631    2.640
## adjCV        2.827    2.731    2.706    2.615    2.629    2.632    2.631
##        7 comps  8 comps  9 comps  10 comps  11 comps
## CV       2.627    2.425    2.341     1.078     1.041
## adjCV    2.628    2.423    2.347     1.072     1.039
## 
## TRAINING: % variance explained
##        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
## X       17.517   33.191    46.44    56.91    66.10    74.54    82.82
## Sales    6.978    8.991    15.67    16.36    16.63    19.94    20.28
##        8 comps  9 comps  10 comps  11 comps
## X        90.11    94.09     97.58    100.00
## Sales    32.24    38.55     86.55     87.36
#Vemos que ocupamos practicamente todos los componentes
validationplot(pcr.fit)

pcr.fit<-pcr(Sales~.,data=Carseats_new,subset=train,scale=T,validation="CV") #esta funci?n realiza la regresion por componentes
summary(pcr.fit)
## Data:    X dimension: 320 11 
##  Y dimension: 320 1
## Fit method: svdpc
## Number of components considered: 11
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV           2.827    2.733    2.725    2.621    2.634    2.633    2.643
## adjCV        2.827    2.731    2.724    2.619    2.631    2.631    2.632
##        7 comps  8 comps  9 comps  10 comps  11 comps
## CV       2.606    2.419    2.345     1.077     1.040
## adjCV    2.611    2.416    2.353     1.074     1.038
## 
## TRAINING: % variance explained
##        1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
## X       17.517   33.191    46.44    56.91    66.10    74.54    82.82
## Sales    6.978    8.991    15.67    16.36    16.63    19.94    20.28
##        8 comps  9 comps  10 comps  11 comps
## X        90.11    94.09     97.58    100.00
## Sales    32.24    38.55     86.55     87.36
validationplot(pcr.fit,val.type="MSEP")

pcr.pred<-predict(pcr.fit,newdata = Carseats_new[-train,-1],ncomp =10)
erro_pcr<-mean((pcr.pred - Carseats_new$Sales[-train])^2) # se calcula el error de prueba
coef(pcr.fit)
## , , 11 comps
## 
##                        Sales
## CompPrice        1.421678289
## Income           0.435574293
## Advertising      0.818280207
## Population      -0.005471083
## Price           -2.278382091
## ShelveLocGood    1.993481567
## ShelveLocMedium  1.014605061
## Age             -0.778616707
## Education       -0.072352114
## UrbanYes         0.045327464
## USYes           -0.042326169
erro_pcr
## [1] 1.1638

Comparacion de errores MSE de prueba final

res_final<-data.frame(c("lasso_1se","lasso_min","ridge_1se","ridge_min","exhaust","PCR"),c(erro_lasso_1se,erro_lasso_lmin,erro_ridge_1se,erro_ridge_lmin,error_bestsubset,erro_pcr))
colnames(res_final)<-NULL
res_final

Podemos ver que PCR usando 10 componentes No supera a los modelos de Ridge Lasso o Exhaustivo, por lo cual ademas que es menos facil de interpretar su resultado en base a los componentes principales,

El mejor modelo si se quiere mas parsimonia seria el Lasso con un error estandard en lambda y si se quiere mayor precisión el Lasso con el lambda que minimiza

Dado que no se especifica el objetivo considerare que el mejor modelo es el lasso usando un lambda con un errror estandard porque es el mas parsimonioso y no le afecta mucho a la precision del MSE de prueba.

Coeficientes de modelo seleccionado

coef(lasso_1se)
## 12 x 1 sparse Matrix of class "dgCMatrix"
##                          s0
## (Intercept)      6.73364438
## CompPrice        0.07664656
## Income           0.01163099
## Advertising      0.10436310
## Population       .         
## Price           -0.08504838
## ShelveLocGood    4.36087400
## ShelveLocMedium  1.57627189
## Age             -0.04259526
## Education        .         
## UrbanYes         .         
## USYes            .

Error MSE de prueba

erro_lasso_1se
## [1] 1.083014

En este modelo Population, Education, Urban y US, tienen coeficiente cero por lo cual no los ocpuamos esos datos.

Updated: