Introducción

En las últimas semanas he investigado un poco sobre los métodos de regularización en modelos de regresión. Estos métodos mejoran la capacidad de un modelo de regresión para generalizar y evitar el sobreajuste. Además, su costo computacional es mucho menor en comparación con otros modelos de ML.

Dos técnicas comunes de regularización son la regresión Lasso y Ridge:

  • Lasso:
    • Añade un término de penalización proporcional al valor absoluto de las coeficientes del modelo.
    • Promueve que algunos coeficientes se conviertan exactamente cero.
    • Elimina eficazmente las variables penalizadas.
    • Puede ser un método automatizado de selección de variables.
  • Ridge:
    • Añade un término de penalización proporcional al cuadrado de los coeficientes del modelo.
    • No reduce algunos coeficientes a cero como la regresión Lasso.
    • Puede reducir eficazmente el sobreajuste.

Como todo metodo tienen sus ventajas y desventajas, las principal ventaja es que al penalizar los coeficientes ambos metodos evitan el sobre ajuste. Sin embargo, esto mismo genera que no puedan ser utilizados para realizar inferencias si no para tener una mejor capacidad predictiva.

La regresión Lasso es preferible cuando se tiene un número limitado de variables predictoras que son significativas dentro de un modelo de regresión. Esto se debe a que la regresión Lasso puede eliminar de manera efectiva las variables no significativas del modelo. Por otro lado, cuando se tiene un gran número de variables significativas en el modelo, es preferible utilizar el método Ridge.

Otro modelo existente es la Regresión Lineal Elástica la cual utiliza las penalizaciones del método Lasso y Ridge para regularizar el modelo de regresión.

A continuación, se emplearan las tres técnicas en dataset que previamente ya fue analizado y al cual se le ajusto un modelo de Bosque Aleatorio. Los datos corresponde al costo de un seguro de salud para un población y el objetivo es desarrollar un modelo de aprendizaje supervisado para poder predecir dicho costo.

Análisis exploratorio

Estos datos ya fueron previamente explorados y se encontro una relación importante entre el costo del seguro y si la persona en cuestión era fumador activo o no. A continuación, se visualizan algunas variables para demostrarlo:

select(.data = data, charges) %>% 
  mutate(charges_log = log(charges)) %>% 
  pivot_longer(cols = charges:charges_log) %>% 
  ggplot(aes(x = value, fill = name)) +
  geom_density(alpha = 1/2) +
  facet_wrap(vars(name), scales = "free") +
  labs(title = "Distribución de Densidad Costo del Seguro") +
  tidyquant::theme_tq() +
  theme(legend.title = element_blank())

Se puede observar que la distribución del costo del seguro tiene un sesgo positivo o hacia la derecha. Para este tipo de escenarios, es común aplicar la transformación logarítmica para ajustar los datos a una distribución normal.

mutate_if(.tbl = data,
          .predicate = is.character,
          .funs = as_factor) %>% 
  mutate(children_status = if_else(children == 0, "Not Have", "Have") %>% as_factor) %>% 
  ggplot(aes(x = bmi,
             y = charges, 
             colour = smoker)) + 
  geom_point(alpha = 1/3, size = 2) +
  geom_smooth(method = "lm", se = F, linewidth = .7) +
  scale_y_continuous(labels = dollar) +
  labs(x = "Índice de Masa Corporal", y = "Costos de Seguros", color = "¿Es Fumador?") +
  facet_wrap(vars(region), scales = "fixed") +
  tidyquant::theme_tq()

mutate_if(.tbl = data,
          .predicate = is.character,
          .funs = as_factor) %>% 
  mutate(children_status = if_else(children == 0, "Not Have", "Have") %>% as_factor) %>% 
  ggplot(aes(x = age,
             y = charges, 
             colour = smoker)) + 
  geom_point(alpha = 1/3, size = 2) +
  geom_smooth(method = "lm", se = F, linewidth = .7) +
  scale_y_continuous(labels = dollar) +
  labs(x = "Edad del Asegurado", y = "Costos de Seguros", color = "¿Es Fumador?") +
  facet_wrap(vars(region), scales = "fixed") +
  tidyquant::theme_tq()

Partición de los datos

Entrenamiento y prueba

Se particionan los datos en dos conjuntos. El primero para ser utilizado en el entrenamiento de los modelos y el segundo para probar el resultado de dichos modelos.

set.seed(1234)
split <- 
  mutate_at(.tbl = data, 
            vars(charges),
            .funs = log) %>% 
  initial_split(prop = 0.80, 
                strata = charges)

test <- testing(split)
train <- training(split)

Los modelos que se utilizarán contemplan hiperparámetros que serán ajustados por ende. También es necesario crear una validación cruzada del conjunto de entrenamiento para este proceso.

set.seed(1234)
folds <- vfold_cv(data = train, strata = smoker)

Ingeniería de variables

En el análisis previo se diseño un recipe de preprocesamiento de los datos que funciono bastante bien.

recipe <- 
  recipe(charges ~., data = train) %>% 
  update_role(index, new_role = "ID") %>% 
  step_mutate(children_status = if_else(children == 0, "Not Have", "Have"),
              is_obsessive = if_else(bmi > 30, "Yes", "No")) %>% 
  step_string2factor(all_nominal_predictors()) %>%
  step_dummy(all_nominal_predictors()) %>% 
  step_interact(terms = ~ smoker_yes:bmi) %>% 
  step_interact(terms = ~ smoker_yes:age) %>% 
  step_normalize(all_numeric_predictors()) 

Especificación de modelos

Ya teniendo la receta con el preprocesamiento de los datos, sólo falta establecer el flujo de trabajo de los modelos que se utilizarán. Como previamente fue mencionado se crearán worflows para modelos Lasso, Ridge y Elactic Net, adicionalmente a modo de comparativo se utilizarán también un arbol de decisión y un modelo de k vecinos.

ridge_spec <- 
  linear_reg(penalty = tune(), mixture = 0) %>% 
  set_mode("regression") %>% 
  set_engine("glmnet")

lasso_spec <- 
  linear_reg(penalty = tune(), mixture = 1) %>% 
  set_mode("regression") %>% 
  set_engine("glmnet") 

glmnet_spec <- 
  linear_reg(penalty = tune(), 
             mixture = tune()) %>% 
  set_engine("glmnet") %>% 
  set_mode("regression")

knn_spec <- 
  nearest_neighbor(neighbors = tune(), 
                   weight_func = tune(), 
                   dist_power = tune()) %>% 
  set_engine("kknn") %>% 
  set_mode("regression")

cart_spec <- 
  decision_tree(cost_complexity = tune(),
                min_n = tune(), 
                tree_depth = tune()) %>% 
  set_engine("rpart") %>% 
  set_mode("regression")

Ajuste de hiperparámetros

A continuación, se realiza el ajuste de los hiperparámetros para cada uno de los modelos. Aun cuando el dataset es bastante pequeño, tenemos 5 modelos y un total de 10 conjuntos de validación cruzada para este proceso. Esto implica que si deseamos probar 50 combinaciones de hiperparámetros, se deben ajustar 5 * 10 * 50 = 2500 modelos distintos. Por esto, se utilizará cómputo paralelizado para acelerar el proceso.

workflows <- 
  workflow_set(
    preproc = list(recipe = recipe), 
    models = list(ridge = ridge_spec, lasso = lasso_spec, glmnet = glmnet_spec, knn = knn_spec,
                  tree = cart_spec))
    
grid_ctrl <-
  control_grid(
    save_pred = TRUE,
    parallel_over = "everything",
    save_workflow = TRUE
  )
cores <- parallel::detectCores(logical = FALSE)
cl <- parallel::makePSOCKcluster(cores)
doParallel::registerDoParallel(cl)

grid_results <-
  workflows %>%
  workflow_map(
    fn = "tune_grid",
    seed = 1503,
    resamples = folds,
    grid = 50,
    control = grid_ctrl,
    metrics = metric_set(rmse, rsq)
  )

doParallel::stopImplicitCluster()

Resultado del ajuste de los hiperparámetros:

grid_results %>% 
  rank_results(rank_metric = "rmse") %>% 
  ggplot(aes(x = rank, y = mean, color = wflow_id)) +
  geom_errorbar(aes(
    ymin = mean - std_err,
    ymax = mean + std_err),
    alpha = 0.5) +
  geom_point(size = 0.5) +
  facet_wrap(vars(.metric), scales = "free_y", ncol = 2) +
  tidyquant::theme_tq() +
  theme(legend.position = "top",
        legend.text = element_text(size = 10),
        text = element_text(color = "black"))

Se puede observar que el modelo de regresión elástica disminuye la raíz del error cuadrático medio y maximiza el coeficiente de determinación, por lo tanto, es el mejor modelo de todos los ajustes. Empleado las funciones extract_workflow_set_result y select_best se puede extraer la mejor combinación de hiperparámetros.

extract_workflow_set_result(
  x = grid_results, 
  id = "recipe_glmnet") %>%
  select_best(metric = "rmse") %>% 
  select(-.config) %>%
  gt::gt(caption = "Mejor Combinación de Hiperparámetros")
Mejor Combinación de Hiperparámetros
penalty mixture
2.314584e-07 0.08361681

Validación en conjunto de prueba

Por ultimó, solo queda finalizar el mejor modelo con la mejor combinación de hiperparámetros y emplear dicho modelo para predecir el conjunto de prueba y verificar el rendimiento del modelo en datos nunca antes vistos.

best_results <- 
  grid_results %>% 
  extract_workflow_set_result("recipe_glmnet") %>% 
  select_best(metric = "rmse") 

final_model <- 
  grid_results %>% 
  extract_workflow("recipe_glmnet") %>% 
  finalize_workflow(best_results) %>% 
  last_fit(split = split)
collect_metrics(x = final_model) %>% 
  select(-all_of(c(".config", ".estimator"))) %>% 
  mutate_at(vars(.estimate), ~ round(., 3)) %>% 
  gt::gt(caption = "Métricas") %>% 
  gt::cols_label(.metric = "Métricas", .estimate = "Valor")
Métricas
Métricas Valor
rmse 0.379
rsq 0.834

La forma más sencilla de evaluar visualmente la concordancia entre los valores observados y los predichos es con un gráfico de dispersión.

Se puede observar que se obtuvo un modelo bastante decente empleando los métodos de regularización en modelos de regresión multivariante. Además, es importante destacar que dicha metodología supero ampliamente a modelos de arboles de decisión y k vecinos.