MACHINE LEARNING: H2O, IML y DALEX


Versión PDF: Github

Más sobre ciencia de datos: cienciadedatos.net

Introducción


En este documento se presenta un ejemplo práctico de cómo entrenar modelos de machine learning con la librería H2O y de cómo compararlos e interpretarlos con Dalex e IML.

H2O


H2O es un producto creado por la compañía H2O.ai con el objetivo de combinar los principales algoritmos de machine learning con Big Data. Gracias a su forma de comprimir y almacenar los datos, H2O es capaz de trabajar con grandes volúmenes de registros en un único ordenador (emplea todos sus cores) o en un cluster de muchos ordenadores. Internamente, H2O está escrito en Java y sigue el paradigma Key/Value para almacenar los datos y Map/Reduce para sus algoritmos. Gracias a sus API, es posible acceder a todas sus funciones desde R, Python o Scala, así como por una interfaz web llamada Flow.

Las características que más destacan de esta librería son:

  • Incorpora los principales algoritmos de machine learning:

    • Cox Proportional Hazards (CoxPH)
    • Deep Learning (Neural Networks)
    • Distributed Random Forest (DRF)
    • Generalized Linear Model (GLM)
    • Gradient Boosting Machine (GBM)
    • Naïve Bayes Classifier
    • Stacked Ensembles
    • XGBoost
    • Aggregator
    • Generalized Low Rank Models (GLRM)
    • K-Means Clustering
    • Isolation Forest
    • Principal Component Analysis (PCA)
    • Word2vec
  • Sus algoritmos permites paralelización para hacer uso de todos los cores del ordenador o incluso de un cluster.

  • Incorpora en los propios modelos las principales trasformaciones de preprocesado (escalado, encoding de variables cualitativas, eliminación de predictores con varianza constante e imputación de valores ausentes). De esta forma, todas las transformaciones se aprenden durante el entrenamiento y se aplican automáticamente a las nuevas predicciones.

  • Permite la búsqueda de hiperparámetros por grid search y random search.

  • Todos los algoritmos incluyen criterios de parada temprana para agilizar el entrenamiento.

Todas estas características hacen de H2O una herramientas excelente aun cuando el volumen de datos es limitado. Para conocer más detalles sobre esta librería y sus modelos consultar Machine Learning con H2O y R.

Dalex y IML


En términos generales, los modelos de machine learning que consiguen mejores resultados en las predicciones, lo hacen gracias a su capacidad para encontrar relaciones complejas en los datos (interacciones entre predictores, relaciones no lineales…). Sin embargo, una de las desventajas asociadas es que su interpretabilidad suele ser baja. No es fácil comprender cómo está participando cada predictor en el modelo y en sus predicciones.

A medida que ha avanzado el desarrollo de modelos predictivos, se han ido mejorando las estrategias para entender su comportamiento. Algunas de ellas son intrínsecas al algoritmo (los coeficientes de regresión de un modelo lineal, la importancia de los predictores en un random forest…) y otras son agnósticas al tipo de modelo. Los paquetes Dalex e IML tienen implementados la mayoría de estos métodos y son compatibles con modelos de H2O.

Otros paquetes


Paquetes utilizados a lo largo del documento.

library(tidyverse)
library(h2o)
library(DALEX)
library(DALEXtra)
library(iml)
library(skimr)
library(DataExplorer)
library(ggpubr)
library(univariateML)
library(recipes)



Datos


El set de datos rent, disponible en el paquete gamlss.data, contiene información sobre el precio del alquiler de 1969 viviendas situadas en Munich en el año 1993. Además del precio, incluye 9 variables adicionales:

  • R: precio del alquiler.

  • Fl: metros cuadrados de la vivienda.

  • A: año de construcción.

  • B: si tiene cuarto de baño (1) o no (0).

  • H: si tiene calefacción central (1) o no (0).

  • L: si la cocina está equipada (1) o no (0).

  • Sp: si la calidad del barrio donde está situada la vivienda es superior la media (1) o no (0).

  • Sm: si la calidad del barrio donde está situada la vivienda es inferior la media (1) o no (0).

  • loc: combinación de Sp y Sm indicando si la calidad del barrio donde está situada la vivienda es inferior (1), igual (2) o superior (3) a la media.

data("rent", package = "gamlss.data")
datos <- rent

# Se descartan las variables SP y Sm ya que su información está recogida en la
# variable loc.
datos <- datos %>% select(-c(Sp, Sm))

# Se renombran las columnas para que sean más descriptivas.
colnames(datos) <- c("precio", "metros", "anyo", "banyo", "calefaccion", "cocina",
                     "situacion")



Análisis exploratorio


Antes de entrenar un modelo predictivo, o incluso antes de realizar cualquier cálculo con un nuevo conjunto de datos, es muy importante realizar una exploración descriptiva de los mismos. Este proceso permite entender mejor qué información contiene cada variable, así como detectar posibles errores.

Los paquetes skimr, DataExplorer y GGally facilitan mucho esta tarea gracias a sus funciones preconfiguradas.

Tabla resumen

skim(datos)
Data summary
Name datos
Number of rows 1969
Number of columns 7
_______________________
Column type frequency:
factor 4
numeric 3
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
banyo 0 1 FALSE 2 0: 1925, 1: 44
calefaccion 0 1 FALSE 2 0: 1580, 1: 389
cocina 0 1 FALSE 2 0: 1808, 1: 161
situacion 0 1 FALSE 3 2: 1247, 3: 550, 1: 172

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
precio 0 1 811.88 379.00 101.7 544.2 737.8 1022 3000 ▇▇▂▁▁
metros 0 1 67.73 20.86 30.0 52.0 67.0 82 120 ▆▇▇▅▂
anyo 0 1 1948.48 29.02 1890.0 1934.0 1957.0 1972 1988 ▃▁▃▇▇
head(datos)

Todas las columnas tienen el tipo adecuado.

Valores ausentes


Junto con el estudio del tipo de variables, es básico conocer el número de observaciones disponibles y si todas ellas están completas. Los valores ausentes son importantes a la hora de crear modelos, algunos algoritmos no aceptan observaciones incompletas o bien se ven muy influenciados por ellas.

# Número de datos ausentes por variable
datos %>% map_dbl(.f = function(x){sum(is.na(x))})
##      precio      metros        anyo       banyo calefaccion      cocina 
##           0           0           0           0           0           0 
##   situacion 
##           0
plot_missing(
  data    = datos, 
  title   = "Porcentaje de valores ausentes",
  ggtheme = theme_bw(),
  theme_config = list(legend.position = "none")
)



Variables respuesta


Cuando se crea un modelo, conviene estudiar la distribución de la variable respuesta, ya que, a fin de cuentas, es lo que interesa predecir. La variable precio tiene una distribución asimétrica con una cola positiva debido a que unas pocas viviendas tienen un precio de alquiler muy superior a la media. Este tipo de distribución suele visualizarse mejor tras aplicar el logarítmica o la raíz cuadrada.

p1 <- ggplot(data = datos, aes(x = precio)) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Distribución original", x = "gastos_totales") +
      theme_bw() 

p2 <- ggplot(data = datos, aes(x = sqrt(precio))) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Transformación raíz cuadrada", x = "gastos_totales") +
      theme_bw() 

p3 <- ggplot(data = datos, aes(x = log10(precio))) +
      geom_density(fill = "steelblue", alpha = 0.8) +
      geom_rug(alpha = 0.1) +
      scale_x_continuous(labels = scales::comma) +
      labs(title = "Transformación logarítmica") +
      theme_bw()

ggarrange(p1, p2, p3, ncol = 1, align = "v")

# Tabla de estadísticos principales 
summary(datos$precio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   101.7   544.2   737.8   811.9  1022.0  3000.0



Algunos modelos de machine learning y aprendizaje estadístico requieren que la variable respuesta se distribuya de una forma determinada. Por ejemplo, para los modelos de regresión lineal (LM), la distribución tiene que ser de tipo normal. Para los modelos lineales generalizados (GLM), la distribución tiene que ser de la familia exponencial. Existen varios paquetes en R que permiten identificar a qué distribución se ajustan mejor los datos, uno de ellos es univariateML. Para conocer más sobre cómo identificar distribuciones consultar Ajuste de distribuciones con R.

# Se comparan únicamente las distribuciones con un dominio [0, +inf)
# Cuanto menor el valor AIC mejor el ajuste
comparacion_aic <- AIC(
                    mlbetapr(datos$precio),
                    mlexp(datos$precio),
                    mlinvgamma(datos$precio),
                    mlgamma(datos$precio),
                    mllnorm(datos$precio),
                    mlrayleigh(datos$precio),
                    mlinvgauss(datos$precio),
                    mlweibull(datos$precio),
                    mlinvweibull(datos$precio),
                    mllgamma(datos$precio)
                   )
comparacion_aic %>% rownames_to_column(var = "distribucion") %>% arrange(AIC)

La distribución con mejor ajuste acorde al valor AIC es la gamma.

Variables continuas


plot_density(
  data    = datos %>% select(-precio),
  ncol    = 3,
  title   = "Distribución variables continuas",
  ggtheme = theme_bw(),
  theme_config = list(
                  plot.title = element_text(size = 16, face = "bold"),
                  strip.text = element_text(colour = "black", size = 12, face = 2)
                 )
  )



Como el objetivo del estudio es predecir el precio de alquiler de las viviendas, el análisis de cada variable se hace también en relación a la variable respuesta precio. Analizando los datos de esta forma, se pueden extraer ideas sobre qué variables están más relacionadas con el precio y de qué forma.

custom_corr_plot <- function(variable1, variable2, df, alpha=0.3){
  p <- df %>%
       mutate(
         # Truco para que se ponga el título estilo facet
        title = paste(toupper(variable2), "vs", toupper(variable1))
       ) %>%
       ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) + 
       geom_point(alpha = alpha) +
       # Tendencia no lineal
       geom_smooth(se = FALSE, method = "gam", formula =  y ~ splines::bs(x, 3)) +
       # Tendencia lineal
       geom_smooth(se = FALSE, method = "lm", color = "firebrick") +
       facet_grid(. ~ title) +
       theme_bw() +
       theme(strip.text = element_text(colour = "black", size = 10, face = 2),
             axis.title = element_blank())
  return(p)
}
variables_continuas <- c("anyo", "metros")

plots <- map(
            .x = variables_continuas,
            .f = custom_corr_plot,
            variable2 = "precio",
            df = datos
         )

ggarrange(plotlist = plots, ncol = 2, nrow = 1) %>%
  annotate_figure(
    top = text_grob("Correlación con el precio de alquiler", face = "bold", size = 16,
                    x = 0.4)
  )



Correlación variables continuas


Algunos modelos (LM, GLM, …) se ven perjudicados si incorporan predictores altamente correlacionados. Por esta razón, es conveniente estudiar el grado de correlación entre las variables disponibles.

plot_correlation(
  data     = datos,
  cor_args = list(use = "pairwise.complete.obs"),
  type     = "continuous",
  title    = "Matriz de correlación variables continuas",
  theme_config = list(legend.position = "none",
                      plot.title = element_text(size = 16, face = "bold"),
                      axis.title = element_blank(),
                      axis.text.x = element_text(angle = -45, hjust = +0.1)
                     )
)

GGally::ggscatmat(
  data = datos %>% select_if(is.numeric),
  alpha = 0.1) +
theme_bw() +
labs(title = "Correlación por pares") +
theme(
  plot.title = element_text(size = 16, face = "bold"),
  #axis.text = element_blank(),
  strip.text = element_text(colour = "black", size = 10, face = 2)
)



Variables cualitativas


plot_bar(
  datos,
  ncol    = 2,
  title   = "Número de observaciones por grupo",
  ggtheme = theme_bw(),
  theme_config = list(
                   plot.title = element_text(size = 16, face = "bold"),
                   strip.text = element_text(colour = "black", size = 12, face = 2),
                   legend.position = "none"
                  )
)

custom_box_plot <- function(variable1, variable2, df, alpha=0.3){
  p <- df %>%
       mutate(
         # Truco para que se ponga el título estilo facet
        title = paste(toupper(variable2), "vs", toupper(variable1))
       ) %>%
       ggplot(aes(x = !!sym(variable1), y = !!sym(variable2))) + 
       geom_violin(alpha = alpha) +
       geom_boxplot(width = 0.1, outlier.shape = NA) +
       facet_grid(. ~ title) +
       theme_bw() +
       theme(strip.text = element_text(colour = "black", size = 10, face = 2),
             axis.title = element_blank())
  return(p)
}
variables_cualitativas <- c("banyo", "calefaccion", "cocina", "situacion")

plots <- map(
            .x = variables_cualitativas,
            .f = custom_box_plot,
            variable2 = "precio",
            df = datos
         )

ggarrange(plotlist = plots, ncol = 2, nrow = 2) %>%
  annotate_figure(
    top = text_grob("Correlación con precio", face = "bold", size = 16,
                    x = 0.28)
  )

Si alguno de los niveles de una variable cualitativa tiene muy pocas observaciones en comparación a los otros niveles, puede ocurrir que, durante la validación cruzada o bootstrapping, algunas particiones no contengan ninguna observación de dicha clase (varianza cero), lo que puede dar lugar a errores. En este caso hay que tener precaución con la variable banyo.

table(datos$banyo)
## 
##    0    1 
## 1925   44



Modelos


El objetivo es crear un modelo capaz de predecir el precio del alquiler. En los siguientes apartados se emplean los principales algoritmos disponibles en H2O y se comparan para identificar el que mejor resultados consigue.

Iniciación del cluster


# Creación de un cluster local con todos los cores disponibles.
h2o.init(ip = "localhost",
         # -1 indica que se empleen todos los cores disponibles.
         nthreads = -1,
         # Máxima memoria disponible para el cluster.
         max_mem_size = "6g")

# Se eliminan los datos del cluster por si ya había sido iniciado.
h2o.removeAll()
h2o.no_progress()
datos_h2o <- as.h2o(x = datos, destination_frame = "datos_h2o")



Train-Test


set.seed(123)
particiones     <- h2o.splitFrame(data = datos_h2o, ratios = c(0.8), seed = 123)
datos_train_h2o <- h2o.assign(data = particiones[[1]], key = "datos_train_h2o")
datos_test_h2o  <- h2o.assign(data = particiones[[2]], key = "datos_test_h2o")

Se almacenan en formato data.frame para las funciones de diagnóstico de dalexy iml.

datos_train <- as.data.frame(datos_train_h2o)
datos_test  <- as.data.frame(datos_test_h2o)

Se comprueba que la variable respuesta en similar en los dos conjuntos y que también la variable banyos.

summary(datos_train$precio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   101.7   544.7   741.6   813.9  1026.5  3000.0
summary(datos_test$precio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   127.1   537.8   718.6   804.0  1000.0  2764.5
summary(datos_train$banyo)/nrow(datos_train)
##          0          1 
## 0.97963081 0.02036919
summary(datos_test$banyo)/nrow(datos_test)
##          0          1 
## 0.96984925 0.03015075



GLM


Entrenamiento


# Se define la variable respuesta y los predictores.
var_respuesta <- "precio"
# Para este modelo se emplean todos los predictores disponibles.
predictores   <- setdiff(h2o.colnames(datos_train_h2o), var_respuesta)



Optimización de hiperparámetros

# Valores de alpha que se van a comparar.
hiperparametros <- list(alpha = seq(0, 1, 0.1))
grid_glm <- h2o.grid(
    # Algoritmo y parámetros
    algorithm      = "glm",
    family         = "gamma",
    # Variable respuesta y predictores
    y              = var_respuesta,
    x              = predictores, 
    # Datos de entrenamiento
    training_frame = datos_train_h2o,
    # Preprocesado
    standardize    = TRUE,
    missing_values_handling = "Skip",
    ignore_const_cols = TRUE,
    # Hiperparámetros
    hyper_params    = hiperparametros,
    # Tipo de búsqueda
    search_criteria = list(strategy = "Cartesian"),
    lambda_search   = TRUE,
    # Selección automática del solver adecuado
    solver          = "AUTO",
    # Estrategia de validación para seleccionar el mejor modelo
    seed            = 123,
    nfolds          = 3,
    keep_cross_validation_predictions = TRUE,
    grid_id         = "grid_glm"
)

# Se muestran los modelos ordenados de mayor a menor rmse
resultados_grid_glm <- h2o.getGrid(
                         grid_id = "grid_glm",
                         sort_by = "rmse",
                         decreasing = FALSE
                       )

as.data.frame(resultados_grid_glm@summary_table)



Mejor modelo

# Se reentrena el modelo con los mejores hiperparámetros
mejores_hiperparam <- h2o.getModel(resultados_grid_glm@model_ids[[1]])@parameters

modelo_glm <- h2o.glm(
    # Variable respuesta y predictores
    y              = var_respuesta,
    x              = predictores,
    # Datos de entrenamiento
    training_frame = datos_train_h2o,
    # Preprocesado
    standardize    = TRUE,
    missing_values_handling = "Skip",
    ignore_const_cols = TRUE,
    # Hiperparámetros
    alpha = mejores_hiperparam$alpha,
    lambda_search   = TRUE,
    # Selección automática del solver adecuado
    solver          = "AUTO",
    # Estrategia de validación (para comparar con otros modelos)
    seed            = 123,
    nfolds          = 10,
    keep_cross_validation_predictions = TRUE,
    model_id        = "modelo_glm"
)



Influencia predictores


Importancia predictores

h2o.varimp(modelo_glm)
h2o.varimp_plot(modelo_glm)



Diagnóstico de los residuos


explainer_glm <- DALEXtra::explain_h2o(
                  model = modelo_glm,
                  data  = datos_train[, predictores],
                  y     = datos_train[[var_respuesta]],
                  label = "modelo_glm"
                )

plot(model_performance(explainer_glm))

predicciones_train <- h2o.predict(
                        modelo_glm,
                        newdata = datos_train_h2o
                      )
predicciones_train <- h2o.cbind(
                       datos_train_h2o["precio"],
                       predicciones_train
                      )
predicciones_train <- as.data.frame(predicciones_train)
predicciones_train <- predicciones_train %>%
                      mutate(
                        residuo = predict - precio
                      )

p1 <- ggplot(predicciones_train, aes(x = precio, y  = predict)) +
      geom_point(alpha = 0.1) +
      geom_smooth(method = "gam", color = "red", size = 1, se = FALSE) +
      labs(title = "Predicciones vs valor real") +
      theme_bw()

p2 <- ggplot(predicciones_train, aes(1:nrow(predicciones_train), y  = residuo)) +
      geom_point(alpha = 0.1) +
      geom_hline(yintercept = 0, color = "red", size = 1) +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(predicciones_train, aes(x  = residuo)) +
      geom_density() +
      geom_rug() +
      labs(title = "Residuos del modelo") +
      theme_bw()

p4 <- ggplot(predicciones_train, aes(sample  = predict)) +
      stat_qq() +
      stat_qq_line(color = "red", size = 1) +
      labs(title = "QQ-plot residuos del modelo") +
      theme_bw()

ggpubr::ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) %>%
ggpubr::annotate_figure(
  top = ggpubr::text_grob("Diagnóstico residuos entrenamiento",
                          color = "black", face = "bold", size = 14)
)



Predicción test


Predicciones

predicciones_test <- h2o.predict(
                        object  = modelo_glm,
                        newdata = datos_test_h2o
                     )
predicciones_test     <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_test



Error test

rmse_test_glm <- MLmetrics::RMSE(
                    y_pred = datos_test$precio,
                    y_true = datos_test$prediccion
                 )
paste("Error de test (rmse) del modelo GLM:", rmse_test_glm)
## [1] "Error de test (rmse) del modelo GLM: 378.70675549663"



Escritura del modelo


# Se guarda el modelo en el directorio actual en la carpeta modelos
h2o.saveModel(object = modelo_glm, path = "./modelos", force = TRUE)
file.rename(file.path("./modelos", modelo_glm@model_id), "./modelos/modelo_glm.h2o")



GBM


Entrenamiento


Optimización de hiperparámetros

# Hiperparámetros que se quieren comparar (random search)
hiperparametros <- list(
                     ntrees      = c(500, 1000, 2000),
                     learn_rate  = seq(0.01, 0.1, 0.01),
                     max_depth   = seq(2, 15, 1),
                     sample_rate = seq(0.5, 1.0, 0.2),
                     col_sample_rate = seq(0.1, 1.0, 0.3)
                    )
# Criterios de parada para la búsqueda
search_criteria <- list(
                    strategy = "RandomDiscrete",
                    max_models = 50, # Número máximo de combinaciones
                    max_runtime_secs = 60*10, # Tiempo máximo de búsqueda
                    stopping_tolerance = 0.001, # Mejora mínima
                    stopping_rounds = 5,
                    seed = 123
                   )


grid_gbm <- h2o.grid(
    # Algoritmo y parámetros
    algorithm = "gbm",
    distribution = "gaussian",
    # Variable respuesta y predictores
    y = var_respuesta,
    x = predictores,
    # Datos de entrenamiento
    training_frame = datos_train_h2o,
    # Preprocesado
    ignore_const_cols = TRUE,
    # Parada temprana
    score_tree_interval = 100,
    stopping_rounds = 3,
    stopping_metric = "rmse",
    stopping_tolerance = 0.01,
    # Hiperparámetros optimizados
    hyper_params = hiperparametros,
    # Estrategia de validación para seleccionar el mejor modelo
    seed   = 123,
    nfolds = 3,
    # Tipo de búsqueda
    search_criteria = search_criteria,
    keep_cross_validation_predictions = TRUE,
    grid_id         = "grid_gbm"
)

# Se muestran los modelos ordenados de mayor a menor rmse
resultados_grid_gbm <- h2o.getGrid(
                        grid_id = "grid_gbm",
                        sort_by = "rmse",
                        decreasing = FALSE
                       )

as.data.frame(resultados_grid_gbm@summary_table)



Mejor modelo

# Se reentrena el modelo con los mejores hiperparámetros
mejores_hiperparam <- h2o.getModel(resultados_grid_gbm@model_ids[[1]])@parameters

modelo_gbm <- h2o.gbm(
    # Variable respuesta y predictores
    y = var_respuesta,
    x = predictores,
    # Datos de entrenamiento
    training_frame = datos_train_h2o,
    # Preprocesado
    ignore_const_cols = TRUE,
    # Hiperparámetros
    learn_rate  = mejores_hiperparam$learn_rate,
    max_depth   = mejores_hiperparam$max_depth,
    ntrees      =  mejores_hiperparam$ntrees,
    sample_rate = mejores_hiperparam$sample_rate,
    # Estrategia de validación (para comparar modelos)
    seed            = 123,
    nfolds          = 10,
    keep_cross_validation_predictions = TRUE,
    model_id        = "modelo_gbm"
)



Influencia predictores


Importancia predictores

h2o.varimp(modelo_gbm)
h2o.varimp_plot(modelo_gbm)



Gráficos PDP

par(mfrow = c(3, 2))
pdp_plots <- h2o.partialPlot(
                object = modelo_gbm,
                data   = datos_train_h2o,
                cols   = predictores,
                nbins  = 20,
                plot   = TRUE,
                plot_stddev = TRUE
             )

par(mfrow = c(1, 1))



Curvas ICE

custom_ice <- function(modelo_h2o, data, y, predictores = NA) {

  predictor <- Predictor$new(
                model = modelo_h2o, 
                data  = data, 
                y     = y, 
                class = "regression"
              )
  
  if(is.na(predictores)) {
    predictores <- colnames(data)
  }
  
  graficos <- list()
  
  for (i in seq_along(predictores)) {
    ice <- FeatureEffect$new(
                  predictor = predictor,
                  feature   =  predictores[i],
                  method    = "pdp+ice",
                  grid.size = 20
           )
    plot_ice <- plot(ice) + ggtitle(predictores[i])
    graficos[[i]] <- plot_ice
  }
  
  return(graficos)
}

graficos_ice <- custom_ice(
                  modelo_h2o = modelo_gbm,
                  data       = datos_train[, predictores],
                  y          = datos_train[[var_respuesta]]
                )

ggarrange(plotlist = graficos_ice, ncol = 2, nrow = 3)



Diagnóstico de los residuos


explainer_gbm <- DALEXtra::explain_h2o(
                  model = modelo_gbm,
                  data  = datos_train[, predictores],
                  y     = datos_train[[var_respuesta]],
                  label = "modelo_gbm" 
                )
## Preparation of a new explainer is initiated
##   -> model label       :  modelo_gbm 
##   -> data              :  1571  rows  6  cols 
##   -> target variable   :  1571  values 
##   -> model_info        :  package h2o , ver. 3.28.1.2 , task regression (  default  ) 
##   -> predict function  :  yhat.H2ORegressionModel  will be used (  default  )
##   -> predicted values  :  numerical, min =  261.1752 , mean =  814.0984 , max =  1798.123  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -789.5428 , mean =  -0.2181171 , max =  1629.115  
##   A new explainer has been created! 
plot(model_performance(explainer_gbm))

predicciones_train <- h2o.predict(
                        modelo_gbm,
                        newdata = datos_train_h2o
                      )
predicciones_train <- h2o.cbind(
                       datos_train_h2o["precio"],
                       predicciones_train
                      )
predicciones_train <- as.data.frame(predicciones_train)
predicciones_train <- predicciones_train %>%
                      mutate(
                        residuo = predict - precio
                      )

p1 <- ggplot(predicciones_train, aes(x = precio, y  = predict)) +
      geom_point(alpha = 0.1) +
      geom_smooth(method = "gam", color = "red", size = 1, se = FALSE) +
      labs(title = "Predicciones vs valor real") +
      theme_bw()

p2 <- ggplot(predicciones_train, aes(1:nrow(predicciones_train), y  = residuo)) +
      geom_point(alpha = 0.1) +
      geom_hline(yintercept = 0, color = "red", size = 1) +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(predicciones_train, aes(x  = residuo)) +
      geom_density() +
      geom_rug() +
      labs(title = "Residuos del modelo") +
      theme_bw()

p4 <- ggplot(predicciones_train, aes(sample  = predict)) +
      stat_qq() +
      stat_qq_line(color = "red", size = 1) +
      labs(title = "QQ-plot residuos del modelo") +
      theme_bw()

ggpubr::ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) %>%
ggpubr::annotate_figure(
  top = ggpubr::text_grob("Diagnóstico residuos entrenamiento",
                          color = "black", face = "bold", size = 14)
)



Predicción test


Predicciones

predicciones_test <- h2o.predict(
                        object  = modelo_gbm,
                        newdata = datos_test_h2o
                     )
predicciones_test     <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_test

Error test

rmse_test_gbm <- MLmetrics::RMSE(
                  y_pred = datos_test$precio,
                  y_true = datos_test$prediccion
                )
paste("Error de test (rmse) del modelo GBM:", rmse_test_gbm)
## [1] "Error de test (rmse) del modelo GBM: 279.859014376654"



Escritura del modelo


# Se guarda el modelo en el directorio actual en la carpeta modelos
h2o.saveModel(object = modelo_gbm, path = "./modelos", force = TRUE)
file.rename(file.path("./modelos", modelo_gbm@model_id), "./modelos/modelo_gbm.h2o")



RF


Entrenamiento


Optimización de hiperparámetros

# Hiperparámetros que se quieren comparar (random search)
hiperparametros <- list(
                     ntrees      = c(500, 1000, 2000),
                     max_depth   = seq(2, 15, 1),
                     sample_rate = seq(0.5, 1.0, 0.2)
                    )
# Criterios de parada para la búsqueda
search_criteria <- list(
                    strategy = "RandomDiscrete",
                    max_models = 50, # Número máximo de combinaciones
                    max_runtime_secs = 60*10, # Tiempo máximo de búsqueda
                    stopping_tolerance = 0.001, # Mejora mínima
                    stopping_rounds = 5,
                    seed = 123
                   )

grid_rf <- h2o.grid(
    # Algoritmo y parámetros
    algorithm = "drf",
    # Variable respuesta y predictores
    y = var_respuesta,
    x = predictores,
    # Datos de entrenamiento
    training_frame = datos_train_h2o,
    # Preprocesado
    ignore_const_cols = TRUE,
    # Parada temprana
    score_tree_interval = 100,
    stopping_rounds = 5,
    stopping_metric = "rmse",
    stopping_tolerance = 0.01,
    # Hiperparámetros optimizados
    hyper_params = hiperparametros,
    # Estrategia de validación para seleccionar el mejor modelo
    seed   = 123,
    nfolds = 3,
    # Tipo de búsqueda
    search_criteria = search_criteria,
    keep_cross_validation_predictions = TRUE,
    grid_id         = "grid_rf"
)
# Se muestran los modelos ordenados de mayor a menor rsme
resultados_grid_rf <- h2o.getGrid(
                        grid_id = "grid_rf",
                        sort_by = "rmse",
                        decreasing = FALSE
                      )

as.data.frame(resultados_grid_rf@summary_table)



Mejor modelo

# Se reentrena el modelo con los mejores hiperparámetros
mejores_hiperparam <- h2o.getModel(resultados_grid_rf@model_ids[[1]])@parameters

modelo_rf <- h2o.randomForest(
    # Variable respuesta y predictores
    y = var_respuesta,
    x  = predictores,
    # Datos de entrenamiento
    training_frame = datos_train_h2o,
    # Preprocesado
    ignore_const_cols = TRUE,
    # Hiperparámetros
    max_depth   = mejores_hiperparam$max_depth,
    ntrees      =  mejores_hiperparam$ntrees,
    sample_rate = mejores_hiperparam$sample_rate,
    # Estrategia de validación para seleccionar el mejor modelo
    seed            = 123,
    nfolds          = 10,
    keep_cross_validation_predictions = TRUE,
    model_id        = "modelo_rf"
)



Influencia predictores


Importancia predictores

h2o.varimp(modelo_rf)
h2o.varimp_plot(modelo_rf)



Gráficos PDP

par(mfrow = c(3, 2))
pdp_plots <- h2o.partialPlot(
                object = modelo_rf,
                data   = datos_train_h2o,
                cols   = predictores,
                nbins  = 20,
                plot   = TRUE,
                plot_stddev   = TRUE
              )

par(mfrow = c(1, 1))



Curvas ICE

custom_ice <- function(modelo_h2o, data, y, predictores = NA) {

  predictor <- Predictor$new(
                model = modelo_h2o, 
                data  = data, 
                y     = y, 
                class = "regression"
              )
  
  if(is.na(predictores)) {
    predictores <- colnames(data)
  }
  
  graficos <- list()
  
  for (i in seq_along(predictores)) {
    ice <- FeatureEffect$new(
                  predictor = predictor,
                  feature   =  predictores[i],
                  method    = "pdp+ice",
                  grid.size = 20
           )
    plot_ice <- plot(ice) + ggtitle(predictores[i])
    graficos[[i]] <- plot_ice
  }
  
  return(graficos)
}

graficos_ice <- custom_ice(
                  modelo_h2o = modelo_gbm,
                  data       = datos_train[, predictores],
                  y          = datos_train[[var_respuesta]]
                )

ggarrange(plotlist = graficos_ice, ncol = 2, nrow = 3)



Diagnóstico de los residuos


explainer_rf <- DALEXtra::explain_h2o(
                  model = modelo_rf,
                  data  = datos_train[, predictores],
                  y     = datos_train[[var_respuesta]],
                  label = "modelo_rf" 
                )
## Preparation of a new explainer is initiated
##   -> model label       :  modelo_rf 
##   -> data              :  1571  rows  6  cols 
##   -> target variable   :  1571  values 
##   -> model_info        :  package h2o , ver. 3.28.1.2 , task regression (  default  ) 
##   -> predict function  :  yhat.H2ORegressionModel  will be used (  default  )
##   -> predicted values  :  numerical, min =  304.4805 , mean =  813.403 , max =  1722.642  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -797.0642 , mean =  0.477308 , max =  1331.144  
##   A new explainer has been created! 
plot(model_performance(explainer_rf))

predicciones_train <- h2o.predict(
                        modelo_rf,
                        newdata = datos_train_h2o
                      )
predicciones_train <- h2o.cbind(
                       datos_train_h2o["precio"],
                       predicciones_train
                      )
predicciones_train <- as.data.frame(predicciones_train)
predicciones_train <- predicciones_train %>%
                      mutate(
                        residuo = predict - precio
                      )

p1 <- ggplot(predicciones_train, aes(x = precio, y  = predict)) +
      geom_point(alpha = 0.5) +
      geom_smooth(method = "gam", color = "red", size = 1, se = FALSE) +
      labs(title = "Predicciones vs valor real") +
      theme_bw()

p2 <- ggplot(predicciones_train, aes(1:nrow(predicciones_train), y  = residuo)) +
      geom_point(alpha = 0.5) +
      geom_hline(yintercept = 0, color = "red", size = 1) +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(predicciones_train, aes(x  = residuo)) +
      geom_density() +
      geom_rug() +
      labs(title = "Residuos del modelo") +
      theme_bw()

p4 <- ggplot(predicciones_train, aes(sample  = predict)) +
      stat_qq() +
      stat_qq_line(color = "red", size = 1) +
      labs(title = "QQ-plot residuos del modelo") +
      theme_bw()

ggpubr::ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) %>%
ggpubr::annotate_figure(
  top = ggpubr::text_grob("Diagnóstico residuos entrenamiento",
                          color = "black", face = "bold", size = 14)
)



Predicción test


Predicciones

predicciones_test <- h2o.predict(
                        object  = modelo_rf,
                        newdata = datos_test_h2o
                     )
predicciones_test <- as.vector(predicciones_test)

datos_test$prediccion <- predicciones_test

Error test

rmse_test_rf <- MLmetrics::RMSE(
                  y_pred = datos_test$precio,
                  y_true = datos_test$prediccion
                )
paste("Error de test (rmse) del modelo Random Forest:", rmse_test_rf)
## [1] "Error de test (rmse) del modelo Random Forest: 285.701289143036"



Escritura del modelo


# Se guarda el modelo en el directorio actual en la carpeta modelos
h2o.saveModel(object = modelo_rf, path = "./modelos", force = TRUE)
file.rename(file.path("./modelos", modelo_rf@model_id), "./modelos/modelo_rf.h2o")



XGBOOST


Entrenamiento


Optimización de hiperparámetros

# Hiperparámetros que se quieren comparar (random search)
hiperparametros <- list(
                     learn_rate  = c(0.01, 0.1, 0.01),
                     ntrees      = c(500, 1000, 2000),
                     max_depth   = seq(2, 15, 1),
                     reg_lambda  = c(0, 0.5, 1),
                      reg_alpha  = c(0, 0.5, 1),
                     sample_rate = seq(0.5, 1.0, 0.2)
                    )
# Criterios de parada para la búsqueda
search_criteria <- list(
                    strategy = "RandomDiscrete",
                    max_models = 50, # Número máximo de combinaciones
                    max_runtime_secs = 60*10, # Tiempo máximo de búsqueda
                    stopping_tolerance = 0.001, # Mejora mínima
                    stopping_rounds = 5,
                    seed = 123
                   )

grid_xgboost <- h2o.grid(
    # Algoritmo y parámetros
    algorithm = "xgboost",
    booster   = "gblinear",
    # Variable respuesta y predictores
    y = var_respuesta,
    x = predictores,
    # Datos de entrenamiento.
    training_frame = datos_train_h2o,
    # Preprocesado
    ignore_const_cols = TRUE,
    # Parada temprana
    score_tree_interval = 100,
    stopping_rounds = 5,
    stopping_metric = "rmse",
    stopping_tolerance = 0.01,
    # Hiperparámetros optimizados
    hyper_params = hiperparametros,
    # Estrategia de validación para seleccionar el mejor modelo
    seed   = 123,
    nfolds = 3,
    # Tipo de búsqueda
    search_criteria = search_criteria,
    keep_cross_validation_predictions = TRUE,
    grid_id         = "grid_xgboost"
)
# Se muestran los modelos ordenados de mayor a menor rmse
resultados_grid_xgboost <- h2o.getGrid(
                              grid_id = "grid_xgboost",
                              sort_by = "rmse",
                              decreasing = FALSE
                            )

as.data.frame(resultados_grid_xgboost@summary_table)



Mejor modelo

# Se reentrena el modelo con los mejores hiperparámetros
mejores_hiperparam <- h2o.getModel(resultados_grid_xgboost@model_ids[[1]])@parameters

modelo_xgboost <- h2o.xgboost(
    # Variable respuesta y predictores
    y              = var_respuesta,
    x              = predictores,
    distribution   = "gaussian",
    booster        = "gblinear",
    # Datos de entrenamiento.
    training_frame = datos_train_h2o,
    # Preprocesado
    ignore_const_cols = TRUE,
    # Hiperparámetros.
    learn_rate  = mejores_hiperparam$learn_rate,
    max_depth   = mejores_hiperparam$max_depth,
    ntrees      =  mejores_hiperparam$ntrees,
    sample_rate = mejores_hiperparam$sample_rate,
    reg_lambda  = mejores_hiperparam$reg_lambda,
    reg_alpha   = mejores_hiperparam$reg_alpha,
    # Estrategia de validación para seleccionar el mejor modelo.
    seed            = 123,
    nfolds          = 10,
    keep_cross_validation_predictions = TRUE,
    model_id        = "modelo_xgboost"
)



Influencia predictores


Importancia predictores

h2o.varimp(modelo_xgboost)
h2o.varimp_plot(modelo_xgboost)



Gráficos PDP

pdp_plots <- h2o.partialPlot(
                object = modelo_xgboost,
                data   = datos_train_h2o,
                cols   = predictores,
                nbins  = 20,
                plot   = TRUE,
                plot_stddev = TRUE
              )



Curvas ICE

custom_ice <- function(modelo_h2o, data, y, predictores = NA) {

  predictor <- Predictor$new(
                model = modelo_h2o, 
                data  = data, 
                y     = y, 
                class = "regression"
              )
  
  if(is.na(predictores)) {
    predictores <- colnames(data)
  }
  
  graficos <- list()
  
  for (i in seq_along(predictores)) {
    ice <- FeatureEffect$new(
                  predictor = predictor,
                  feature   =  predictores[i],
                  method    = "pdp+ice",
                  grid.size = 20
           )
    plot_ice <- plot(ice) + ggtitle(predictores[i])
    graficos[[i]] <- plot_ice
  }
  
  return(graficos)
}

graficos_ice <- custom_ice(
                  modelo_h2o = modelo_gbm,
                  data       = datos_train[, predictores],
                  y          = datos_train[[var_respuesta]]
                )

ggarrange(plotlist = graficos_ice, ncol = 2, nrow = 3)



Diagnóstico de los residuos



explainer_xgboost <- DALEXtra::explain_h2o(
                      model = modelo_xgboost,
                      data  = datos_train[, predictores],
                      y     = datos_train[[var_respuesta]],
                      label = "modelo_xgboost" 
                    )
## Preparation of a new explainer is initiated
##   -> model label       :  modelo_xgboost 
##   -> data              :  1571  rows  6  cols 
##   -> target variable   :  1571  values 
##   -> model_info        :  package h2o , ver. 3.28.1.2 , task regression (  default  ) 
##   -> predict function  :  yhat.H2ORegressionModel  will be used (  default  )
##   -> predicted values  :  numerical, min =  196.1371 , mean =  813.9111 , max =  1484.282  
##   -> residual function :  difference between y and yhat (  default  )
##   -> residuals         :  numerical, min =  -934.4101 , mean =  -0.03079092 , max =  1784.308  
##   A new explainer has been created! 
plot(model_performance(explainer_xgboost))

predicciones_train <- h2o.predict(
                        modelo_xgboost,
                        newdata = datos_train_h2o
                      )
predicciones_train <- h2o.cbind(
                       datos_train_h2o["precio"],
                       predicciones_train
                      )
predicciones_train <- as.data.frame(predicciones_train)
predicciones_train <- predicciones_train %>%
                      mutate(
                        residuo = predict - precio
                      )

p1 <- ggplot(predicciones_train, aes(x = precio, y  = predict)) +
      geom_point(alpha = 0.5) +
      geom_smooth(method = "gam", color = "red", size = 1, se = FALSE) +
      labs(title = "Predicciones vs valor real") +
      theme_bw()

p2 <- ggplot(predicciones_train, aes(1:nrow(predicciones_train), y  = residuo)) +
      geom_point(alpha = 0.5) +
      geom_hline(yintercept = 0, color = "red", size = 1) +
      labs(title = "Residuos del modelo") +
      theme_bw()

p3 <- ggplot(predicciones_train, aes(x  = residuo)) +
      geom_density() +
      geom_rug() +
      labs(title = "Residuos del modelo") +
      theme_bw()

p4 <- ggplot(predicciones_train, aes(sample  = predict)) +
      stat_qq() +
      stat_qq_line(color = "red", size = 1) +
      labs(title = "QQ-plot residuos del modelo") +
      theme_bw()

ggpubr::ggarrange(p1, p2, p3, p4, ncol = 2, nrow = 2) %>%
ggpubr::annotate_figure(
  top = ggpubr::text_grob("Diagnóstico residuos entrenamiento",
                          color = "black", face = "bold", size = 14)
)



Predicción test


Predicciones

predicciones_test <- h2o.predict(
                        object  = modelo_xgboost,
                        newdata = datos_test_h2o
                     )
predicciones_test  <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_test

Error test

rmse_test_xgboost <- MLmetrics::RMSE(
                       y_pred = datos_test$precio,
                       y_true = datos_test$prediccion
                     )
paste("Error de test (rmse) del modelo XGBOOST:", rmse_test_xgboost)
## [1] "Error de test (rmse) del modelo XGBOOST: 295.302739166193"



Escritura del modelo


# Se guarda el modelo en el directorio actual en la carpeta modelos
h2o.saveModel(object = modelo_xgboost, path = "./modelos", force = TRUE)
file.rename(file.path("./modelos", modelo_xgboost@model_id), "./modelos/modelo_xgboost.h2o")



ENSEMBLE


Entrenamiento


Para poder entrenar el ensemble, los modelos tienen que ser entrenados en la misma sesión de H2O, no pueden estar cargados de disco.

  • Todos los modelos deben de emplear el mismo número de particiones en la validación cruzada (mismo valor en el argumento nfolds).

  • Deben de almacenarse las predicciones hechas durante la validación cruzada (keep_cross_validation_predictions = True).

  • Todos los modelos deben ser entrenados sobre el mismo conjunto de entrenamiento. Las filas deben de ser las mismas aunque cada modelo puede utilizar distintos predictores (columnas).

modelo_ensemble <- h2o.stackedEnsemble(
                     y = var_respuesta,
                     x = predictores,
                     training_frame = datos_train_h2o,
                     base_models = list(modelo_glm,
                                        modelo_rf,
                                        modelo_xgboost,
                                        modelo_gbm),
                     model_id = "modelo_ensemble"
                    )
modelo_ensemble
## Model Details:
## ==============
## 
## H2ORegressionModel: stackedensemble
## Model ID:  modelo_ensemble 
## Number of Base Models: 4
## 
## Base Models (count by algorithm type):
## 
##     drf     gbm     glm xgboost 
##       1       1       1       1 
## 
## Metalearner:
## 
## Metalearner algorithm: glm
## 
## 
## H2ORegressionMetrics: stackedensemble
## ** Reported on training data. **
## 
## MSE:  74942.17
## RMSE:  273.7557
## MAE:  210.8262
## RMSLE:  0.3612048
## Mean Residual Deviance :  74942.17



Predicción test


Predicciones

predicciones_test <- h2o.predict(
                        object  = modelo_ensemble,
                        newdata = datos_test_h2o
                     )
predicciones_test     <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_test



Error test

rmse_test_ensemble <- MLmetrics::RMSE(
                        y_pred = datos_test$precio,
                        y_true = datos_test$prediccion
                      )
paste("Error de test (rmse) del modelo ensemble:", rmse_test_ensemble)
## [1] "Error de test (rmse) del modelo ensemble: 284.220732184386"



Escritura del modelo


# Se guarda el modelo en el directorio actual en la carpeta modelos
h2o.saveModel(object = modelo_ensemble, path = "./modelos", force = TRUE)
file.rename(file.path("./modelos", modelo_ensemble@model_id), "./modelos/modelo_ensemble.h2o")



Comparación de modelos


Residuos train


p1 <- plot(
         model_performance(explainer_glm),
         model_performance(explainer_gbm),
         model_performance(explainer_rf),
         model_performance(explainer_xgboost)
      )
       
       
p2 <- plot(
         model_performance(explainer_glm),
         model_performance(explainer_gbm),
         model_performance(explainer_rf),
         model_performance(explainer_xgboost),
         geom = "boxplot"
      )
ggpubr::ggarrange(p1, p2, nrow = 2)



Error validación cruzada


extraer_metricas_cv <- function(modelo_h2o, metrica) {
  metricas_cv <- modelo_h2o@model$cross_validation_metrics_summary %>%
                 data.frame() %>%
                 rownames_to_column(var = "metrica") %>%
                 filter(metrica == {{metrica}}) %>%
                 select(starts_with("cv_")) %>%
                 as.numeric()
  return(metricas_cv)
}

metricas_cv_glm <- extraer_metricas_cv(modelo_h2o = modelo_glm, metrica = "rmse")
metricas_cv_gbm <- extraer_metricas_cv(modelo_h2o = modelo_gbm, metrica = "rmse")
metricas_cv_rf <- extraer_metricas_cv(modelo_h2o = modelo_rf, metrica = "rmse")
metricas_cv_xgboost <- extraer_metricas_cv(modelo_h2o = modelo_xgboost, metrica = "rmse")

df_metricas_cv <- data.frame(
                    glm = metricas_cv_glm,
                    gbm = metricas_cv_gbm,
                    rf  = metricas_cv_rf,
                    xgboost = metricas_cv_xgboost
                  ) %>%
                  pivot_longer(
                    dplyr::everything(),
                    names_to = "modelo",
                    values_to = "rmse"
                  )

ggplot(data = df_metricas_cv, aes(x = modelo, y = rmse, color = modelo)) +
  geom_violin() +
  geom_boxplot(outlier.shape = NA, width = 0.2) +
  labs(title = "Comparación errores de validación cruzada") +
  theme_bw() +
  theme(legend.position = "none")



Error test


datos_error_test <- data.frame(
                    rmse_test_glm,
                    rmse_test_gbm,
                    rmse_test_rf,
                    rmse_test_xgboost,
                    rmse_test_ensemble
                  ) %>%
                unlist() %>%
                enframe(name = "modelo", value = "rmse_test") %>%
                arrange(rmse_test)
datos_error_test



Creative Commons Licence
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.