Regularización: Lasso, Ridge y Elastic Net

library(dplyr)
library(tidyr)
library(broom)
library(ggplot2)
library(ISLR)
library(GGally)
library(modelr)
library(cowplot)
library(glmnet)
library(rlang)
library(purrr)
library(caret)
set.seed(1992)

Las técnicas de regularización son útiles para trabajar con conjuntos con gran cantidad de variables, las cuales pueden introducir variabilidad en las estimaciones de los parámetros. El problema que vamos a tratar de resolver es predecir el salario de un jugador de la NBA en base a ciertos predictores.

Conjunto de datos

Vamos a utilizar dos conjuntos de datos provenientes de Kaggle:

1. stats cuenta con las estadísticas de los jugadores (https://www.kaggle.com/koki25ando/salary)

2. salary cuenta con los salarios de los jugadores (https://www.kaggle.com/drgilermo/nba-players-stats)

Year: Season

Player: name

Pos: Position

Age: Age

Tm: Team

G: Games

GS: Games Started

MP: Minutes Played

PER: Player Efficiency Rating

TS%: True Shooting % : medida de eficiencia del jugador al lanzar

3PAr: 3-Point Attempt Rate

FTr:Free Throw Rate

ORB%: Offensive Rebound Percentage

DRB%: Defensive Rebound Percentage

TRB%: Total Rebound Percentage

AST%: Assist Percentage

STL%: Steal Percentage

BLK%: Block Percentage

TOV%: Turnover Percentage

USG%: Usage Percentage

OWS: Offensive Win Shares

DWS: Defensive Win Shares

WS: Win Shares

WS/48: Win Shares Per 48 Minutes

OBPM: Offensive Box Plus/Minus

DBPM: Defensive Box Plus/Minus

BPM: Box Plus/Minus : comparaciones contra el “jugador promedio” de la NBA

VORP: Value Over Replacement: Medida vinculada al BPM

FG: Field Goals

FGA: Field Goal Attempts

FG%: Field Goal Percentage

3P: 3-Point Field Goals

3PA: 3-Point Field Goal Attempts

3P%: 3-Point Field Goal Percentage

2P: 2-Point Field Goals

2PA: 2-Point Field Goal Attempts

2P%: 2-Point Field Goal Percentage

eFG%: Effective Field Goal Percentage

FT: Free Throws

FTA: Free Throw Attempts

FT%: Free Throw Percentage

ORB: Offensive Rebounds

DRB: Defensive Rebounds

TRB: Total Rebounds

AST: Assists

STL: Steals

BLK: Blocks

TOV: Turnovers

PF: Personal Fouls

PTS:Points

# Los datos de salario son para la temporada 2017-2018
salarios <- read.csv("NBA_season1718_salary.csv", stringsAsFactors = F) %>% rename(salary=season17_18)
# Filtramos las estadisticas para quedarnos con la temporada 2017
estadisticas <- read.csv("Seasons_Stats.csv", stringsAsFactors = F) %>% filter(Year>= 2017) %>% select(-c(blanl, blank2))
glimpse(salarios)

Observations: 573
Variables: 4
$ X      <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 3...
$ Player <chr> "Stephen Curry", "LeBron James", "Paul Millsap", "Gordon Hayward", "Blake Griffin", "Kyle Lowry", "Russell Westbrook", "Mike Con...
$ Tm     <chr> "GSW", "CLE", "DEN", "BOS", "DET", "TOR", "OKC", "MEM", "HOU", "TOR", "BOS", "OKC", "POR", "NOP", "MIA", "GSW", "WAS", "HOU", "P...
$ salary <dbl> 34682550, 33285709, 31269231, 29727900, 29512900, 28703704, 28530608, 28530608, 28299399, 27739975, 27734405, 26243760, 26153057...

colnames(estadisticas)

 [1] "X"      "Year"   "Player" "Pos"    "Age"    "Tm"     "G"      "GS"     "MP"     "PER"    "TS."    "X3PAr"  "FTr"    "ORB."   "DRB."   "TRB."  
[17] "AST."   "STL."   "BLK."   "TOV."   "USG."   "OWS"    "DWS"    "WS"     "WS.48"  "OBPM"   "DBPM"   "BPM"    "VORP"   "FG"     "FGA"    "FG."   
[33] "X3P"    "X3PA"   "X3P."   "X2P"    "X2PA"   "X2P."   "eFG."   "FT"     "FTA"    "FT."    "ORB"    "DRB"    "TRB"    "AST"    "STL"    "BLK"   
[49] "TOV"    "PF"     "PTS"   

Vemos que salarios contiene el nombre, equipo y salario en dolares de los jugadores.

Estadisticas es una tabla que contiene 51 variables, de las cuales 46 contienen informacion distintas estadisticas de la temporada 2017-2018.

Hacemos el join por jugador y equipo. Un jugador puede ser cambiado/vendido durante la temporada. Asi mantenemos a los jugadores con las estadisticas correspondientes al equipo con quien acordaron el salario.

nba <- inner_join(salarios, estadisticas, by=c('Player', 'Tm')) %>%
  select(-c(X.x, X.y, Year)) %>% 
  drop_na() # Eliminamos los NA

Analisis Exploratorio

Gráfico de la relacion entre la posicion y el salario

ggplot(nba, aes(Pos, salary, fill=Pos)) +
  geom_boxplot() +
  geom_text(aes(label=ifelse(salary>28700000,as.character(Player),'')),hjust=-0.1,vjust=0) +
  theme_bw() +
  labs(title= "Boxplots: salarios y posicion de juego", x="Posicion", y="Salario")

Correlagrama

ggcorr(nba, layout.exp = 2) + labs(title='Correlograma variables cuantitativas')

data in column(s) 'Player', 'Tm', 'Pos' are not numeric and were ignored

GGpairs (algunas variables)

Seleccionamos algunas variables y vemos sus relaciones usando ggpairs.

nba %>% select(salary, Age, PTS, GS, DRB, TRB, AST, BLK) %>% ggpairs() + theme_bw()

Modelo Lineal

Vamos a probar un modelo lineal que incluya todas las variables (excepto al jugador y equipo). Obtengan las estimaciones de los parametros junto a su p-valor e intervalo de confianza.

Coeficientes estimados

Vemos los coeficientes estimados y sus p-valores asociados

# Eliminamos jugador y equipo
nba = nba %>% select(-c(Player, Tm)) 
# Modelo lineal
modelo_lineal = nba %>% lm(formula = salary~., data = .)
#Coeficientes
lineal_coef= modelo_lineal %>% tidy(conf.int=TRUE)
ggplot(lineal_coef, aes(term, estimate))+
  geom_point()+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
  labs(title = "Coeficientes de la regresion lineal", x="", y="Estimacion e Int. Confianza") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90))

ggplot(lineal_coef, aes(reorder(term, -p.value), p.value, fill=p.value))+
  geom_bar(stat = 'identity', aes(fill=p.value))+
  labs(title = "P-valor de los regresores", x="", y="P-valor") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90)) + 
  scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )

¿Qué notamos aqui?

Hay ciertas coeficientes estimados que presentan una gran variabilidad pero la escala de las variables puede ocultarnos la verdadera variabilidad de los estimadores.

# Reescalamos las variables numericas
nba_scaled = nba %>% mutate_at(vars(-Pos), scale)
# Nuevo modelo lineal 
modelo_lineal_scal = nba_scaled %>% lm(formula = salary~., data = .)
lineal_coef_scal = modelo_lineal_scal %>% tidy(conf.int=TRUE)
ggplot(lineal_coef_scal, aes(term, estimate))+
  geom_point()+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high))+
  labs(title = "Coeficientes de la regresion lineal", subtitle="Variables escaladas", x="", y="Estimacion e Int. Confianza") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90))

ggplot(lineal_coef_scal, aes(reorder(term, -p.value), p.value, fill=p.value))+
  geom_bar(stat = 'identity', aes(fill=p.value))+
  labs(title = "P-valor de los regresores",subtitle="Variables escaladas", x="", y="P-valor") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90)) + 
  scale_fill_gradient2(high='firebrick', low = 'forestgreen', mid='yellow2',midpoint = 0.5 )

Evaluacion de los modelos

Obtemos la evaluacion de ambos modelos ¿Cómo esperan que sean los valores de diagnóstico para ambos modelos?

modelo_lineal = modelo_lineal %>% glance() %>% select(r.squared, adj.r.squared, p.value) 
modelo_lineal_scal = modelo_lineal_scal %>% glance() %>% select(r.squared, adj.r.squared, p.value)
bind_rows(modelo_lineal, modelo_lineal_scal) %>% mutate(modelo= c('lineal', 'lineal_escalado'))

Partición Train y Testing

Realizamos una partición entre dataset de entrenamiento y testeo usando la función resample_partition del paquete modelr

train_test <- nba %>% resample_partition(c(train=0.7,test=0.3))
nba <- train_test$train %>% as_tibble()
test <- train_test$test %>% as_tibble()

Regularizacion

La libreria glmnet es permite trabajar con modelos ridge, lasso y elastic net. La funcion que vamos a utilizar es glmnet. Es necesario que le pasemos un objeto matriz con los regresores y un vector con la variable a explicar (en este caso los salarios)

Con el parametro \(\alpha\) indicamos con que tipo de modelo deseamos trabajar:

• Ridge: \(\alpha=0\)

• Lasso: \(\alpha=1\)

• Elastic Net: \(0<\alpha<1\)

Lasso

En este caso vamos a trabajar con \(\alpha=1\).

1. ¿Cuál es la penalización que introduce el modelo Lasso?

2. ¿Cómo impacta esto en las variables?

# Vector con los salarios
nba_salary = nba$salary
# Matriz con los regresores
nba_mtx = model.matrix(salary~., data = nba)
# Modelo Lasso
lasso.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = F) # Que esta haciendo este parametro?
                 
lasso_coef = lasso.mod %>% tidy()
lasso_coef

Gráficos de analisis

El comando plot nos permite realizar dos graficos relevantes.

Grafico de coeficientes en funcion del lambda

plot(lasso.mod, 'lambda')

Grafico de coeficientes en funcion de la norma de penalizacion

plot(lasso.mod)

¿Qué muestra cada uno de estos graficos?

Podemos realizar los graficos para los valores de lambda en ggplot.

g1=lasso_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Lasso con Intercepto",  y="Coeficientes")
g2=lasso_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Lasso sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)

Veamos un poco mejor aquellas variables que sobreviven para mayores valores de lambda ¿Qué tienen en común todas estas variables?

lasso_coef %>% filter(term %in% c("G", "GS", "MP", "FGA", "X2PA", "TRB", "PF", "PTS")) %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line(size=1) + theme_bw() +
  labs(title="Lasso sin Intercepto", y="Coeficientes", subtitle= "\"Mejores\" variables") +
  scale_color_brewer(palette = 'Set1')

Vemos que las variables que “sobreviven” para mayores valores de lambda son las que están medidas con una escala mayor.

Estandarizacion en glmnet

Existen dos maneras de poder estandarizar las variables en glmnet.

1. Setear standardize = TRUE. Con esto se estandariza las regresoras, los coeficientes estimados estan en la escala original de la variable

2. Pasar los conjuntos de datos estandarizados.

Lasso estandarizado

# Modelo lasso
lasso.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = TRUE) # Estandarizamos
                 
lasso_coef = lasso.mod %>% tidy()
lasso_coef

Gráficos de analisis

El comando plot nos permite realizar dos graficos relevantes.

Grafico de coeficientes en funcion del lambda

plot(lasso.mod, 'lambda')

Grafico de coeficientes en funcion de la norma de penalizacion

plot(lasso.mod)

g1=lasso_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Lasso con Intercepto",  y="Coeficientes")
g2=lasso_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Lasso sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)

¿Podemos decidir cuál es el valor óptimo de lambda?

Elección lambda óptimo

Para elegir el valor óptimo de lambda, lo común es realizar cross-validation. La función cv.glmnet nos permite realizar esto de manera sencilla.

Al igual que para la función glmnet cuenta con los parametros:

• x: matriz de variables

• y: vector de la variable a predecir

• alpha: tipo de modelo

• standardize: flag logico para estandarizar las variables

Nuevo parametro

• type.measure: funcion de perdida/error que se va a utilizar en CV. Para los modelos de regularizacion el default es MSE

Salida Base

lasso_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=1, standardize = T)
lasso_cv

$lambda
 [1] 5823214.261 5305895.956 4834534.784 4405048.038 4013715.712 3657148.270 3332257.298 3036228.744 2766498.551 2520730.445 2296795.699 2092754.700
[13] 1906840.140 1737441.718 1583092.184 1442454.638 1314310.944 1197551.182 1091164.036  994228.030  905903.552  825425.578  752097.044  685282.816
[25]  624404.179  568933.833  518391.319  472338.863  430377.581  392144.022  357307.028  325564.857  296642.572  270289.664  246277.875  224399.227
[37]  204464.218  186300.181  169749.786  154669.683  140929.255  128409.488  117001.943  106607.813   97137.069   88507.679   80644.901   73480.630
[49]   66952.814   61004.910   55585.402   50647.348   46147.977   42048.318   38312.862   34909.253   31808.011   28982.275   26407.569   24061.594
[61]   21924.028   19976.358   18201.713   16584.723   15111.382   13768.929   12545.735   11431.206   10415.690    9490.388    8647.288    7879.087
[73]    7179.131    6541.357    5960.240    5430.749    4948.296    4508.703    4108.163    3743.205

$cvm
 [1] 6.363118e+13 5.926281e+13 5.433217e+13 5.000961e+13 4.628932e+13 4.319022e+13 4.048677e+13 3.815649e+13 3.616044e+13 3.447761e+13 3.311936e+13
[12] 3.202560e+13 3.113311e+13 3.039201e+13 2.980580e+13 2.925888e+13 2.857799e+13 2.795221e+13 2.745739e+13 2.703128e+13 2.669242e+13 2.646202e+13
[23] 2.627007e+13 2.611997e+13 2.600910e+13 2.591685e+13 2.585481e+13 2.585990e+13 2.593769e+13 2.603345e+13 2.610286e+13 2.616221e+13 2.623036e+13
[34] 2.630096e+13 2.635939e+13 2.649800e+13 2.666127e+13 2.682751e+13 2.702130e+13 2.718456e+13 2.729957e+13 2.740067e+13 2.748624e+13 2.753810e+13
[45] 2.755868e+13 2.759864e+13 2.765643e+13 2.774204e+13 2.787288e+13 2.798757e+13 2.815954e+13 2.838111e+13 2.857094e+13 2.879562e+13 2.903841e+13
[56] 2.929185e+13 2.959957e+13 2.987821e+13 3.014961e+13 3.039210e+13 3.060300e+13 3.079236e+13 3.093575e+13 3.105150e+13 3.117758e+13 3.131139e+13
[67] 3.145076e+13 3.159369e+13 3.174168e+13 3.187316e+13 3.200487e+13 3.211199e+13 3.223191e+13 3.239179e+13 3.254625e+13 3.268166e+13 3.277286e+13
[78] 3.294040e+13 3.304629e+13 3.314435e+13

$cvsd
 [1] 7.736798e+12 7.388669e+12 6.649727e+12 6.006497e+12 5.473897e+12 5.048968e+12 4.723632e+12 4.465990e+12 4.262311e+12 4.107387e+12 3.998177e+12
[12] 3.925714e+12 3.881921e+12 3.851532e+12 3.828443e+12 3.785963e+12 3.727056e+12 3.692428e+12 3.659408e+12 3.607617e+12 3.543249e+12 3.487758e+12
[23] 3.443244e+12 3.424927e+12 3.431368e+12 3.443013e+12 3.459304e+12 3.493783e+12 3.519312e+12 3.539710e+12 3.566708e+12 3.599479e+12 3.623850e+12
[34] 3.640486e+12 3.652095e+12 3.679458e+12 3.710503e+12 3.747061e+12 3.795228e+12 3.825423e+12 3.837981e+12 3.855719e+12 3.880706e+12 3.904739e+12
[45] 3.932997e+12 3.964577e+12 3.991326e+12 4.010399e+12 4.018100e+12 4.009198e+12 4.003415e+12 4.011381e+12 4.009916e+12 4.023713e+12 4.032237e+12
[56] 4.034433e+12 4.050743e+12 4.056274e+12 4.063082e+12 4.071753e+12 4.075389e+12 4.087869e+12 4.095455e+12 4.104105e+12 4.099421e+12 4.096852e+12
[67] 4.105106e+12 4.109080e+12 4.116089e+12 4.131878e+12 4.167937e+12 4.195578e+12 4.213773e+12 4.237172e+12 4.258212e+12 4.273418e+12 4.291434e+12
[78] 4.304215e+12 4.302191e+12 4.305480e+12

$cvup
 [1] 7.136798e+13 6.665148e+13 6.098189e+13 5.601611e+13 5.176322e+13 4.823919e+13 4.521040e+13 4.262248e+13 4.042275e+13 3.858499e+13 3.711754e+13
[12] 3.595131e+13 3.501503e+13 3.424354e+13 3.363424e+13 3.304484e+13 3.230504e+13 3.164464e+13 3.111679e+13 3.063889e+13 3.023567e+13 2.994977e+13
[23] 2.971332e+13 2.954490e+13 2.944046e+13 2.935986e+13 2.931412e+13 2.935368e+13 2.945700e+13 2.957316e+13 2.966957e+13 2.976169e+13 2.985421e+13
[34] 2.994145e+13 3.001148e+13 3.017746e+13 3.037178e+13 3.057457e+13 3.081652e+13 3.100998e+13 3.113755e+13 3.125639e+13 3.136694e+13 3.144283e+13
[45] 3.149167e+13 3.156322e+13 3.164776e+13 3.175244e+13 3.189098e+13 3.199676e+13 3.216296e+13 3.239249e+13 3.258086e+13 3.281933e+13 3.307065e+13
[56] 3.332628e+13 3.365032e+13 3.393448e+13 3.421269e+13 3.446385e+13 3.467839e+13 3.488023e+13 3.503120e+13 3.515561e+13 3.527700e+13 3.540824e+13
[67] 3.555586e+13 3.570277e+13 3.585777e+13 3.600504e+13 3.617281e+13 3.630757e+13 3.644568e+13 3.662896e+13 3.680446e+13 3.695507e+13 3.706429e+13
[78] 3.724462e+13 3.734848e+13 3.744983e+13

$cvlo
 [1] 5.589438e+13 5.187414e+13 4.768244e+13 4.400312e+13 4.081543e+13 3.814125e+13 3.576313e+13 3.369050e+13 3.189813e+13 3.037022e+13 2.912118e+13
[12] 2.809988e+13 2.725119e+13 2.654048e+13 2.597736e+13 2.547291e+13 2.485093e+13 2.425979e+13 2.379798e+13 2.342366e+13 2.314917e+13 2.297426e+13
[23] 2.282683e+13 2.269504e+13 2.257773e+13 2.247384e+13 2.239551e+13 2.236612e+13 2.241837e+13 2.249374e+13 2.253615e+13 2.256273e+13 2.260651e+13
[34] 2.266048e+13 2.270729e+13 2.281855e+13 2.295077e+13 2.308045e+13 2.322607e+13 2.335914e+13 2.346159e+13 2.354495e+13 2.360553e+13 2.363336e+13
[45] 2.362568e+13 2.363407e+13 2.366510e+13 2.373164e+13 2.385478e+13 2.397837e+13 2.415613e+13 2.436972e+13 2.456102e+13 2.477190e+13 2.500618e+13
[56] 2.525742e+13 2.554883e+13 2.582194e+13 2.608653e+13 2.632035e+13 2.652761e+13 2.670449e+13 2.684029e+13 2.694740e+13 2.707816e+13 2.721453e+13
[67] 2.734565e+13 2.748461e+13 2.762559e+13 2.774129e+13 2.783694e+13 2.791641e+13 2.801814e+13 2.815462e+13 2.828804e+13 2.840824e+13 2.848142e+13
[78] 2.863619e+13 2.874410e+13 2.883887e+13

$nzero
 s0  s1  s2  s3  s4  s5  s6  s7  s8  s9 s10 s11 s12 s13 s14 s15 s16 s17 s18 s19 s20 s21 s22 s23 s24 s25 s26 s27 s28 s29 s30 s31 s32 s33 s34 s35 s36 
  0   1   2   2   2   3   3   3   3   3   3   4   4   4   5   7   7   7   7   6   6   7   7   7   7   7   7   7   7   8   9  10  10  11  11  13  15 
s37 s38 s39 s40 s41 s42 s43 s44 s45 s46 s47 s48 s49 s50 s51 s52 s53 s54 s55 s56 s57 s58 s59 s60 s61 s62 s63 s64 s65 s66 s67 s68 s69 s70 s71 s72 s73 
 16  19  19  22  23  22  23  24  22  22  22  22  22  24  27  28  30  31  33  34  34  36  37  37  38  38  38  38  38  42  42  42  43  42  39  40  39 
s74 s75 s76 s77 s78 s79 
 39  40  41  39  39  41 

$name
                 mse 
"Mean-Squared Error" 

$glmnet.fit

Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 1, standardize = T) 

       Df    %Dev    Lambda
  [1,]  0 0.00000 5823000.0
  [2,]  1 0.09126 5306000.0
  [3,]  2 0.16950 4835000.0
  [4,]  2 0.23850 4405000.0
  [5,]  2 0.29570 4014000.0
  [6,]  3 0.34420 3657000.0
  [7,]  3 0.38710 3332000.0
  [8,]  3 0.42260 3036000.0
  [9,]  3 0.45210 2766000.0
 [10,]  3 0.47670 2521000.0
 [11,]  3 0.49700 2297000.0
 [12,]  4 0.51430 2093000.0
 [13,]  4 0.52960 1907000.0
 [14,]  4 0.54230 1737000.0
 [15,]  5 0.55540 1583000.0
 [16,]  7 0.56950 1442000.0
 [17,]  7 0.58310 1314000.0
 [18,]  7 0.59440 1198000.0
 [19,]  7 0.60370 1091000.0
 [20,]  6 0.61000  994200.0
 [21,]  6 0.61510  905900.0
 [22,]  7 0.61950  825400.0
 [23,]  7 0.62460  752100.0
 [24,]  7 0.62880  685300.0
 [25,]  7 0.63230  624400.0
 [26,]  7 0.63520  568900.0
 [27,]  7 0.63760  518400.0
 [28,]  7 0.63960  472300.0
 [29,]  7 0.64130  430400.0
 [30,]  8 0.64370  392100.0
 [31,]  9 0.64810  357300.0
 [32,] 10 0.65240  325600.0
 [33,] 10 0.65630  296600.0
 [34,] 11 0.65980  270300.0
 [35,] 11 0.66290  246300.0
 [36,] 13 0.66590  224400.0
 [37,] 15 0.66860  204500.0
 [38,] 16 0.67140  186300.0
 [39,] 19 0.67430  169700.0
 [40,] 19 0.67730  154700.0
 [41,] 22 0.67980  140900.0
 [42,] 23 0.68400  128400.0
 [43,] 22 0.68740  117000.0
 [44,] 23 0.69020  106600.0
 [45,] 24 0.69250   97140.0
 [46,] 22 0.69460   88510.0
 [47,] 22 0.69620   80640.0
 [48,] 22 0.69760   73480.0
 [49,] 22 0.69870   66950.0
 [50,] 22 0.69970   61000.0
 [51,] 24 0.70070   55590.0
 [52,] 27 0.70220   50650.0
 [53,] 28 0.70410   46150.0
 [54,] 30 0.70570   42050.0
 [55,] 31 0.70750   38310.0
 [56,] 33 0.70910   34910.0
 [57,] 34 0.71070   31810.0
 [58,] 34 0.71230   28980.0
 [59,] 36 0.71380   26410.0
 [60,] 37 0.71530   24060.0
 [61,] 37 0.71660   21920.0
 [62,] 38 0.71770   19980.0
 [63,] 38 0.71890   18200.0
 [64,] 38 0.71990   16580.0
 [65,] 38 0.72070   15110.0
 [66,] 38 0.72130   13770.0
 [67,] 42 0.72190   12550.0
 [68,] 42 0.72240   11430.0
 [69,] 42 0.72330   10420.0
 [70,] 43 0.72450    9490.0
 [71,] 42 0.72620    8647.0
 [72,] 39 0.72760    7879.0
 [73,] 40 0.72760    7179.0
 [74,] 39 0.72810    6541.0
 [75,] 39 0.72840    5960.0
 [76,] 40 0.72860    5431.0
 [77,] 41 0.72890    4948.0
 [78,] 39 0.72960    4509.0
 [79,] 39 0.73010    4108.0
 [80,] 41 0.73060    3743.0
 [81,] 43 0.73120    3411.0
 [82,] 43 0.73180    3108.0
 [83,] 42 0.73240    2832.0
 [84,] 42 0.73290    2580.0
 [85,] 42 0.73430    2351.0
 [86,] 42 0.73560    2142.0
 [87,] 43 0.73650    1952.0
 [88,] 44 0.73730    1778.0
 [89,] 44 0.73800    1620.0
 [90,] 45 0.73860    1476.0
 [91,] 45 0.73910    1345.0
 [92,] 45 0.73960    1226.0
 [93,] 45 0.74000    1117.0
 [94,] 45 0.74030    1018.0
 [95,] 45 0.74060     927.2
 [96,] 45 0.74090     844.8
 [97,] 45 0.74110     769.8
 [98,] 45 0.74130     701.4
 [99,] 46 0.74140     639.1
[100,] 48 0.74140     582.3

$lambda.min
[1] 518391.3

$lambda.1se
[1] 1442455

attr(,"class")
[1] "cv.glmnet"

Brinda muchisima informacion:

• lambda

• cvm (Cross-validation mean): es la media del MSE (error)

• cvsd (Cross-validation Standard Error): desvio estandar del MSE (error)

• cvup y cvlo: Limite superior e inferior

• nzero: Coeficientes distintos de cero

• lambda.min: lambda para el cual el MSE (error) es minimo

En glm.fit tenemos la cantidad de variables, el valor de lambda y el porcentaje de deviance explicada por el modelo

Grafico Base

plot(lasso_cv)

El gráfico nos muestra la media del MSE con su limite superior e inferior y la cantidad de varaibles que sobreviven para cada valor de lambda.

Modelr

lasso_cv %>% tidy()

lasso_cv %>% glance()

Seleccionamos el lambda optimo para crear el modelo final

lasso_lambda_opt = lasso_cv$lambda.min
lasso_opt = glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = lasso_lambda_opt)
# Salida estandar
lasso_opt

Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 1, lambda = lasso_lambda_opt,      standardize = TRUE) 

     Df   %Dev Lambda
[1,]  7 0.6376 518400

# Tidy
lasso_opt %>% tidy()

# Glance
lasso_opt %>% glance()

Ridge

En este caso vamos a trabajar con \(\alpha=0\). Vamos a replicar lo que ya realizamos para Lasso.

1. ¿Cuál es la penalización que introduce el modelo Ridge?

2. ¿Cómo impacta esto en las variables?

ridge.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0, # Indicador del tipo de regularizacion
                 standardize = TRUE)
                 
ridge_coef= ridge.mod %>% tidy()
ridge_coef 

¿Qué ven de distinto en los coeficientes estimados del modelo respecto a Lasso?

Gráficos de analisis

Grafico de coeficientes en funcion del lambda

plot(ridge.mod, 'lambda')

¿Qué ven de distinto en este gráfico respecto al que obtuvimos con la regresión Lasso?

Grafico de coeficientes en funcion de la norma de penalizacion

plot(ridge.mod)

g1=ridge_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Ridge con Intercepto",  y="Coeficientes")
g2=ridge_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Ridge sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)

Elección lambda óptimo

ridge_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=0, standardize = T)

Grafico Base

plot(ridge_cv)

Seleccionamos el lambda óptimo para crear el modelo final

ridge_lambda_opt = ridge_cv$lambda.min
ridge_opt = glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = ridge_lambda_opt)
# Salida estandar
ridge_opt

Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 0, lambda = ridge_lambda_opt,      standardize = TRUE) 

     Df   %Dev  Lambda
[1,] 50 0.6392 6541000

# Tidy
ridge_opt %>% tidy()

Elastic Net

El modelo Elastic Net incorpora los dos tipos de penalización: Lasso (Norma L1) y Ridge (Norma L2). El parámetro \(\alpha\) regula la importancia de cada penalización, cuanto más cerca de cero será más importante la penalización del tipo Ridge y más cerca de 1, la tipo Lasso.

En este caso vamos a trabajar con \(\alpha=0.5\). Vamos a replicar lo que ya realizamos para Lasso y Ridge

elastic.mod=glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0.5, # Indicador del tipo de regularizacion
                 standardize = TRUE)
                 
elastic_coef= elastic.mod %>% tidy()
elastic_coef 

¿Qué ven de distinto en los coeficientes estimados del modelo respecto a Lasso y Ridge?

Gráficos de analisis

Grafico de coeficientes en funcion del lambda

plot(elastic.mod, 'lambda')

¿Qué ven en este gráfico de distinto a los dos anteriores?

g1=elastic_coef  %>% ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Elastic Net con Intercepto",  y="Coeficientes")
g2=elastic_coef %>% filter(term!='(Intercept)') %>% 
  ggplot(., aes(log(lambda), estimate, group=term, color=term)) + geom_line() + theme_bw()  + theme(legend.position = 'none') +
  labs(title="Elastic Net sin Intercepto", y="Coeficientes")
plot_grid(g1,g2)

Elección lambda óptimo

elastic_cv=cv.glmnet(x=nba_mtx,y=nba_salary,alpha=0.5, standardize = T)

Grafico Base

Presten especial atención al eje superior ¿Qué está sucediendo?

plot(elastic_cv)

Seleccionamos el lambda optimo para crear el modelo final

elastic_lambda_opt = elastic_cv$lambda.min
elastic_opt = glmnet(x=nba_mtx, # Matriz de regresores
                 y=nba_salary, #Vector de la variable a predecir
                 alpha=0.5, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = elastic_lambda_opt)
# Salida estandar
elastic_opt

Call:  glmnet(x = nba_mtx, y = nba_salary, alpha = 0.5, lambda = elastic_lambda_opt,      standardize = TRUE) 

     Df   %Dev Lambda
[1,] 25 0.6949 161300

# Tidy
elastic_opt %>% tidy()

Breve comparacion entre modelos

ridge_dev = ridge_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
  ggplot(., aes(log(lambda), dev.ratio)) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = log(ridge_lambda_opt), color='steelblue', size=1.5) +
  labs(title='Ridge: Deviance') +
  theme_bw() 
lasso_dev = lasso_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
  ggplot(., aes(log(lambda), dev.ratio)) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = log(lasso_lambda_opt), color='firebrick', size=1.5) +
  labs(title='Lasso: Deviance') +
  theme_bw()
elastic_dev = elastic_coef %>% select(lambda, dev.ratio) %>% distinct() %>%
  ggplot(., aes(log(lambda), dev.ratio)) +
  geom_point() +
  geom_line() +
  geom_vline(xintercept = log(elastic_lambda_opt), color='forestgreen', size=1.5) +
  labs(title='Elastic Net: Deviance') +
  theme_bw()
plot_grid(ridge_dev, lasso_dev, elastic_dev)

Testing

Con los modelos optimos que encontramos pueden probar cual es el MSE y RMSE en los datasets de training y testing, para decidir cual es el modelo que minimiza el error en las predicciones.