Regularización

El objetivo de la regularización es la reducción de la varianza del modelo agregando penalizaciones a los coeficientes estimados.

library("rsample")
library("dplyr")
library("glmnet")
library("dplyr")
library("ggplot2")

# Creamos las muestras train test y plantamos una semilla pra reproducibildiad
set.seed(123)
ames_split <- initial_split(AmesHousing::make_ames(),prop = 0.7,
                            strata = "Sale_Price")

ames_training <- training(ames_split)
ames_test <- testing(ames_split)

¿ Por qué regularizar?

La meotodología de mínimos cuadrados ordinarios funciona bien cuando se cumplen una serie de supuestos:

• Relación lineal entre las variables
• Normalidad multivariada
• No autocorrelación
• Homocedasticidad (Varianza constante del error)
• Hay más observaciones que variables (n > p)
• No multicolinealidad

Por tal razón debería haber alternativas cuando los datos son muy distintos a aquellos que se ajustan al modelo clásico de regresión.

Una buena opción es la regresión regularizada o también llamada regresión de penelización o métodos de contracción para controlar los parametros estimados. La idea de la regresión regularizada es poner restricciones de magnitud a los coeficientes y contraerlos progresivamente hacia cero. Esta restricción ayuda a reducir la varianza del modelo .

Regresión regularizada

La función objetivo de la regresión regularizada es similar a MCO pero con una penalización adicional:

\[ min\{RSS + p\} \]

Al aplicar la penalización a la suma cuadrada de los errores se restringe el tamaño de los coeficientes tal que la unica manera en que los coeficientes puedan aumentar es que la inclusión de una variables sea comparable con una disminución significativa del error.

Regresión Ridge

Para la estimación de los coeficientes en mínimos cuadrados debemos minimizar la suma de los errores al cuadrado. Para generar una regresión tipo rigde agregamos la penalización y mínimizados la expresión:

\[ \sum_{i=1}^n{(y_i - \beta_o - \sum_{j=1}^p{\beta_jx_{ij}})^2} + \lambda\sum_{j=1}^p{\beta_{j}^2} \] Donde la primera expresión es la suma de los errores al cuadrado y es un parametro que debe ser tuneado.

# Partimos las tablas en variables depentientes de independientes


#La función model.matrix expande los datos a una matrix transformando factores a dummies
ames_train_x <- model.matrix(Sale_Price ~ ., ames_training)[, -1]
ames_train_y <- log(ames_training$Sale_Price)

ames_test_x <- model.matrix(Sale_Price ~ ., ames_test)[, -1]
ames_test_y <- log(ames_test$Sale_Price)

1. Esencial que los predictores estén reescalados o estandarizados antes de realizar la regresión regularizada

• alpha=0 es para regresión ridge
• alpha=1 es para regresión lasso
• 0<alpha<1 es para elastic net

# Por defecto aplica el lambda 100 veces
ames_ridge <- glmnet(
  x = ames_train_x,
  y = ames_train_y,
  alpha = 0
  
)

plot(ames_ridge, xvar = "lambda")

# Lambdas
ames_ridge$lambda[1:5]

## [1] 285.8055 260.4153 237.2807 216.2014 196.9946

# coeficientes más pequeños
coef(ames_ridge)[c("Gr_Liv_Area", "TotRms_AbvGrd"), 100]

##   Gr_Liv_Area TotRms_AbvGrd 
##  0.0001100975  0.0095201380

#Coeficientes más grandes
coef(ames_ridge)[c("Gr_Liv_Area", "TotRms_AbvGrd"), 1]

##   Gr_Liv_Area TotRms_AbvGrd 
##  5.838870e-40  1.332454e-37

Sin embargo aún no sabemos con seguridad cuál lambda garantiza que tengamos el mejor modelo. Por ello realizados cross-validation con 10 pliegues

# REgresión ridge con 10-fold cross validation
ames_ridge_cv <- cv.glmnet(
  x = ames_train_x,
  y = ames_train_y,
  alpha = 0
)

plot(ames_ridge_cv)

min(ames_ridge_cv$cvm)       # minimo MSE

## [1] 0.01752896

ames_ridge_cv$lambda.min     # lambda que minimiza la MSE

## [1] 0.1051301

ames_ridge_min <- glmnet(
  x = ames_train_x,
  y = ames_train_y,
  alpha = 0
)

plot(ames_ridge_min, xvar = "lambda")
abline(v = log(ames_ridge_cv$lambda.1se), col = "red", lty = "dashed")

Variables que mayor influcencia tienen en el modelo (Por el tamaño de sus coeficiente)

coef(ames_ridge_cv, s = "lambda.1se") %>%
  broom::tidy() %>%
  filter(row != "(Intercept)") %>%
  top_n(25, wt = abs(value)) %>%
  ggplot(aes(value, reorder(row, value))) +
  geom_point() +
  ggtitle("Top 25 variables con mayor influencia") +
  xlab("Coeficiente") +
  ylab(NULL) + theme_classic()

Regresión Lasso ( least absolute shrinkage and selection operator)

A diferencia de de ridge, matemáticamente el único cambio es que ahora los coeficientes de la penalización están en valor absoluto en vez de elevados al cuadrado. Esto tiene efectos distintos a la ridge:

• No castiga de la misma manera a los coeficientes muy grandes
• En la regresión ridge los coeficientes decaian hacia cero, es Lasso los coeficientes puede volverse cero!, lo que implica que la regresión lasso tiene otro efecto y es que automaticamente depura las variables que on agregan poder preditivo al modelo.

\[ \sum_{i=1}^n{(y_i - \beta_o - \sum_{j=1}^p{\beta_jx_{ij}})^2} + \lambda\sum_{j=1}^p{|\beta_{j}|} \]

## Apply lasso regression to ames data
ames_lasso <- glmnet(
  x = ames_train_x,
  y = ames_train_y,
  alpha = 1
)

plot(ames_lasso, xvar = "lambda")

Al igual que antes aplicamos cross validation

# Apply CV Ridge regression to ames data
ames_lasso_cv <- cv.glmnet(
  x = ames_train_x,
  y = ames_train_y,
  alpha = 1
)
# plot results
plot(ames_lasso_cv)

min(ames_lasso_cv$cvm)       # MSE minimo

## [1] 0.01908478

ames_lasso_cv$lambda.min     # lambda que minimiza el MSE

## [1] 0.002994143

ames_lasso_min <- glmnet(
  x = ames_train_x,
  y = ames_train_y,
  alpha = 1
)

plot(ames_lasso_min, xvar = "lambda")
abline(v = log(ames_lasso_cv$lambda.min), col = "red", lty = "dashed")
abline(v = log(ames_lasso_cv$lambda.1se), col = "red", lty = "dashed")

Variables de mayor influencia en la regresión Lasso

coef(ames_lasso_cv, s = "lambda.1se") %>%
  tidy() %>%
  filter(row != "(Intercept)") %>%
    top_n(25, wt = abs(value)) %>%
  ggplot(aes(value, reorder(row, value), color = value > 0)) +
  geom_point(show.legend = FALSE) +
  ggtitle("Variables de influencia en el Modelo Lasso") +
  xlab("Coeficiente") +
  ylab(NULL) +theme_classic()

## Warning: 'tidy.dgCMatrix' is deprecated.
## See help("Deprecated")

## Warning: 'tidy.dgTMatrix' is deprecated.
## See help("Deprecated")

Sin embargo al remover variables perdemos precisión

# minimum Ridge MSE
min(ames_ridge_cv$cvm)

## [1] 0.01752896

## [1] 0.02147691

# minimum Lasso MSE
min(ames_lasso_cv$cvm)

## [1] 0.01908478

## [1] 0.02275227

Regresión Elastic Net

Es una combinación de Ridge y Lasso. Se decide entonces qué peso se le da a cada método de penalización y se implementa la regresión

\[ \sum_{i=1}^n{(y_i - \beta_o - \sum_{j=1}^p{\beta_jx_{ij}})^2} + \lambda_{1}\sum_{j=1}^p{\beta_{j}^2} +\lambda_{2}\sum_{j=1}^p{|\beta_{j}|} \]

lasso    <- glmnet(ames_train_x, ames_train_y, alpha = 1.0) 
elastic1 <- glmnet(ames_train_x, ames_train_y, alpha = 0.25) 
elastic2 <- glmnet(ames_train_x, ames_train_y, alpha = 0.75) 
ridge    <- glmnet(ames_train_x, ames_train_y, alpha = 0.0)

par(mfrow = c(2, 2), mar = c(6, 4, 6, 2) + 0.1)
plot(lasso, xvar = "lambda", main = "Lasso (Alpha = 1)\n\n\n")
plot(elastic1, xvar = "lambda", main = "Elastic Net (Alpha = .25)\n\n\n")
plot(elastic2, xvar = "lambda", main = "Elastic Net (Alpha = .75)\n\n\n")
plot(ridge, xvar = "lambda", main = "Ridge (Alpha = 0)\n\n\n")

Tuning

En modelos Ridge y Lasso debemos tunear el parametro \(\lambda\). En elastic Net hay que tunear también el parametro alpha, que nos indica el paso correcto que se le debe dedat a cada penalización.

# maintain the same folds across all models
fold_id <- sample(1:10, size = length(ames_train_y), replace=TRUE)

# search across a range of alphas
tuning_grid <- tibble::tibble(
  alpha      = seq(0, 1, by = .1),
  mse_min    = NA,
  mse_1se    = NA,
  lambda_min = NA,
  lambda_1se = NA
)

for(i in seq_along(tuning_grid$alpha)) {
  
  # fit CV model for each alpha value
  fit <- cv.glmnet(ames_train_x, ames_train_y, alpha = tuning_grid$alpha[i], foldid = fold_id)
  
  # extract MSE and lambda values
  tuning_grid$mse_min[i]    <- fit$cvm[fit$lambda == fit$lambda.min]
  tuning_grid$mse_1se[i]    <- fit$cvm[fit$lambda == fit$lambda.1se]
  tuning_grid$lambda_min[i] <- fit$lambda.min
  tuning_grid$lambda_1se[i] <- fit$lambda.1se
}

tuning_grid %>%
  mutate(se = mse_1se - mse_min) %>%
  ggplot(aes(alpha, mse_min)) +
  geom_line(size = 2) +
  geom_ribbon(aes(ymax = mse_min + se, ymin = mse_min - se), alpha = .25) +
  ggtitle("MSE ± un error estandar") + theme_minimal()

Predicciones

Para finalizar, las predicciones de los modelos solo requiren que se especifique el modelo entrenado

# Modelo Lasso con cross validation
cv_lasso   <- cv.glmnet(ames_train_x, ames_train_y, alpha = 1.0)
min(cv_lasso$cvm)

## [1] 0.01912533

# Predicción del modelo con el lambda seleccionado
pred <- predict(cv_lasso, s = cv_lasso$lambda.min, ames_test_x)

# Error del modelo
mean((ames_test_y - pred)^2)

## [1] 0.02684147

Nota final

Este documento es una copia del artículo publicado por la unversidad de Cincinnati: http://uc-r.github.io/regularized_regression#req