Como vimos en el ejemplo anterior, en algunos casos podemos construir modelos lineales de complejidad considerable (por ejemplo, transformando variables, incluyendo interacciones). De manera que aún cuando muchas veces se considera un modelo lineal como “simple” o de “baja complejidad”, es posible que la variabilidad de las estimaciones sea grande y sobreajustar.
En este caso, minimizar el error cuadrático medio sobre la muestra de entrenamiento como vimos en la sección anterior puede no ser muy buena idea. En esta parte veremos una estrategia básica para “limitar” el grado de ajuste que nuestro modelo lineal puede alcanzar con datos de entrenamiento, esperando que esto mejore el desempeño predictivo. En otras palabras, modificamos nuestra función objetivo para que las \(f(x)\) no sobreajusten a las \(y\). Este primera estrategia es llamada regularización, y consiste agregar una penalización por tamaño de los coeficientes a la función objetivo.
Veamos primero un ejemplo simulado.
Consideremos un problema donde tenemos unas 100 entradas con 120 casos. Supondremos que la función verdadera es
\[f(x) = \sum_{j=1}^{100} \beta_j x_j\] con ciertos pesos o coeficientes \(\beta\) que fijaremos.
library(tidyverse)
library(gt)
set.seed(28015)
beta_vec <- rnorm(100, 0, 0.2)
p <- length(beta_vec)
beta <- tibble(term = str_c('V', 1:p), valor = beta_vec)
head(beta)# A tibble: 6 × 2
term valor
<chr> <dbl>
1 V1 -0.121
2 V2 0.0374
3 V3 -0.129
4 V4 0.240
5 V5 -0.00962
6 V6 -0.0443 Simulamos datos:
sim_datos <- function(n, beta){
p <- nrow(beta)
mat_x <- matrix(rnorm(n * p, 0, 1), n, p) + rnorm(n, 0, 5)
colnames(mat_x) <- beta |> pull(term)
beta_vec <- beta |> pull(valor)
f_x <- mat_x %*% beta_vec
y <- as.numeric(f_x) + rnorm(n, 0, 1)
datos <- as_tibble(mat_x)
datos |> mutate(y = y)
}
datos <- sim_datos(n = 4000, beta = beta)Separamos datos de entrenamiento y prueba y definimos y ajustamos un predictor lineal:
library(tidymodels)
set.seed(994)
n_entrena <- nrow(datos) * 0.03
separacion <- initial_split(datos, 0.03)
dat_ent <- training(separacion)
modelo <- linear_reg() |> set_engine("lm")
receta <- recipe(y ~ ., dat_ent)
flujo <- workflow() |>
add_model(modelo) |>
add_recipe(receta)
flujo_ajustado <- fit(flujo, dat_ent)
mod_1 <- flujo_ajustado |> extract_fit_engine()Extraemos los coeficientes y graficamos ajustados contra verdaderos:
coefs_1 <- tidy(mod_1) |>
left_join(beta, by = "term")
ggplot(coefs_1 |> filter(term != "(Intercept)"),
aes(x = valor, y = estimate)) +
geom_point() +
xlab('Coeficientes verdaderos') +
ylab('Coeficientes estimados') +
geom_abline() 
Y notamos que las estimaciones no son buenas. Podemos hacer otra simulación para confirmar que el problema es que las estimaciones son muy variables.
Con otra muestra de entrenamiento, vemos que las estimaciones tienen varianza alta.
datos_ent_2 <- sim_datos(n = 120, beta = beta)
mod_2 <- fit(flujo, datos_ent_2) |> extract_fit_engine()
coefs_2 <- tidy(mod_2)
qplot(coefs_1$estimate, coefs_2$estimate) + xlab('Coeficientes mod 1') +
ylab('Coeficientes mod 2') +
geom_abline(intercept=0, slope =1)Warning: `qplot()` was deprecated in ggplot2 3.4.0.
En la práctica, nosotros tenemos una sola muestra de entrenamiento. Así que, con una muestra de tamaño
\(n=120\) como en este ejemplo, obtendremos típicamente resultados no muy buenos. Estos coeficientes ruidosos afectan nuestras predicciones de manera negativa, aún cuando el modelo ajustado parece reproducir razonablemente bien la variable respuesta:
library(patchwork)
dat_pr <- testing(separacion)
preds_entrena <- predict(flujo_ajustado, dat_ent) |>
bind_cols(dat_ent |> select(y))
preds_prueba <- predict(flujo_ajustado, dat_pr) |>
bind_cols(dat_pr |> select(y))
g_1 <- ggplot(preds_entrena, aes(x = .pred, y = y)) +
geom_abline(colour = "red") +
geom_point() +
xlab("Predicción") + ylab("y") +
labs(subtitle = "Muestra de entrenamiento")
g_2 <- ggplot(preds_prueba, aes(x = .pred, y = y)) +
geom_abline(colour = "red") +
geom_point() +
xlab("Predicción") + ylab("y") +
labs(subtitle = "Muestra de prueba")
g_1 + g_2
Como el problema es la variabilidad de los coeficientes (en este ejemplo sabemos que no hay sesgo pues conocemos el modelo verdadero), podemos atacar este problema poniendo restricciones a los coeficientes, de manera que caigan en rangos más aceptables.
Una manera de hacer esto es restringir el rango de los coeficientes cambiando la función que minimizamos para ajustar el modelo lineal. Recordamos que la cantidad que queremos minimizar es
\[D(\beta) = D(a_0, \beta_1, \ldots, \beta_p) = \sum_{i=1}^N (y^{(i)} - f_\beta (x^{(i)}))^2 = \sum_{i=1}^N (y^{(i)} - \beta_0 - \beta_1 x_1^{(i)}-\beta_2x_2^{(i)} - \cdots - \beta_px_p^{(i)})^2\]
donde la suma es sobre los datos de entrenamiento. Queremos encontrar \(a =(\beta_0, \beta_1, \ldots, \beta_p)\) para resolver
\[\min_\beta D(\beta)\]
En el ejemplo que estamos considerando, vemos que existe mucha variación en los coeficientes obtenidos de muestra de entrenamiento a muestra de entrenamiento, y que algunos de ellos toman valores muy grandes positivos o negativos. Podemos entonces intentar resolver mejor el problema penalizado
\[\min_\beta \left\{ D(\beta) + \lambda \sum_{j=1}^p \beta_j^2 \right\}\]
Si escogemos un valor relativamente grande de \(\lambda > 0\), entonces terminaremos con una solución donde los coeficientes
no pueden alejarse mucho de 0, y esto previene parte del sobreajuste que observamos en nuestro primer ajuste. Otra manera de decir esto es: intentamos minimizar cuadrados, pero no permitimos que los coeficientes se alejen demasiado de cero, o ponemos un costo a soluciones que intentan “mover” mucho los coeficientes para ajustar mejor al conjunto de entrenamiento.
En este caso obtenemos:
modelo_reg <- linear_reg(mixture = 0, penalty = 0.1) |>
set_engine("glmnet", lambda.min.ratio = 0)
flujo_reg <- workflow() |>
add_model(modelo_reg) |>
add_recipe(receta)
flujo_reg <- fit(flujo_reg, dat_ent)
mod_reg <- flujo_reg |> extract_fit_parsnip()Los coeficientes del modelo penalizado son:
coefs_penalizado <- tidy(mod_reg)
Attaching package: 'Matrix'The following objects are masked from 'package:tidyr':
expand, pack, unpackLoaded glmnet 4.1-8coefs_penalizado# A tibble: 101 × 3
term estimate penalty
<chr> <dbl> <dbl>
1 (Intercept) -0.236 0.1
2 V1 -0.0774 0.1
3 V2 0.0340 0.1
4 V3 -0.0405 0.1
5 V4 0.151 0.1
6 V5 0.150 0.1
7 V6 0.114 0.1
8 V7 -0.207 0.1
9 V8 0.0148 0.1
10 V9 0.191 0.1
# ℹ 91 more rowsNótese que efectivamente la suma de cuadrados de los coeficientes penalizados es considerablemente más chica que las del modelo no penalizado:
sum(coefs_penalizado$estimate[-1]^2)[1] 1.957447sum(coefs_1$estimate[-1]^2)[1] 12.59754Los nuevos coeficientes estimados tienen menor variación, y están más cercanos a los valores reales:
qplot(x = beta$valor, coefs_penalizado$estimate[-1]) +
xlab('Coeficientes') +
ylab('Coeficientes estimados') +
geom_abline()
preds_prueba_2 <- predict(mod_reg, dat_pr) |>
bind_cols(dat_pr |> select(y))
preds_prueba_ambas <- bind_rows(
preds_prueba |> mutate(tipo = "sin penalizar (prueba)"),
preds_prueba_2 |> mutate(tipo = "penalizado (prueba)"))
ggplot(preds_prueba_ambas, aes(x = .pred, y = y)) +
geom_abline(colour = "red") +
geom_point(alpha = 0.3) +
xlab("Predicción") + ylab("y") +
facet_wrap(~ tipo, nrow = 1) +
labs(subtitle = "Muestra de prueba") 
metricas <- metric_set(mae, rmse)
res_1 <- metricas(preds_prueba, truth = y, estimate = .pred) |>
mutate(tipo = "no penalizado")
res_2 <- metricas(preds_prueba_2, truth = y, estimate = .pred) |>
mutate(tipo = "penalizado")
bind_rows(res_1, res_2) |>
arrange(.metric) |>
gt() |> fmt_number(.estimate, decimals = 2)| .metric | .estimator | .estimate | tipo |
|---|---|---|---|
| mae | standard | 2.34 | no penalizado |
| mae | standard | 1.27 | penalizado |
| rmse | standard | 2.95 | no penalizado |
| rmse | standard | 1.60 | penalizado |
Y vemos que los errores de predicción se reducen considerablemente.
Cuando agregamos el término de penalización tipo ridge al error de entrenamiento como objetivo a minimizar en el ajuste, los coeficientes de la solución penalizada están encogidos con respecto a los no penalizados.
Regularizar reduce la varianza de los coeficientes a lo largo de distintas muestras de entrenamiiento, lo que reduce la posibilidad de sobreajuste.
Utilizamos regularización para reducir el error de predicción cuando el problema es variabilidad grande de los coeficientes (coeficientes ruidosos) en modelos relativamente grandes o con pocos datos de entrenamiento.
En general, a métodos donde restringimos el espacio de modelos o penalizamos ajustes complejos en la función de pérdida que nos interesa se llaman métodos con regularización. Un ejemplos es todos los modelos donde en lugar de considerar la función de perdida \(L\) solamente, consideramos minimizar
\[L(f) + \Omega(f),\] donde \(f\) es una medida de la complejidad, como puede ser: que la función \(f\) tiene oscilaciones grandes o pendientes grandes, tiene un número grande de discontinuidades, etc.
Consideremos los siguientes datos clásicos de Radiación Solar, Temperatura, Velocidad del Viento y Ozono para distintos días en Nueva York (Chambers et al. (1983)):
air_data <- airquality |>
mutate(Wind_cat = cut(Wind, quantile(Wind, c(0, 0.33, 0.66, 1)),
include.lowest = T)) |>
filter(!is.na(Ozone) & !is.na(Solar.R))
air <- air_data
ggplot(air, aes(x = Solar.R, y = Ozone, colour = Temp)) +
geom_point() +
facet_wrap(~Wind_cat, ncol = 3) +
scale_colour_gradientn(colours = rainbow(2, rev = TRUE))
La gráfica muestra algunas interacciones y relaciones no lineales. Formulamos un modelo lineal como sigue:
receta_ozono <- recipe(Ozone ~ Temp + Wind + Solar.R,
data = air) |>
step_ns(Temp, Wind, Solar.R, deg_free = 2) |>
step_interact(terms = ~ starts_with("Temp_ns"):starts_with("Wind_ns")) |>
step_interact(terms = ~ starts_with("Temp_ns"):starts_with("Solar.R_ns"))
ajuste_ozono <- workflow() |>
add_recipe(receta_ozono) |>
add_model(linear_reg() |> set_engine("lm")) |>
fit(air)Y el ajuste se ve como sigue:
pred_grid <- expand_grid(Wind = c(5,10,15),
Temp = seq(60, 90, 10),
Solar.R = seq(20, 300, by = 10)) |>
mutate(Wind_cat = cut(Wind,
quantile(airquality$Wind, c(0, 0.33, 0.66, 1)),
include.lowest = T))
pred_grid <- pred_grid |>
bind_cols(predict(ajuste_ozono, pred_grid))
g_lineal <- ggplot(air, aes(x = Solar.R, colour = Temp)) +
geom_point(aes(y = Ozone)) +
facet_wrap( ~ Wind_cat) +
scale_colour_gradientn(colours = rainbow(2, rev = TRUE)) +
geom_line(data = pred_grid,
aes(y = .pred, group = interaction(Temp, Wind_cat)), size = 1) +
labs(subtitle = "Curvas de modelo lineal, para viento = 5, 10, 15") Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.g_lineal
Nótese que algunos aspectos de este modelo parecen ser muy ruidosos: por ejemplo, el comportamiento de las curvas para el primer pánel (donde hay pocos datos de temperatura baja), el hecho de que en algunos casos parece haber curvaturas decrecientes e incluso predicciones negativas. No deberíamos dar mucho crédito a las predicciones de este modelo, y tiene peligro de producir predicciones desastrosas.
Sin embargo, si usamos algo de regularización:
ajuste_ozono <- workflow() |>
add_recipe(receta_ozono) |>
add_model(linear_reg(mixture = 0, penalty = 3.0) |>
set_engine("glmnet", lambda.min.ratio = 0)) |>
fit(air)
# nota: normalmente no es necesario usar lambda.min.ratioY el ajuste se ve como sigue:
pred_grid <- expand_grid(Wind = c(5,10,15),
Temp = seq(60, 90, 10),
Solar.R = seq(10, 320, by = 10)) |>
mutate(Wind_cat =
cut(Wind, quantile(airquality$Wind, c(0, 0.33, 0.66, 1)),
include.lowest = T))
pred_grid <- pred_grid |>
bind_cols(predict(ajuste_ozono, pred_grid))
g_lineal <- ggplot(air, aes(x = Solar.R, colour = Temp)) +
geom_point(aes(y = Ozone)) +
facet_wrap( ~ Wind_cat) +
scale_colour_gradientn(colours = rainbow(2, rev = TRUE)) +
geom_line(data = pred_grid, aes(y = .pred, group = interaction(Temp, Wind_cat)), size = 1) +
labs(subtitle = "Curvas de modelo lineal, para viento = 5, 10, 15")
g_lineal
Este ajuste se ve mucho más razonable.
Como vimos antes, no es posible seleccionar el parámetro \(\lambda\) usando la muestra de entrenamiento (¿con qué \(\lambda\) cómo se obtiene el menor error cuadrático medio sobre la muestra de entrenamiento). Usaremos un conjunto de validación relativamente grande
set.seed(191)
source("../R/casas_traducir_geo.R")Rows: 1460 Columns: 81
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (43): MSZoning, Street, Alley, LotShape, LandContour, Utilities, LotConf...
dbl (38): Id, MSSubClass, LotFrontage, LotArea, OverallQual, OverallCond, Ye...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Rows: 27 Columns: 3
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): Neighborhood
dbl (2): lat, long
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.# esta proporción es para ejemplificar
casas_split <- initial_split(casas, prop = 0.25)
casas_entrena <- training(casas_split)
receta_casas <-
recipe(precio_miles ~ calidad_gral +
area_hab_m2 + area_garage_m2 + area_sotano_m2 +
area_2o_piso_m2 +
año_construccion + año_venta + condicion_venta +
nombre_zona +
condicion_gral +
condicion_exteriores +
tipo_sotano + calidad_sotano +
baños_completos + num_coches +
aire_acondicionado +
tipo_edificio + estilo,
data = casas_entrena) |>
step_filter(condicion_venta == "Normal") |>
step_select(-condicion_venta, skip = TRUE) |>
step_cut(calidad_gral, breaks = c(3, 5, 7, 8),
include_outside_range = TRUE) |>
step_cut(condicion_gral, breaks = c(3, 5, 7, 8),
include_outside_range = TRUE) |>
step_mutate(sin_piso_2 = as.numeric(area_2o_piso_m2 == 0)) |>
step_novel(tipo_sotano, calidad_sotano) |>
step_novel(condicion_exteriores, tipo_edificio, estilo) |>
step_unknown(tipo_sotano, calidad_sotano, new_level = "sin_sotano") |>
step_unknown(condicion_exteriores) |>
step_other(nombre_zona, threshold = 0.05) |>
step_dummy(calidad_gral, condicion_gral,
nombre_zona, aire_acondicionado,
calidad_sotano, tipo_sotano, condicion_exteriores,
tipo_edificio, estilo) |>
step_interact(terms = ~ c(area_hab_m2, area_garage_m2,
area_sotano_m2, area_2o_piso_m2):starts_with("calidad_gral")) |>
step_interact(terms = ~ area_sotano_m2:starts_with("calidad_sotano")) |>
step_interact(terms = ~ c(area_hab_m2, area_garage_m2,
area_sotano_m2, area_2o_piso_m2):starts_with("condicion_gral")) |>
step_interact(terms = ~ c(area_hab_m2, area_garage_m2,
area_sotano_m2, area_2o_piso_m2):starts_with("nombre_zona")) |>
step_zv(all_predictors())Para ver el número de entradas de este modelo:
prep(receta_casas) |> juice() |> dim()[1] 289 108modelo_penalizado <- linear_reg(mixture = 0.0, penalty = tune()) |>
set_engine("glmnet", lambda.min.ratio = 0)
flujo_casas <- workflow() |>
add_recipe(receta_casas) |>
add_model(modelo_penalizado)Construimos manualmente el conjunto de validación:
# creamos un objeto con datos de entrenamiento y de prueba
val_split <- manual_rset(casas_split |> list(), "validación")
lambda_params <- parameters(penalty(range = c(-3, 3),
trans = log10_trans()))
lambda_grid <- grid_regular(lambda_params, levels = 20)
lambda_grid# A tibble: 20 × 1
penalty
<dbl>
1 0.001
2 0.00207
3 0.00428
4 0.00886
5 0.0183
6 0.0379
7 0.0785
8 0.162
9 0.336
10 0.695
11 1.44
12 2.98
13 6.16
14 12.7
15 26.4
16 54.6
17 113.
18 234.
19 483.
20 1000 mis_metricas <- metric_set(rmse)
eval_tbl <- tune_grid(flujo_casas,
resamples = val_split,
grid = lambda_grid,
metrics = mis_metricas)
ridge_ajustes_tbl <- eval_tbl |>
unnest(cols = c(.metrics)) |>
select(id, penalty, .metric, .estimate)
ridge_ajustes_tbl |> gt()| id | penalty | .metric | .estimate |
|---|---|---|---|
| validación | 1.000000e-03 | rmse | 35.36721 |
| validación | 2.069138e-03 | rmse | 35.36721 |
| validación | 4.281332e-03 | rmse | 35.36721 |
| validación | 8.858668e-03 | rmse | 35.36721 |
| validación | 1.832981e-02 | rmse | 35.36721 |
| validación | 3.792690e-02 | rmse | 35.36721 |
| validación | 7.847600e-02 | rmse | 34.88798 |
| validación | 1.623777e-01 | rmse | 33.64199 |
| validación | 3.359818e-01 | rmse | 32.41837 |
| validación | 6.951928e-01 | rmse | 31.22570 |
| validación | 1.438450e+00 | rmse | 30.26859 |
| validación | 2.976351e+00 | rmse | 29.62744 |
| validación | 6.158482e+00 | rmse | 29.45135 |
| validación | 1.274275e+01 | rmse | 29.65517 |
| validación | 2.636651e+01 | rmse | 30.23900 |
| validación | 5.455595e+01 | rmse | 31.52315 |
| validación | 1.128838e+02 | rmse | 34.17121 |
| validación | 2.335721e+02 | rmse | 38.82727 |
| validación | 4.832930e+02 | rmse | 45.67700 |
| validación | 1.000000e+03 | rmse | 54.01752 |
ggplot(ridge_ajustes_tbl,
aes(x = penalty, y = .estimate, colour = .metric)) +
geom_point() + geom_line() + scale_x_log10() 
Y vemos que con una penalización alrededor de \(\lambda = 1\) podemos obtener mejor desempeño que con el modelo no regularizado.
Pregunta: en qué partes de la gráfica es relativamente grande la varianza? ¿en qué parte es relativamente grande el sesgo?
Se puede controlar la varianza de mínimos cuadrados de otras maneras. Cuando la varianza proviene también de la inclusión de variables que no necesariamente están relacionadas con la respuesta, podemos usar métodos de selección de variables, como en stepwise regression, por ejemplo.
Otra manera interesante de lograr mejor desempeño predictivo con modelos más parsimoniosos resulta de usar un término de penalización distinto al de ridge. En ridge, el problema que resolvemos es minimizar el objetivo
\[D(\beta) + \lambda \sum_{j=1}^p \beta_j^2\]
En regresión lasso, usamos una penalización de tipo \(L_1\):
\[D(a) + \lambda \sum_{j=1}^p |\beta_j|\] En un principio, no parece ser muy diferente a ridge. Veremos sin embargo que usar esta penalización también se puede ver como un proceso de selección de variables.
Consideramos cómo predecir el porcentaje de grasa corporal a partir de distintas mediciones de dimensiones corporales:
dat_grasa <- read_csv(file = '../datos/bodyfat.csv')
set.seed(183)
grasa_particion <- initial_split(dat_grasa, 0.7)
grasa_ent <- training(grasa_particion)
grasa_pr <- testing(grasa_particion)# nota: con glmnet no es necesario normalizar, pero aquí lo hacemos
# para ver los coeficientes en términos de las variables estandarizadas:
grasa_receta <- recipe(grasacorp ~ ., grasa_ent) |>
update_role(cadera, cuello, muñeca,
tobillo, rodilla, new_role = "ninguno") |>
step_normalize(all_predictors())
modelo_2 <- linear_reg(mixture = 0, penalty = 0) |>
set_engine("glmnet", lambda.min.ratio = 1e-20)
flujo_2 <- workflow() |>
add_model(modelo_2) |>
add_recipe(grasa_receta)
flujo_2 <- flujo_2 |> fit(grasa_ent)
modelo_2 <- extract_fit_parsnip(flujo_2)
coefs <- modelo_2 |> pluck("fit") |>tidy() |>
filter(term != "(Intercept)")
g_l2 <- ggplot(coefs, aes(x = lambda, y = estimate, colour = term)) +
geom_line(size = 1.4) + scale_x_log10() +
scale_colour_manual(values = cbb_palette) +
labs(subtitle = "Regularizacion L2")
g_l2
Estas gráfica se llama traza de los coeficientes, y nos muestra cómo cambian los coefi´cientes conforme cambiamos la regularización. Nótese que cuando la regularización es chica, obtenemos algunos resultados contra-intuitivos como que el coeficiente de peso es negativo para predecir el nivel de grasa corporal. Cuando regularizamos más, este coeficiente es positivo. La razón de esto tiene qué ver con la correlación fuerte de las variables de entrada, por ejemplo:
cor(grasa_ent |> select(peso, abdomen, biceps, muslo)) |>
round(2) peso abdomen biceps muslo
peso 1.00 0.91 0.81 0.88
abdomen 0.91 1.00 0.71 0.81
biceps 0.81 0.71 1.00 0.75
muslo 0.88 0.81 0.75 1.00Ahora probemos con regularización lasso:
## mixture = 1 es regresión lasso
modelo_1 <- linear_reg(mixture = 1, penalty = 0) |>
set_engine("glmnet", lambda.min.ratio = 0)
flujo_1 <- workflow() |>
add_model(modelo_1) |>
add_recipe(grasa_receta)flujo_1 <- flujo_1 |> fit(grasa_ent)
modelo_1 <- extract_fit_parsnip(flujo_1)
coefs <- modelo_1 |> pluck("fit") |> tidy() |>
filter(term != "(Intercept)")
ggplot(coefs,
aes(x = lambda, y = estimate, colour = term)) +
geom_line(size = 1.4) + scale_x_log10() +
scale_colour_manual(values = cbb_palette) +
labs(subtitle = "Regularizacion L1")
Y nótese que conforme aumentamos la penalización, algunas variables salen del modelo (sus coeficientes son cero). Por ejemplo, para un valor de \(lambda\) intermedio, obtenemos un modelo simple de la forma:
coefs |> filter(step == 21) |>
select(term, estimate, lambda)# A tibble: 3 × 3
term estimate lambda
<chr> <dbl> <dbl>
1 edad 0.160 0.434
2 estatura -0.452 0.434
3 abdomen 6.67 0.434Y nótese que este modelo solo incluye 3 variables.La traza confirma que la regularización lasso, además de encoger coeficientes, saca variables del modelo conforme el valor de regularización aumenta.
La razón de esta diferencia cualitativa entre cómo funciona lasso y ridge se puede entender considerando que los problemas de penalización mostrados arriba puede escribirse en forma de problemas de restricción. Por ejemplo, lasso se puede reescribir como el problema de resolver
\[\min_a D(\beta)\] sujeto a \[\sum_{i=1}^p |\beta_j| < t\] En la gráfica siguiente (tomada de Hastie, Tibshirani, y Friedman (2017)), lasso está a la izquierda y ridge está a la derecha, las curvas rojas son curvas de nivel de la suma de cuadrados \(D(a)\), y \(\hat{\beta}\) es el estimador usual de mínimos cuadrados de los coeficientes (sin penalizar). En azul está la restricción:
Nota: es posible también utilizar una penalización que mezcla ridge y lasso:
\[\lambda \left (\alpha \sum_j |a_j| + (1-\alpha)\sum_j a_j^2 \right )\]
y \(\alpha\) es un parámetro que podemos afinar:
# elastic net = ridge + lasso
# mixture es alpha y penalty es lambda
modelo_enet <- linear_reg(mixture = 0.5, penalty = 0.05)
# y si queremos afinar:
modelo_enet <- linear_reg(mixture = tune(), penalty = tune())Otra forma de hacer regularización que se utiliza comunmente se basa en el método de minimización que usamos para obtener nuestra función \(\hat{f}\) para hacer predicciones. La idea es utilizar un método iterativo que comience con una \(f_0\) simple, y luego iterar a una nueva \(f_1\) que se adapta mejor a los datos pero no es muy diferente a \(f_0\). En lugar de seguir iterando hasta llegar a un mínimo de \(L(f),\) evaluamos con una muestra de prueba para encontrar un lugar apropiado para detenernos (early stopping). También podemos modificar \(L(f)\) en cada paso para evitar atorarnos en un mínimo sobreajustado.
Una manera de hacer esto es usando el método de descenso estocástico, (ver apéndices Apéndice A y Apéndice B) que consiste en:
Aunque este método es más útil en casos como redes neuronales o métodos basados en árboles, podemos comenzar por un ejemplo en regresión lineal para entender su efecto regularizador:
library(keras3)
x_grasa <- grasa_receta |> prep() |> juice() |>
select(abdomen, edad, antebrazo, biceps, estatura, muslo, pecho, peso) |>
as.matrix()
vars_nombres <- colnames(x_grasa)
y_grasa <- grasa_receta |> prep() |> juice() |> pull(grasacorp)
## keras tiene distintos algos de optimización
modelo_reg <- keras_model_sequential() |>
layer_dense(units = 1,
kernel_initializer = initializer_constant(0))
modelo_reg |> compile(
loss = "mse",
optimizer = optimizer_sgd(learning_rate = 0.01)
)
# esto es más eficiente hacerlo con callbacks en general:
pesos_tbl <- map_dfr(1:400, function(epoca){
modelo_reg |> fit(
x = x_grasa, y = y_grasa,
epochs = 1,
verbose = FALSE)
pesos_tbl <- get_weights(modelo_reg)[[1]] |> t() |>
as_tibble()
names(pesos_tbl) <- vars_nombres
pesos_tbl |> mutate(epoca = epoca)
}
)library(patchwork)
g_dest <- pesos_tbl |> pivot_longer(cols = -contains(c("grasacorp", "epoca"))) |>
ggplot(aes(x = epoca, y = value, colour = name)) +
geom_line(linewidth = 1.1) + scale_colour_manual(values = cbb_palette) +
scale_x_continuous(trans = compose_trans("log10", "reverse")) +
labs(subtitle = "Descenso estocástico")
g_dest + g_l2
Y vemos que si inicializamos el proceso de minimización con valores chicos, pararnos en una época (iteración completa sobre los datos) nos permite tener un efecto similar al de utilizar regularización tipo L2.
El método de descenso estocástico (usualmente por minilotes) nos permite resolver problemas de optimización, y muchas veces actúa también como regularizador (al cambiar en cada paso la función de pérdida, y utilizando early stopping). Sus ventajas son:
Al cambiar la función de pérdida en cada paso, es posible escapar de puntos estacionarios subóptimos (si el problema tiene varios puntos estacionarios, es decir, donde el gradiente es cero). Evitamos mínimos locales sobreajustados: estos desaparecen cuando consideramos “lotes” o “batches” de datos, en lugar del conjunto completo para cada iteración.
Generalmente detenemos las iteraciones cuando el error de validación deja de disminuir, lo cual es una forma de regularización que mantiene los parámetros en valores relativamente chicos.
Es eficiente en el sentido de que no es necesario utilizar todo los datos para hacer un paso suficientemente bueno, y es escalable a grandes conjuntos de datos.
Es crucial escoger un tamaño de paso adecuado para cada problema. Generalmente se considera un parámetro que debe ser afinado, de manera similar al parámetro de regularización L2 que vimos arriba.
Finalmente, notamos que este tipo de regularización resulta en comportamientos a veces inesperados del desempeño predictivo en función del tamaño del modelo utilizado, por ejemplo, puede ocurrir el fenómeno de doble descenso (ver James et al. (2014)):
Esto se puede deber, por ejemplo, pues cuando existen pocos grados de libertad, la solución es escencialmente única (del paso 1 a 2 de arriba). Sin embargo, con modelos más grandes aparecen muchas otras soluciones posibles: las mejores en entrenamiento son peores en generalización, pero las que son “regularizadas” por descenso estocástico son mejores en desempeño predictivo. Puedes pensar en un polinomio de grado muy alto: la solución de grado alto que encontramos con descenso estocástico puede ser muy suave, en contraste con lo que encontraríamos si llegáramos a un mínimo global o cercano al global.