MACHINE LEARNING: H2O, IML y DALEX
Joaquín Amat Rodrigo
Abril, 2020
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 deSpySmindicando 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)| 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
## 0plot_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.0Algunos 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 44Modelos
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.0summary(datos_test$precio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 127.1 537.8 718.6 804.0 1000.0 2764.5summary(datos_train$banyo)/nrow(datos_train)## 0 1
## 0.97963081 0.02036919summary(datos_test$banyo)/nrow(datos_test)## 0 1
## 0.96984925 0.03015075GLM
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)