Versión PDF: Github

Más sobre ciencia de datos: cienciadedatos.net

Machine Learning con H2O y R

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	complete_rate	ordered	n_unique	top_counts
banyo	1	FALSE	2	0: 1925, 1: 44
calefaccion	1	FALSE	2	0: 1580, 1: 389
cocina	1	FALSE	2	0: 1808, 1: 161
situacion	1	FALSE	3	2: 1247, 3: 550, 1: 172

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
precio	1	811.88	379.00	101.7	544.2	737.8	1022	3000	▇▇▂▁▁
metros	1	67.73	20.86	30.0	52.0	67.0	82	120	▆▇▇▅▂
anyo	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 ( [33m default [39m ) 
##   -> predict function  :  yhat.H2ORegressionModel  will be used ( [33m default [39m )
##   -> predicted values  :  numerical, min =  261.1752 , mean =  814.0984 , max =  1798.123  
##   -> residual function :  difference between y and yhat ( [33m default [39m )
##   -> residuals         :  numerical, min =  -789.5428 , mean =  -0.2181171 , max =  1629.115  
##  [32m A new explainer has been created! [39m

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 ( [33m default [39m ) 
##   -> predict function  :  yhat.H2ORegressionModel  will be used ( [33m default [39m )
##   -> predicted values  :  numerical, min =  304.4805 , mean =  813.403 , max =  1722.642  
##   -> residual function :  difference between y and yhat ( [33m default [39m )
##   -> residuals         :  numerical, min =  -797.0642 , mean =  0.477308 , max =  1331.144  
##  [32m A new explainer has been created! [39m

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)

MACHINE LEARNING: H2O, IML y DALEX

Joaquín Amat Rodrigo

Abril, 2020

Introducción

H2O

Dalex y IML

Otros paquetes

Datos

Análisis exploratorio

Valores ausentes

Variables respuesta

Variables continuas

Correlación variables continuas

Variables cualitativas

Modelos

Iniciación del cluster

Train-Test

GLM

Entrenamiento

Influencia predictores

Diagnóstico de los residuos

Predicción test

Escritura del modelo

GBM

Entrenamiento

Influencia predictores

Diagnóstico de los residuos

Predicción test

Escritura del modelo

RF

Entrenamiento

Influencia predictores

Diagnóstico de los residuos

Predicción test

Escritura del modelo

XGBOOST

Entrenamiento

Influencia predictores