Publicacion_lasso_ridge
Regresion Ridge VS Lasso
Adrián Alejandro Rodríguez Villarreal
22 de diciembre de 2018
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-16data("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.8691000summ_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.000000summ_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] 7which.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.048199n_predictores<-which.min(val.errors)
error_bestsubset<-min(val.errors)
cbind(n_predictores,error_bestsubset)## n_predictores error_bestsubset
## [1,] 10 1.046329Si 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.08633875Ahora 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.09135486lasso_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.01878724lasso_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.1488857ridge_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.0191305325cv.out_ridge$lambda.1se## [1] 0.2370681ridge_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.0653223192erro_lasso_1se<-mean( (predict(lasso_1se,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_lasso_1se## [1] 1.083014erro_lasso_lmin<-mean( (predict(lasso_lmin,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_lasso_lmin## [1] 1.041594erro_ridge_1se<-mean( (predict(ridge_1se,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_ridge_1se## [1] 1.156963erro_ridge_lmin<-mean( (predict(ridge_lmin,newx=x[-train,])-Carseats$Sales[-train])^2 )
erro_ridge_lmin## [1] 1.083586plot(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.08633875Otra 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':
##
## loadingsset.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.36validationplot(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.042326169erro_pcr## [1] 1.1638Comparacion 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_finalPodemos 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.083014En este modelo Population, Education, Urban y US, tienen coeficiente cero por lo cual no los ocpuamos esos datos.