6.7 Métodos de regularización
Como ya se comentó, el procedimiento habitual para ajustar un modelo de regresión lineal es emplear mínimos cuadrados, es decir, utilizar como criterio de error la suma de cuadrados residual \[\mbox{RSS} = \sum\limits_{i=1}^{n}\left( y_{i} - \beta_0 - \boldsymbol{\beta}^t \mathbf{x}_{i} \right)^{2}\]
Si el modelo lineal es razonablemente adecuado, utilizar \(\mbox{RSS}\) va a dar lugar a estimaciones con poco sesgo, y si además \(n\gg p\), entonces el modelo también va a tener poca varianza (bajo las hipótesis estructurales, la estimación es insesgada y además de varianza mínima entre todas las técnicas insesgadas). Las dificultades surgen cuando \(p\) es grande o cuando hay correlaciones altas entre las variables predictoras: tener muchas variables dificulta la interpretación del modelo, y si además hay problemas de colinealidad o se incumple \(n\gg p\), entonces la estimación del modelo va a tener muchas varianza y el modelo estará sobreajustado. La solución pasa por forzar a que el modelo tenga menos complejidad para así reducir su varianza. Una forma de conseguirlo es mediante la regularización (regularization o shrinkage) de la estimación de los parámetros \(\beta_1, \beta_2,\ldots, \beta_p\) que consiste en considerar todas las variables predictoras pero forzando a que algunos de los parámetros se estimen mediante valores muy próximos a cero, o directamente con ceros. Esta técnica va a provocar un pequeño aumento en el sesgo pero a cambio una notable reducción en la varianza y una interpretación más sencilla del modelo resultante.
Hay dos formas básicas de lograr esta simplificación de los parámetros (con la consiguiente simplificación del modelo), utilizando una penalización cuadrática (norma \(L_2\)) o en valor absoluto (norma \(L_1\)):
Ridge regression (Hothorn et al., 2006) \[\mbox{min}_{\beta_0, \boldsymbol{\beta}} \mbox{RSS} + \lambda\sum_{j=1}^{p}\beta_{j}^{2}\]
Equivalentemente, \[\mbox{min}_{\beta_0, \boldsymbol{\beta}} \mbox{RSS}\] sujeto a \[\sum_{j=1}^{p}\beta_{j}^{2} \le s\]
LASSO (least absolute shrinkage and selection operator, Tibshirani, 1996) \[\mbox{min}_{\beta_0, \boldsymbol{\beta}} RSS + \lambda\sum_{j=1}^{p}|\beta_{j}|\]
Equivalentemente, \[\mbox{min}_{\beta_0, \boldsymbol{\beta}} \mbox{RSS}\] sujeto a \[\sum_{j=1}^{p}|\beta_{j}| \le s\]
Una formulación unificada consiste en considerar el problema \[\mbox{min}_{\beta_0, \boldsymbol{\beta}} RSS + \lambda\sum_{j=1}^{p}|\beta_{j}|^d\]
Si \(d=0\), la penalización consiste en el número de variables utilizadas, por tanto se corresponde con el problema de selección de variables; \(d=1\) se corresponde con LASSO y \(d=2\) con ridge.
La ventaja de utilizar LASSO es que va a forzar a que algunos parámetros sean cero, con lo cual también se realiza una selección de las variables más influyentes. Por el contrario, ridge regression va a incluir todas las variables predictoras en el modelo final, si bien es cierto que algunas con parámetros muy próximos a cero: de este modo va a reducir el riesgo del sobreajuste, pero no resuelve el problema de la interpretabilidad. Otra posible ventaja de utilizar LASSO es que cuando hay variables predictoras correlacionadas tiene tendencia a seleccionar una y anular las demás (esto también se puede ver como un inconveniente, ya que pequeños cambios en los datos pueden dar lugar a distintos modelos), mientras que ridge tiende a darles igual peso.
Dos generalizaciones de LASSO son least angle regression (LARS, Efron et al., 2004) y elastic net (Zou y Hastie, 2005). Elastic net combina las ventajas de ridge y LASSO, minimizando \[\mbox{min}_{\beta_0, \boldsymbol{\beta}} \ \mbox{RSS} + \lambda \left( \frac{1 - \alpha}{2}\sum_{j=1}^{p}\beta_{j}^{2} + \alpha \sum_{j=1}^{p}|\beta_{j}| \right)\] siendo \(0 \leq \alpha \leq 1\), un hiperparámetro adicional que determina la combinación lineal de ambas penalizaciones.
Es muy importante estandarizar (centrar y reescalar) las variables predictoras antes de realizar estas técnicas. Fijémonos en que, así como \(\mbox{RSS}\) es insensible a los cambios de escala, la penalización es muy sensible. Previa estandarización, el término independiente \(\beta_0\) (que no interviene en la penalización) tiene una interpretación muy directa, ya que \[\widehat \beta_0 = \bar y =\sum_{i=1}^n \frac{y_i}{n}\]
Los dos métodos de regularización comentados dependen del hiperparámetro \(\lambda\) (equivalentemente, \(s\)). Es muy importante seleccionar adecuadamente el valor del hiperparámetro, por ejemplo utilizando validación cruzada. Hay algoritmos muy eficientes que permiten el ajuste, tanto de ridge regression como de LASSO de forma conjunta (simultánea) para todos los valores de \(\lambda\).
6.7.1 Implementación en R
Hay varios paquetes que implementan estos métodos: h2o, elasticnet, penalized, lasso2, biglasso, etc., pero el paquete glmnet utiliza una de las más eficientes.
library(glmnet)El paquete glmnet no emplea formulación de modelos, hay que establecer la respuesta y y la matriz numérica x correspondiente a las variables explicativas.
Por tanto no se pueden incluir directamente predictores categóricos, habrá que codificarlos empleando variables auxiliares numéricas.
Se puede emplear la función model.matrix() (o Matrix::sparse.model.matrix() si el conjunto de datos es muy grande) para construir la matriz de diseño x a partir de una fórmula (alternativamente se pueden emplear la herramientas implementadas en el paquete caret).
Además, tampoco admite datos faltantes.
La función principal es glmnet():
glmnet(x, y, family, alpha = 1, lambda = NULL, ...)family: familia del modelo lineal generalizado (ver Sección 6.9); por defecto"gaussian"(modelo lineal con ajuste cuadrático), también admite"binomial","poisson","multinomial","cox"o"mgaussian"(modelo lineal con respuesta multivariante).alpha: parámetro \(\alpha\) de elasticnet \(0 \leq \alpha \leq 1\). Por defectoalpha = 1penalización LASSO (alpha = 0para ridge regression).lambda: secuencia (opcional) de valores de \(\lambda\); si no se especifica se establece una secuencia por defecto (en base a los argumentos adicionalesnlambdaylambda.min.ratio). Se devolverán los ajustes para todos los valores de esta secuencia (también se podrán obtener posteriormente para otros valores).
Entre los métodos genéricos disponibles del objeto resultante, coef() y predict() permiten obtener los coeficientes y las predicciones para un valor concreto de \(\lambda\), que se debe especificar mediante el argumento s = valor (“For historical reasons we use the symbol ‘s’ rather than ‘lambda’”).
Aunque para seleccionar el un valor “óptimo” del hiperparámetro \(\lambda\) (mediante validación cruzada) se puede emplear cv.glmnet():
cv.glmnet(x, y, family, alpha, lambda, type.measure = "default", nfolds = 10, ...)Esta función también devuelve los ajustes con toda la muestra de entrenamiento (en la componente $glmnet.fit) y se puede emplear el resultado directamente para predecir o obtener los coeficientes del modelo.
Por defecto seleccionando \(\lambda\) mediante la regla de “un error estándar” de Breiman et al. (1984) (componente $lambda.1se), aunque también calcula el valor óptimo (componente $lambda.min; que se puede seleccionar con estableciendo s = "lambda.min").
Para más detalles consultar la vignette del paquete “An Introduction to glmnet”.
Continuaremos con el ejemplo de los datos de clientes de la compañía de distribución industrial HBAT empleado en la Sección 6.1 (donde solo se consideraban predictores numéricos y tampoco había datos faltantes):
load("data/hbat.RData")
df <- hbat[c(23, 7:19)]
set.seed(1)
nobs <- nrow(df)
itrain <- sample(nobs, 0.8 * nobs)
train <- df[itrain, ]
test <- df[-itrain, ]
x <- as.matrix(train[-1])
y <- train$fidelida6.7.2 Ejemplo: Ridge Regression
Podemos ajustar modelos de regresión ridge (con la secuencia de valores de \(\lambda\) por defecto) con la función glmnet() con alpha=0 (ridge penalty).
Con el método plot() podemos representar la evolución de los coeficientes en función de la penalización (etiquetando las curvas con el índice de la variable si label = TRUE; ver Figura 6.11).
fit.ridge <- glmnet(x, y, alpha = 0)
plot(fit.ridge, xvar = "lambda", label = TRUE)
Figura 6.11: Gráfico de perfil de la evolución de los coeficientes en función del logaritmo de la penalización del ajuste ridge.
Podemos obtener el modelo o predicciones para un valor concreto de \(\lambda\):
coef(fit.ridge, s = 2) # lambda = 2## 14 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 3.56806743
## calidadp 2.41027431
## web 0.94414628
## soporte -0.22183509
## quejas 1.08417665
## publi 0.20121976
## producto 1.41018809
## imgfvent 0.21140360
## precio 0.26171759
## garantia 0.07110803
## nprod 0.04859325
## facturac 0.22695054
## flexprec 0.37732748
## velocida 3.11101217Para seleccionar el parámetro de penalización por validación cruzada empleamos cv.glmnet() (normalmente emplearíamos esta función en lugar de glmnet()).
El correspondiente método plot() muestra la evolución de los errores de validación cruzada en función de la penalización, incluyendo las bandas de un error estándar de Breiman (ver Figura 6.12).
set.seed(1)
cv.ridge <- cv.glmnet(x, y, alpha = 0)
plot(cv.ridge)
Figura 6.12: Error cuadrático medio de validación cruzada en función del logaritmo de la penalización del ajuste ridge, junto con los intervalos de un error estándar. Las líneas verticales se corresponden con lambda.min y lambda.1se.
En este caso el parámetro óptimo, según la regla de un error estándar de Breiman, sería37:
cv.ridge$lambda.1se## [1] 3.413705y el correspondiente modelo contiene todas las variables explicativas:
coef(cv.ridge) # s = "lambda.1se"## 14 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 8.38314266
## calidadp 2.06713538
## web 0.84771656
## soporte -0.17674893
## quejas 1.08099022
## publi 0.25926570
## producto 1.34198207
## imgfvent 0.21510001
## precio 0.15194226
## garantia 0.05417865
## nprod 0.08252518
## facturac 0.45964418
## flexprec 0.24646749
## velocida 2.70697234Finalmente evaluamos la precisión en la muestra de test:
obs <- test$fidelida
newx <- as.matrix(test[, -14])
pred <- predict(cv.ridge, newx = newx) # s = "lambda.1se"
accuracy(pred, obs)## me rmse mae mpe mape r.squared
## -108.0747 108.9324 108.0747 -186.9472 186.9472 -115.30466.7.3 Ejemplo: Lasso
También podríamos ajustar modelos LASSO con la opción por defecto de glmnet() (alpha = 1, lasso penalty).
Pero en este caso lo haremos al mismo tiempo que seleccionamos el parámetro de penalización por validación cruzada (ver Figura 6.13):
set.seed(1)
cv.lasso <- cv.glmnet(x,y)
plot(cv.lasso)
Figura 6.13: Error cuadrático medio de validación cruzada en función del logaritmo de la penalización del ajuste LASSO, junto con los intervalos de un error estándar. Las líneas verticales se corresponden con lambda.min y lambda.1se.
También podemos generar el gráfico con la evolución de los componentes a partir del ajuste almacenado en la componente $glmnet.fit:
plot(cv.lasso$glmnet.fit, xvar = "lambda", label = TRUE)
abline(v = log(cv.lasso$lambda.1se), lty = 2)
abline(v = log(cv.lasso$lambda.min), lty = 3)
Figura 6.14: Evolución de los coeficientes en función del logaritmo de la penalización del ajuste LASSO. Las líneas verticales se corresponden con lambda.min y lambda.1se.
Como podemos observar en la Figura 6.14, la penalización LASSO tiende a forzar que las estimaciones de los coeficientes sean exactamente cero cuando el parámetro de penalización \(\lambda\) es suficientemente grande. En este caso, el modelo resultante (empleando la regla oneSE) solo contiene 4 variables explicativas:
coef(cv.lasso) # s = "lambda.1se"## 14 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) 12.0485398
## calidadp 2.4673862
## web 0.3498592
## soporte .
## quejas .
## publi .
## producto 0.3227830
## imgfvent .
## precio .
## garantia .
## nprod .
## facturac .
## flexprec .
## velocida 6.1011015Por tanto este método también podría ser empleando para la selección de variables (puede hacerse automáticamente; estableciendo relax = TRUE, en la llamada a glmnet() o cv.glmnet(), devolverá los modelos ajustados sin regularización en la componente $relaxed).
Finalmente evaluamos también la precisión en la muestra de test:
pred <- predict(cv.lasso, newx = newx)
accuracy(pred, obs)## me rmse mae mpe mape r.squared
## -129.5991 130.8208 129.5991 -223.8971 223.8971 -166.74006.7.4 Ejemplo: Elastic Net
Podemos ajustar modelos elastic net para un valor concreto de alpha empleando la función glmnet(), pero las opciones del paquete no incluyen la selección de este hiperparámetro.
Aunque se podría implementar fácilmente (como se muestra en help(cv.glmnet)), resulta mucho más cómodo emplear el método "glmnet" de caret:
library(caret)
modelLookup("glmnet") ## model parameter label forReg forClass probModel
## 1 glmnet alpha Mixing Percentage TRUE TRUE TRUE
## 2 glmnet lambda Regularization Parameter TRUE TRUE TRUEset.seed(1)
# Se podría emplear train(fidelida ~ ., data = train, ...)
caret.glmnet <- train(x, y, method = "glmnet",
preProc = c("zv", "center", "scale"),
trControl = trainControl(method = "cv", number = 5),
tuneLength = 5)
caret.glmnet## glmnet
##
## 160 samples
## 13 predictor
##
## Pre-processing: centered (13), scaled (13)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 129, 129, 127, 127, 128
## Resampling results across tuning parameters:
##
## alpha lambda RMSE Rsquared MAE
## 0.100 0.005410605 4.581364 0.7148069 3.414825
## 0.100 0.025113806 4.576940 0.7153275 3.408862
## 0.100 0.116567961 4.545239 0.7187940 3.361951
## 0.100 0.541060544 4.474562 0.7284099 3.295198
## 0.100 2.511380580 4.704072 0.7187452 3.594686
## 0.325 0.005410605 4.573738 0.7157479 3.408931
## 0.325 0.025113806 4.564560 0.7167890 3.397543
## 0.325 0.116567961 4.500833 0.7241961 3.326005
## 0.325 0.541060544 4.438653 0.7349191 3.306102
## 0.325 2.511380580 4.881621 0.7184709 3.757854
## 0.550 0.005410605 4.573800 0.7157344 3.411370
## 0.550 0.025113806 4.552473 0.7182118 3.386635
## [ reached getOption("max.print") -- omitted 13 rows ]
##
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were alpha = 1 and lambda = 0.116568.Los resultados de la selección de los hiperparámetros \(\alpha\) y \(\lambda\) de regularización se muestran en la Figura 6.15:
ggplot(caret.glmnet, highlight = TRUE)
Figura 6.15: Errores RMSE de validación cruzada de los modelos elastic net en función de los hiperparámetros de regularización.
Finalmente, se evalúan las predicciones en la muestra de test del modelo ajustado (que en esta ocasión mejoran los resultados del modelo LASSO ajustando en la sección anterior):
pred <- predict(caret.glmnet, newdata = test)
accuracy(pred, obs)## me rmse mae mpe mape r.squared
## 0.49843131 4.28230543 3.43805155 -0.02851825 6.15711131 0.82026277References
Para obtener el valor óptimo global podemos emplear
cv.ridge$lambda.min, y añadir el argumentos = "lambda.min"a los métodoscoef()ypredict()para obtener los correspondientes coeficientes y predicciones.↩︎