Lasso - Ridge y Elastic Net

Ejemplo

Se bajo el contenido de 816 sitios web de argentina
El mismo se proceso para generar una matriz de documentos y terminos
Estos se encuentran clasificados en categorias Shopping y Sports
Se busca generar un modelo para clasificar nuevas instancias en base al texto de los sitios
Al ser datos con muchos ceros (el 97% de la matriz), tiene sentido usar metodos penalizados.

Leo el archivo y Cargo paquetes

# Pongo un seed.
set.seed(1)

# Leo el archivo
d <- read.table('web_content.txt')

## Dimension (Matriz dispersa)
dim(d)

## [1]  816 2078

# Que se busca predecir?
# Si la pagina web pertenece a la categoria sports o shopping en base a su contenido (texto en la pagina)
table(d$clase)

## 
## shopping   sports 
##      497      319

## Nombres de las variables
names(d)[sample(2000, 100)]

##   [1] "definición"    "etapa"         "logrado"       "siet"         
##   [5] "conferencia"   "seleccionado"  "tandil"        "newel"        
##   [9] "miembro"       "armstrong"     "consecuencia"  "client"       
##  [13] "olímpico"      "ezeiza"        "posición"      "inform"       
##  [17] "pasa"          "unidad"        "existen"       "predio"       
##  [21] "suel"          "constitución"  "motor"         "cambiar"      
##  [25] "defensiva"     "extremo"       "adecuado"      "existencia"   
##  [29] "ríos"          "empezó"        "horario"       "manten"       
##  [33] "importancia"   "comenzó"       "quilm"         "negocio"      
##  [37] "princip"       "bosnia"        "participaron"  "florencio"    
##  [41] "pura"          "miramar"       "posteriorment" "lanú"         
##  [45] "javier"        "prensa"        "alcanzar"      "hicieron"     
##  [49] "pasión"        "ocho"          "herzegovina"   "región"       
##  [53] "garantía"      "cualidad"      "así"           "bella"        
##  [57] "disputa"       "intento"       "mostrar"       "fiesta"       
##  [61] "sarmiento"     "descuento"     "grecia"        "edición"      
##  [65] "mina"          "dado"          "hermosa"       "pinamar"      
##  [69] "azul"          "representant"  "elegido"       "queremo"      
##  [73] "empató"        "trabajar"      "hago"          "rivadavia"    
##  [77] "reconquista"   "europeo"       "plantel"       "sri"          
##  [81] "futbol"        "olx"           "facilidad"     "división"     
##  [85] "perder"        "compra"        "oficina"       "búsqueda"     
##  [89] "crítica"       "capitán"       "costado"       "arbitro"      
##  [93] "médico"        "recurso"       "pintura"       "poner"        
##  [97] "gimnasia"      "federación"    "premio"        "loma"

## top
d[1:5, 1:10]

##      abajo abierta  abierto abril abrió    abrir absoluta acá    acaba
## 1 0.004236       0 0.000000     0     0 0.004698        0   0 0.005393
## 2 0.000000       0 0.000000     0     0 0.000000        0   0 0.000000
## 3 0.000000       0 0.009674     0     0 0.000000        0   0 0.000000
## 4 0.000000       0 0.000000     0     0 0.000000        0   0 0.000000
## 5 0.000000       0 0.000000     0     0 0.000000        0   0 0.000000
##   academia
## 1        0
## 2        0
## 3        0
## 4        0
## 5        0

### Instalar los paquetes que digan FALSE
# En paralelo no anda en windows
library(caret)

## Loading required package: lattice
## Loading required package: ggplot2

library(glmnet)

## Loading required package: Matrix
## Loaded glmnet 1.9-8

library(doParallel)

## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel

library(corrplot)

# Para procesar en paralelo
registerDoParallel(cores=4)

# cor <- cor(d[,names(d)[!names(d) %in% c('clase')]])
# corrplot(cor)

Divido en Training y Testing

# Nombres de variables
predictors <- names(d)[!names(d) %in% "clase"]

# ?createDataPartition
set.seed(123)
inTrainingSet <- createDataPartition(d$clase, p = 0.7, list = FALSE)

train <- d[inTrainingSet, ]
test <- d[-inTrainingSet, ]

x = train[,predictors]
y = train$clase

# glmnet pide una matrix como input, uso as.matrix
# Si hubiera variables categoricas, habria que usar model.matrix()
x = as.matrix(x)

Ejemplo del path usando regresion lasso

set.seed(666)
fit_lasso = glmnet(x, y, family = "binomial", alpha = 1, standardize=TRUE)

# Coeficientes y Devianza
plot(fit_lasso, xvar = "dev", label = TRUE)

plot of chunk unnamed-chunk-4

# Grafico del "camino" (path) de soluciones
# El lambda mas chico, tiene un log(lambda) mas cercano a 0
# Mientras menor es log(lambda), menor es la penalizacion
plot(fit_lasso,  xvar = "lambda", label = TRUE)

plot of chunk unnamed-chunk-4

# Se pueden poner limites inferiores y superiores a los coeficientes
fit_lasso_restringido <- glmnet(x, y, family = "binomial", alpha = 1,
                                standardize=TRUE, lower=-40,upper=40)

plot(fit_lasso_restringido)

plot of chunk unnamed-chunk-4

plot(fit_lasso_restringido, xvar = "dev", label = TRUE)

plot of chunk unnamed-chunk-4

# Lambdas - Devianza
# fit_res <- as.data.frame(cbind(fit_lasso$df, fit_lasso$dev, fit_lasso$lambda))
# names(fit_res) <- c("Q Variables (DF)", "Devianza", "Lambda")
# head(fit_res, 100)

# Predicciones usando todos los lambdas
pred_mat <- predict(fit_lasso_restringido, as.matrix(test[,predictors]), type='class')
accuracy <- function(mat_confusion) sum(diag(mat_confusion)) / sum(mat_confusion)

eval_modelos <- function(x){
  # toma la matriz de confusion como input y calcula metricas
  tbl <- table(x, test$clase)
  if(dim(tbl)[[1]] == 1){
    list(tbl=tbl, acc=0, recall=0, precision=0, specificity=0, F1=0)
  } else {
    accuracy <- sum(diag(tbl)) / sum(tbl)
    sensitividad <- tbl[1,1] / sum(tbl[,1])
    especificidad <- tbl[2,2] / sum(tbl[,2])
    list(tbl=tbl, accuracy=accuracy, sensitividad=sensitividad, especificidad=especificidad)
  } 
}

res_fit_lasso <- apply(pred_mat, 2, eval_modelos)
metricas <- c('sensitividad', 'especificidad', 'accuracy')
met_lasso <- lapply(metricas, function(x) {
  as.vector(unlist(sapply(res_fit_lasso, "[[", x)))
})

met_lasso <- as.data.frame(do.call('cbind', met_lasso))
names(met_lasso) <- metricas
# met_lasso$lambda_fact <- as.factor(row.names(met_lasso))
met_lasso$lambda <- fit_lasso$lambda[-c(1:9)]
melt_met_lasso <- reshape2::melt(met_lasso, id.vars=c('lambda'))

# Grafico
ggplot(melt_met_lasso, aes(x=lambda, y=value, color=variable, group=variable)) +
  geom_line(size=1.2) + scale_x_continuous(breaks=seq(0, 1, 0.03)) +
  labs(title = 'Metricas - Regresion Lasso')

plot of chunk unnamed-chunk-4

Comparación entre lasso, ridge y elastic-net - Validacion cruzada

Elegimos el parametro de regularización por validación cruzada usando 3 folds
Comparamos el AUC de los tres modelos y otras metricas

Regresion Lasso

set.seed(666)
cvfit_lasso = cv.glmnet(x, y, family = "binomial", alpha = 1, nfold=3,
                        parallel = TRUE, standardize=TRUE, type.measur='auc')
plot(cvfit_lasso)

plot of chunk unnamed-chunk-5

Resultados Lasso

pred_lasso = predict(cvfit_lasso, newx = as.matrix(test[,predictors]), 
               s = "lambda.min", type = "class")
confusionMatrix(pred_lasso, test$clase)$table

##           Reference
## Prediction shopping sports
##   shopping      131      4
##   sports         18     91

Ridge - Incluye a todas las variables en la estimación

set.seed(666)
cvfit_ridge = cv.glmnet(x, y, family = "binomial", alpha = 0, nfold=3,
                        parallel=TRUE, standardize=TRUE, type.measur='auc')
plot(cvfit_ridge)

plot of chunk unnamed-chunk-7

Resultados Ridge

pred_ridge = predict(cvfit_ridge, newx = as.matrix(test[,predictors]), 
               s = "lambda.min", type = "class")
confusionMatrix(pred_ridge, test$clase)$table

##           Reference
## Prediction shopping sports
##   shopping      139      4
##   sports         10     91

Elastic - Net

Alpha 0.2 - más cerca de ridge

set.seed(666)
cvfit_elastic_02 = cv.glmnet(x, y, family = "binomial", alpha = 0.2, nfold=3, 
                        parallel=TRUE, standardize=TRUE, type.measur='auc')
plot(cvfit_elastic_02)

plot of chunk unnamed-chunk-9

Resultados Elastic-Net

pred_enet_02 = predict(cvfit_elastic_02, newx = as.matrix(test[,predictors]), s = "lambda.min", type = "class")
confusionMatrix(pred_enet_02, test$clase)$table

##           Reference
## Prediction shopping sports
##   shopping      131      1
##   sports         18     94

Elastic - Net

Alpha 0.8 - más cerca de lasso que ridge

set.seed(666)
cvfit_elastic_08 = cv.glmnet(x, y, family = "binomial", alpha = 0.8, nfold=3, 
                        parallel=TRUE, standardize=TRUE, type.measur='auc')
plot(cvfit_elastic_08)

plot of chunk unnamed-chunk-11

Resultados Elastic-Net

pred_enet_08 = predict(cvfit_elastic_08, newx = as.matrix(test[,predictors]), s = "lambda.min", type = "class")
confusionMatrix(pred_enet_08, test$clase)$table

##           Reference
## Prediction shopping sports
##   shopping      132      4
##   sports         17     91

Comparación de modelos

train_alpha <- function(alpha){
  set.seed(666)
  fit = cv.glmnet(x, y, family = "binomial", alpha = alpha, nfold=3, 
                  parallel=TRUE, standardize=TRUE, type.measur='auc')
  pred = predict(fit, newx = as.matrix(test[,predictors]), 
                 s = "lambda.min", type = "class")
  list(fit=fit, pred=pred)
}

# Defino una secuencia de alphas para modelar
alphas <- seq(0, 1, 0.01)

# Aplico la funcion train_alpha() al vector de alphas
enets <- lapply(alphas, function(a) train_alpha(a))
fits <- lapply(enets, '[[', "fit")
preds <- lapply(enets, '[[', "pred")

# Obtengo AUC y cantidad de variables de cada modelo.
AUC <- sapply(fits, function(x) x$cvm[x$lambda == x$lambda.min])
Q_VARIABLES <- sapply(fits, function(x) sum(coef(x, s='lambda.min') != 0))

# Métricas
resultados <- lapply(preds, function(x){
  tbl <- table(x, test$clase)
  accuracy <- accuracy(tbl)
  sensitividad <- tbl[1,1] / sum(tbl[,1])
  especificidad <- tbl[2,2] / sum(tbl[,2])
  F1 <- (2 * (especificidad * sensitividad) ) / (especificidad + sensitividad)
  list(tbl=tbl, accuracy=accuracy, sensitividad=sensitividad, especificidad=especificidad, F1=F1)
})

# Resultados de todos los modelos
modelos <- data.frame(
  modelo=alphas, 
  accuracy=sapply(resultados, "[[", 'accuracy'),
  sensitividad=sapply(resultados, "[[", 'sensitividad'),
  especificidad = sapply(resultados, "[[", 'especificidad'),
  F1=sapply(resultados, "[[", 'F1'), AUC, Q_VARIABLES
)

melt_modelos <- reshape2::melt(modelos, id.vars="modelo")

# AUC - Cantidad de Variables
ggplot(melt_modelos[melt_modelos$variable %in% c('AUC', 'Q_VARIABLES'),], 
       aes(x=modelo, y=value, color=variable, group=variable)) +
  geom_line(size=1.2) + 
  facet_grid(variable ~ ., scales="free") +
  scale_x_continuous(breaks=seq(0, 1, 0.05)) + 
  labs(list(title = 'Elastic Net - AUC - Cantidad de Variables', 
            x="Alpha", y="AUC - Cantidad de Variables"))

# Especificidad - Sensitividad
ggplot(melt_modelos[melt_modelos$variable %in% c('especificidad', 'sensitividad'),], 
       aes(x=modelo, y=value, color=variable, group=variable)) +
  geom_line(size=1.2) +
  scale_x_continuous(breaks=seq(0, 1, 0.1)) + 
  labs(list(title = 'Elastic Net - Especificidad - Sensitividad', 
          x="Alpha", y="Metricas"))

plot of chunk unnamed-chunk-13

# Todo junto
ggplot(melt_modelos, aes(x=modelo, y=value, color=variable, group=variable)) +
  geom_line(size=1.2) + facet_wrap(~ variable, scales="free")

plot of chunk unnamed-chunk-13

head(modelos, 10)

##    modelo accuracy sensitividad especificidad     F1    AUC Q_VARIABLES
## 1    0.00   0.9426       0.9329        0.9579 0.9452 0.9804        2078
## 2    0.01   0.8607       0.7785        0.9895 0.8714 0.9836         513
## 3    0.02   0.9262       0.8859        0.9895 0.9348 0.9836         478
## 4    0.03   0.9098       0.8591        0.9895 0.9197 0.9836         356
## 5    0.04   0.9467       0.9262        0.9789 0.9518 0.9835         375
## 6    0.05   0.9467       0.9262        0.9789 0.9518 0.9836         369
## 7    0.06   0.9385       0.9128        0.9789 0.9447 0.9838         308
## 8    0.07   0.9467       0.9262        0.9789 0.9518 0.9836         307
## 9    0.08   0.9385       0.9060        0.9895 0.9459 0.9836         258
## 10   0.09   0.9344       0.8993        0.9895 0.9423 0.9837         244

tail(modelos, 10)

##     modelo accuracy sensitividad especificidad     F1    AUC Q_VARIABLES
## 92    0.91   0.9180       0.8859        0.9684 0.9253 0.9755          48
## 93    0.92   0.9180       0.8859        0.9684 0.9253 0.9755          48
## 94    0.93   0.9180       0.8859        0.9684 0.9253 0.9760          46
## 95    0.94   0.9139       0.8792        0.9684 0.9217 0.9759          46
## 96    0.95   0.9139       0.8792        0.9684 0.9217 0.9758          46
## 97    0.96   0.9139       0.8792        0.9684 0.9217 0.9757          46
## 98    0.97   0.9098       0.8792        0.9579 0.9169 0.9757          46
## 99    0.98   0.9098       0.8792        0.9579 0.9169 0.9755          44
## 100   0.99   0.9098       0.8792        0.9579 0.9169 0.9755          44
## 101   1.00   0.9098       0.8792        0.9579 0.9169 0.9755          43

Lasso - Ridge y Elastic Net

EEA

Ejemplo

Leo el archivo y Cargo paquetes

Divido en Training y Testing

Ejemplo del path usando regresion lasso

Comparación entre lasso, ridge y elastic-net - Validacion cruzada

Regresion Lasso

Resultados Lasso

Ridge - Incluye a todas las variables en la estimación

Resultados Ridge

Elastic - Net

Alpha 0.2 - más cerca de ridge

Resultados Elastic-Net

Elastic - Net

Alpha 0.8 - más cerca de lasso que ridge

Resultados Elastic-Net

Comparación de modelos