- 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
- Diagnóstico de los residuos
- Predicción test
- Escritura del modelo
- ENSEMBLE
- Entrenamiento
- Predicción test
- Escritura del modelo
- Comparación de modelos
- Residuos train
- Error validación cruzada
- Error test
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)
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_testError 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! [39mplot(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_testError 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! [39mplot(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_testError 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 ( [33m default [39m )
## -> predict function : yhat.H2ORegressionModel will be used ( [33m default [39m )
## -> predicted values : numerical, min = 196.1371 , mean = 813.9111 , max = 1484.282
## -> residual function : difference between y and yhat ( [33m default [39m )
## -> residuals : numerical, min = -934.4101 , mean = -0.03079092 , max = 1784.308
## [32m A new explainer has been created! [39mplot(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_testError 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.17Predicción test
Predicciones
predicciones_test <- h2o.predict(
object = modelo_ensemble,
newdata = datos_test_h2o
)
predicciones_test <- as.vector(predicciones_test)
datos_test$prediccion <- predicciones_testError 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
Este material, creado por Joaquín Amat Rodrigo, tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.
Se permite:
-
Compartir: copiar y redistribuir el material en cualquier medio o formato.
-
Adaptar: remezclar, transformar y crear a partir del material.
Bajo los siguientes términos:
-
Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.
-
NoComercial: No puedes utilizar el material para fines comerciales.
-
CompartirIgual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.