Capítulo 8 Shrinkage Methods

Cuando el objetivo es encontrar el mejor modelo, es decir, aquel con menor error de predicción en la muestra de test existen distintas alternativas de selección de modelos. La idea de shrinkage es ajustar un modelo que contenga todos los \(p\) predictores usando una técnica que restringe o regulariza las estimaciones de coeficientes hacia cero. Puede que no sea obvio por qué tal restricción debería mejorar el ajuste, pero reducir los valores estimados de los coeficientes puede reducir significativamente su varianza. Las dos técnicas más conocidas para reducir los coeficientes de regresión hacia cero son la regresión ridge y lasso.

La estimación por mínimos cuadrados del Capítulo 7 son “equivalentes en escala”: multiplicar \(X_j\) por una constante \(c\) simplemente modifica el valor del coeficiente estimado por un factor de \(\frac{1}{c}\). En otras palabras, independientemente de cómo se escale el j-ésimo predictor, \(X_j \hat{\beta}_j\) seguirá siendo el mismo. Por el contrario, los coeficientes de regresión ridge y lasso pueden cambiar sustancialmente al multiplicar un predictor por una constante (notar que, en general, las variables tienen distinta escala de medición). Por lo tanto, es mejor estandarizar los predictores, por ejemplo, dividir cada variable por su desvío estándar.

8.1 LASSO

Least Absolute Shrinkage and Selection Operator. Dado:

\(Y = \beta_0 + \beta_1 X_1 +...+ \beta_p X_p + u\)

Suponer que \(M_k\) con \(k = 1,...,K\) es una serie modelos.

Elección de modelos: buscar en \(M_k\) el mejor modelo para predecir fuera de la muestra. La primera alternativa podría ser realizar una búsqueda exhaustiva. En la ocasiones, puede ser impracticable dado que con \(p\) predictores, se pueden construir \(2^p\) modelos (con \(p = 15\) hay \(32.768\) modelos).

Lasso: una manera formal y algorítmica de realizar esa tarea. Para \(\lambda \ge 0\) dado, se considera la siguiente función objetivo (a minimizar):

\[\begin{equation} \tag{8.1} R_l(\beta) = \sum_{i=1}^{n}(y_i - x_i^{'} \beta)^2 + \lambda \sum_{j=1}^{p}|\beta_j| \end{equation}\]

donde \(x_i\) contiene un intercepto.

\(\lambda = 0\) es el estimador de \(MCO\)
\(\lambda = \infty\) el segundo término \(\to \infty\) entonces \(\beta_1,...\beta_p = 0\)
Si \(0 < \lambda < \infty\) Lasso hace algo intermedio entre \(MCO\) y \(0\)

Intuitivamente:

\(\sum_{i=1}^{n}(y_i - x_i^{'} \beta)^2\) penaliza la falta de ajuste
\(\sum_{j=1}^{p}|\beta_j|\) penaliza complejidad

En determinados casos, la ventaja de utilizar lasso es que elige automáticamente qué variables entran en el modelo \((\beta_j \neq 0)\) y cuales no \((\beta_j = 0)\). Por su parte, ridge restringe todos los coeficientes hacia \(0\).

La Figura 8.1 muestra el valor de los coeficientes estandarizados en función de distintos valores del parámetro \(\lambda\) que se obtienen al aplicar lasso a la base de datos de Credit.¹⁹ Cuando \(\lambda = 0\), el resultado coincide con el de mínimos cuadrados, cuando \(\lambda\) es lo suficientemente grande, todos los coeficientes son iguales cero. Entre los dos casos extremos, moviéndose de derecha a izquierda se observa que al principio lasso resulta en un modelo que contiene solo el predictor de rating. Luego, ingresan al modelo student y limit casi simultáneamente, seguidas de cerca por income. Finalmente lo hacen el resto de las variables.

Figura 8.1: Lasso

Como señalamos arriba, si \(\hat{\beta_l} = 0\), al moverse de \(0\) (es decir, incluir la variable) podría pensarse que la función de penalidad de lasso genera costos y beneficios. Por un lado, el costo por incluir la variable sube en el valor de (\(\lambda\)), por el otro, la penalidad de ajuste baja (más variables, mejor ajuste). Si el primero sube más rápido que lo que el segundo baja conviene quedarse en \(0\). En caso contrario conviene moverse afuera de \(0\). Por lo tanto, conviene moverse de \(0\) si la relación es lo suficientemente fuerte (mejora el ajuste), sino conviene evitar la penalización.

8.2 Ridge

Nuevamente, la regresión ridge es muy similar a la de mínimos cuadrados. En particular, la estimación ridge \(\hat{\beta}^{R}\) son los valores que minimizan la ecuación (8.2) donde \(\lambda \ge 0\) es un parámetro de ajuste, que se determinará por separado. La ecuación (8.2) realiza un trade-off entre dos criterios diferentes. Al igual que con los mínimos cuadrados, la regresión ridge busca estimaciones de coeficientes que se ajusten bien a los datos, haciendo que RSS sea pequeña. Sin embargo, el segundo término \(\lambda \sum_{j=1}^{p}\beta_{j}^{2}\), llamado penalización, es pequeño cuando \(\beta_1,...,\beta_p\) están cerca de cero, por lo que tiene el efecto de reducir las estimaciones de \(\beta_j\) hacia cero.

\[\begin{equation} \tag{8.2} R_r(\beta) = \sum_{i=1}^{n}(y_i - x_i^{'} \beta)^2 + \lambda \sum_{j=1}^{p}\beta_{j}^{2} \end{equation}\]

Figura 8.2: Ridge

8.3 Aplicacion practica

Se utiliza el paquete glmnet para realizar la regresión de lasso²⁰ con una sintaxis ligeramente diferente de otras funciones de estimación de modelos. En particular, se debe usar una matriz x y un vector y y no utiliza la sintaxis (y ~ x). Se busca predecir Salary en la base de datos Hitters. Antes de comenzar deben eliminarse lo valores missing.²¹

library(ISLR2)
library(glmnet)
names(Hitters)

##  [1] "AtBat"     "Hits"      "HmRun"     "Runs"      "RBI"       "Walks"    
##  [7] "Years"     "CAtBat"    "CHits"     "CHmRun"    "CRuns"     "CRBI"     
## [13] "CWalks"    "League"    "Division"  "PutOuts"   "Assists"   "Errors"   
## [19] "Salary"    "NewLeague"

dim(Hitters)

## [1] 322  20

sum(is.na(Hitters$Salary))

## [1] 59

La Figura de abajo muestra la base de datos completa para analizar los valores missing. Las variables están en el eje y y las observaciones en el eje x.

hit = Hitters
rownames(hit) = 1:nrow(hit)
hit %>%
  is.na() %>%
  reshape2::melt() %>%
  ggplot(aes(Var2, Var1, fill=value)) + 
  geom_raster() + coord_flip() +
  scale_y_continuous(NULL, expand = c(0, 0)) +
  scale_fill_discrete(name = "", 
                  labels = c("Dato", 
                             "Missing")) +
  labs( x= 'Variable', title= 'Valores missing') +
  theme_minimal() +
  theme(axis.text.y  = element_text(size = 7)) +
  NULL

rm(hit)

Hitters <- na.omit(Hitters)
dim(Hitters)

## [1] 263  20

sum(is.na(Hitters))

## [1] 0

La función model.matrix() es útil para crear la matriz de datos x; no sólo produce una matriz correspondiente a los \(19\) predictores, sino también transforma automáticamente cualquier variable cualitativa en variables dummies (one-hot encoding). La última propiedad es importante porque glmnet() sólo puede tomar variables numéricas y cuantitativas.

x = model.matrix(Salary ~ ., Hitters)[, -1]
y = Hitters$Salary

Por primera vez dividimos la muestra en una base de entrenamiento y otra de test para estimar el error de test de la regresión lasso. Hay varias formas de dividir aleatoriamente una base de datos. Nuestro enfoque elige aleatoriamente un subconjunto de números entre \(1\) y \(n\); estos se pueden usar como índices para las observaciones de entrenamiento.

La función createDataPartition() de la librería caret es útil para crear particiones de train y test. La librería rsample permite realizar divisiones estratificadas (relevante cuando la clase está muy desbalanceada).

Primero establecemos una semilla aleatoria para que los resultados obtenidos sean reproducibles.

set.seed(1)
train <- sample(1:nrow(x), nrow(x) / 2) # selecciona la mitad de las observaciones
test <- (-train)
y.test <- y[test]

Para estimar un modelo lasso usamos la función glmnet() con el argumento alpha = 1 que indica la metodología que queremos usar.²²

También definimos una grilla²³ de \(100\) valores \(\lambda\) que van desde \(\lambda=10^{10}\) a \(\lambda=10^{-2}\), cubriendo esencialmente la gama completa de escenarios desde el modelo nulo que contiene solo la constante, hasta el modelo de mínimos cuadrados.

Podemos ver en el gráfico de abajo que cada curva corresponde a una variable. Muestra la trayectoria de su coeficiente frente a L1 norm (la penalidad de LASSO) del vector de coeficientes completo a medida que varía \(\lambda\). El eje de arriba indica el número de coeficientes distintos de cero para LASSO.

grid <- 10^seq(10, -2, length = 100)
lasso.mod <- glmnet(x[train, ], y[train], alpha = 1, lambda = grid)
plot(lasso.mod, label = TRUE)

Se muestra un resumen de la ruta de glmnet en cada paso:

print(lasso.mod)

## 
## Call:  glmnet(x = x[train, ], y = y[train], alpha = 1, lambda = grid) 
## 
##     Df  %Dev    Lambda
## 1    0  0.00 1.000e+10
## 2    0  0.00 7.565e+09
## 3    0  0.00 5.722e+09
## 4    0  0.00 4.329e+09
## 5    0  0.00 3.275e+09
## 6    0  0.00 2.477e+09
....

Ahora realizamos un ejercicio de cross validation (subsección 6.7.1) y calculamos el error de test asociado. Por defecto el argumento nfolds = 10. Los números en la parte superior de la Figura indican el número de variables en el modelo. La primera y segunda línea discontinua vertical representan el valor de \(\lambda\) con el \(EMC\) mínimo y el mayor valor de \(\lambda\) dentro de un error estándar.

set.seed(1)
cv.out <- cv.glmnet(x[train, ], y[train], alpha = 1)
plot(cv.out)

names(cv.out)

##  [1] "lambda"     "cvm"        "cvsd"       "cvup"       "cvlo"      
##  [6] "nzero"      "call"       "name"       "glmnet.fit" "lambda.min"
## [11] "lambda.1se" "index"

bestlam <- cv.out$lambda.min
lasso.pred <- predict(lasso.mod, s = bestlam, newx = x[test, ])
mean((lasso.pred - y.test)^2) # EMC

## [1] 143673.6

Finalmente, el resultado indica que \(10\) de los \(19\) coeficientes estimados son distintos de \(0\) (notar que hay una contradicción con las últimas dos tablas). Por lo tanto, el modelo de lasso con \(\lambda\) elegido por cross validation contiene solo \(11\) variables (más la constante).²⁴

cv.out$nzero[cv.out$index[[1]]]

## s36 
##  10

out <- glmnet(x, y, alpha = 1, lambda = grid)
lasso.coef <- predict(out, type = "coefficients", s = bestlam)[1:20, ]
lasso.coef

##   (Intercept)         AtBat          Hits         HmRun          Runs 
##    1.27479059   -0.05497143    2.18034583    0.00000000    0.00000000 
##           RBI         Walks         Years        CAtBat         CHits 
##    0.00000000    2.29192406   -0.33806109    0.00000000    0.00000000 
##        CHmRun         CRuns          CRBI        CWalks       LeagueN 
##    0.02825013    0.21628385    0.41712537    0.00000000   20.28615023 
##     DivisionW       PutOuts       Assists        Errors    NewLeagueN 
## -116.16755870    0.23752385    0.00000000   -0.85629148    0.00000000

lasso.coef[lasso.coef != 0]

##   (Intercept)         AtBat          Hits         Walks         Years 
##    1.27479059   -0.05497143    2.18034583    2.29192406   -0.33806109 
##        CHmRun         CRuns          CRBI       LeagueN     DivisionW 
##    0.02825013    0.21628385    0.41712537   20.28615023 -116.16755870 
##       PutOuts        Errors 
##    0.23752385   -0.85629148

Bibliografia

Simon, Noah, Jerome Friedman, Trevor Hastie, y Rob Tibshirani. 2011. «Regularization Paths for Cox’s Proportional Hazards Model via Coordinate Descent». Journal of Statistical Software 39 (5): 1-13. https://www.jstatsoft.org/v39/i05/.

Datos simulados que contiene información sobre 10.000 clientes donde el objetivo es predecir quienes no pagarán su deuda de tarjeta de crédito.↩︎
La función principal de este paquete es glmnet().↩︎
Para más detalles véase ?glmnet::glmnet y (Simon et al. 2011)↩︎
alpha = 0 estima con la metodología ridge.↩︎
También se pueden usar valores específicos.↩︎
Para más detalles de cv.glmnet ver aquí.↩︎