Machine Learning con H2O y R
Joaquín Amat Rodrigo
Julio, 2018 (última actualización Abril 2020)
Más sobre ciencia de datos en cienciadedatos.net
Introducción
En este documento se pretende mostrar cómo crear modelos de machine learning combinando H2O y el lenguaje de programación R.
Aunque son muchos los pasos que preceden al entrenamiento de un modelo (exploración de los datos, transformaciones, selección de predictores, etc.), para no añadir una capa extra de complejidad, se va a asumir que los datos se encuentran prácticamente preparados para que puedan ser aceptados por los algoritmos de aprendizaje. Se puede encontrar una descripción más detallada de las etapas que forman parte el ciclo de vida del modelado en los documentos mlr3 y tidymodels.
¿Qué es H2O?
H2O es un producto creado por la compañía H2O.ai con el objetivo de combinar los principales algoritmos de machine learning y aprendizaje estadístico con el Big Data. Gracias a su forma de comprimir y almacenar los datos, H2O es capaz de trabajar con millones 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.
Aunque la principal ventaja de H2O frente a otras herramientas es su escalabilidad, sus algoritmos son igualmente útiles cuando se trabaja con un volumen de datos reducido.
Comunicación entre H2O y R
El manejo de H2O puede hacerse íntegramente desde R: iniciar el cluster, carga de datos, entrenamiento de modelos, predicción de nuevas observaciones, etc. Es importante tener en cuenta que, aunque los comandos se ejecuten desde R, los datos se encuentran en el cluster de H2O, no en memoria. Solo cuando los datos se cargan en memoria, se les pueden aplicar funciones propias de R.
Las funciones as.data.frame() y as.h2o() permiten transferir los datos de la sesión de R al cluster H2O y viceversa. Hay que tener especial precaución cuando el movimiento de datos se hace desde H2O a R, ya que esto implica cargar en RAM todos los datos y, si son muchos, pueden ocupar toda la memoria. Para evitar este tipo de problemas, conviene realizar todos los pasos posibles (filtrado, agregaciones, cálculo de nuevas columnas…) con las funciones de H2O antes de transferir los datos.
Datos
Para los ejemplos de este documento se emplea el set de datos Adult Data Set disponible en UCI Machine Learning Repository. Este set de datos contiene información sobre 48842 ciudadanos de diferentes países. Para cada ciudadano, se dispone de 15 variables:
age: edad.workclass: tipo de empleo.fnlwgt: importancia relativa de cada observación.education: nivel de estudios.education-num: años de estudios.marital-status: estado civil.occupation: sector de trabajorelationship: estado familiar.race: raza.sex: género.capital-gain: ganancia de capital.capital-loss: perdidas de capital.hours-per-week: horas de trabajo semanal.native-countr: país de origen.salario: si el salario supera o no 50000 dólares anuales.
El objetivo del problema es crear un modelo capaz de predecir correctamente si un ciudadano tiene un salario superior o inferior a 50000 dólares anuales. Se trata de un problema de clasificación binaria supervisado.
Los datos se han modificado para eliminar observaciones con valores ausentes y para recodificar los niveles de algunas variables. El archivo resultante (adult_custom.csv) puede descargarse directamente del repositorio Github.
# EL SIGUIENTE CÓDIGO MUESTRA COMO SE HAN PREPROCESADO LOS DATOS ORIGINALES.
# NO ES NECESARIO EJECUTARLO SI SE HA DESCARGADO EL ARCHIVO adult_custom.csv.
library(tidyverse)
# Datos de train.
datos_train <- read_csv(file = "adult_train.csv", col_names = FALSE)
colnames(datos_train) <- c("age", "workclass", "final_weight", "education",
"education_number", "marital_status", "occupation",
"relationship", "race", "sex", "capital_gain",
"capital_loss", "hours_per_week", "native_country",
"salario")
# Datos de test.
datos_test <- read_csv(file = "adult_test.csv", col_names = FALSE, skip = 1)
colnames(datos_test) <- c("age", "workclass", "final_weight", "education",
"education_number", "marital_status", "occupation",
"relationship", "race", "sex", "capital_gain",
"capital_loss", "hours_per_week", "native_country",
"salario")
# Se combinan los dos conjuntos de datos.
datos <- bind_rows(datos_train, datos_test)
# Se eliminan registros con missing values, que en este set de datos están
# codificados como ?.
remplazar <- function(x){
x[x == "?"] <- NA
return(x)
}
datos <- datos %>% map_df(.f = remplazar)
datos <- drop_na(datos)
# Recodificación de niveles.
datos <- datos %>% mutate(workclass = recode(datos$workclass,
"Without_pay" = "desempleado",
"Self_emp_inc" = "autonomo",
"Self_emp_not_inc" = "autonomo",
"Federal_gov" = "funcionario",
"Local_gov" = "funcionario",
"State_gov" = "funcionario"))
datos <- datos %>% mutate(marital_status = recode(datos$marital_status,
"Married-AF-spouse" = "casado",
"Married-civ-spouse" = "casado",
"Married-spouse-absent"= "casado",
"Widowed" = "separado",
"Never-married" = "soltero",
"Separated" = "separado",
"Divorced" = "separado")
)
norte_america <- c("Canada", "Cuba", "Dominican-Republic", "El-Salvador",
"Guatemala", "Haiti", "Honduras", "Jamaica", "Mexico",
"Nicaragua", "Outlying-US(Guam-USVI-etc)", "Puerto-Rico",
"Trinadad&Tobago", "United-States")
europa <- c("England", "France", "Germany", "Greece", "Holand-Netherlands",
"Hungary", "Ireland", "Italy", "Poland", "Portugal", "Scotland",
"Yugoslavia")
asia <- c("Cambodia", "China", "Hong", "India", "Iran", "Japan", "Laos",
"Philippines", "Taiwan", "Thailand", "Vietnam")
sud_america <- c("Columbia", "Ecuador", "Peru")
otros <- c("South")
datos <- datos %>% mutate(native_country = case_when(
native_country %in% norte_america ~ "Norte America",
native_country %in% asia ~ "Asia",
native_country %in% sud_america ~ "Sud America",
native_country %in% europa ~ "Europa",
native_country %in% otros ~ "Otros",
TRUE ~ as.character(native_country)
))
primaria <- c("Preschool", "1st-4th", "5th-6th", "7th-8th", "9th", "10th",
"11th","12th")
bachillerato <- c("Some-college", "Assoc-acdm", "Assoc-voc")
master <- c("Masters", "Prof-school")
datos <- datos %>% mutate(education = case_when(
education %in% primaria ~ "primaria",
education == "HS-grad" ~ "secundaria",
education %in% bachillerato ~ "bachillerato",
education == "Bachelors" ~ "universidad",
education %in% master ~ "master",
education == "Doctorate" ~ "doctorado",
TRUE ~ as.character(education)
))
datos <- datos %>% mutate(salario = case_when(
salario == "<=50K." ~ "<=50K",
salario == ">50K." ~ ">50K",
TRUE ~ as.character(salario)
))
# Se agrupan las observaciones cuya variable occupation == "Armed-Forces" con
# "Other-service" ya que solo hay 14 observaciones (0.03% del total).
datos <- datos %>% mutate(occupation = case_when(
occupation == "Armed-Forces" ~ "Other-service",
TRUE ~ as.character(occupation)
))
write_csv(x = datos, path = "adult_custom.csv", col_names = TRUE)Iniciación de H2O
Una vez que H2O ha sido instalado (install.packages("h2o")), hay que iniciarlo, bien en todo el cluster o en un solo ordenador. Para este ejemplo, se emplea un único ordenador del que se utilizan todos sus cores en paralelo.
# 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"
)##
## H2O is not running yet, starting it now...
##
## Note: In case of errors look at the following log files:
## /tmp/Rtmp45iDsZ/file37ce5a66a289/h2o_ubuntu_started_from_r.out
## /tmp/Rtmp45iDsZ/file37ce11f8ef69/h2o_ubuntu_started_from_r.err
##
##
## Starting H2O JVM and connecting: .. Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 1 seconds 551 milliseconds
## H2O cluster timezone: Etc/UTC
## H2O data parsing timezone: UTC
## H2O cluster version: 3.32.0.1
## H2O cluster version age: 1 month and 12 days
## H2O cluster name: H2O_started_from_R_ubuntu_rek408
## H2O cluster total nodes: 1
## H2O cluster total memory: 5.33 GB
## H2O cluster total cores: 8
## H2O cluster allowed cores: 8
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## H2O Internal Security: FALSE
## H2O API Extensions: Amazon S3, XGBoost, Algos, AutoML, Core V3, TargetEncoder, Core V4
## R Version: R version 4.0.3 (2020-10-10)# Se eliminan los datos del cluster por si ya había sido iniciado.
h2o.removeAll()
# Para que no se muestre la barra de progreso.
h2o.no_progress()Tras iniciar el cluster (local), se muestran por pantalla sus características, entre las que están: el número de cores activados (4), la memoria total del cluster (5.32 GB), el número de nodos (1 porque se está empleando un único ordenador) y el puerto con el que conectarse a la interfaz web de H2O (http://localhost:54321/flow/index.html).
Si se desea lanzar H2O en un cluster Hadoop ya establecido, solo es necesario especificar la dirección IP y puerto de acceso en h2O.init().
Carga de datos
La carga de datos puede hacerse directamente al cluster H2O, o bien cargándolos primero en memoria en la sesión de R y después transfiriéndolos. La segunda opción no es aconsejable si el volumen de datos es muy grande.
# Carga de datos en el cluster H2O desde url.
url_path <- paste0("https://github.com/JoaquinAmatRodrigo/Estadistica-con-R/raw/",
"master/datos/adult_custom.csv")
datos_h2o <- h2o.importFile(
path = url_path,
header = TRUE,
sep = ",",
destination_frame = "datos_h2o"
)# Carga de datos en el cluster H2O desde local.
datos_h2o <- h2o.importFile(
path = "./adult_custom.csv",
header = TRUE,
sep = ",",
destination_frame = "datos_h2o"
)# Carga de datos en R y transferencia a H2O.
url_path <- paste0("https://github.com/JoaquinAmatRodrigo/Estadistica-con-R/raw/",
"master/datos/adult_custom.csv")
datos_R <- read.csv(file = url_path, header = TRUE)
datos_h2o <- as.h2o(x = datos_R, destination_frame = "datos_h2o")Exploración de los datos
Si bien el set de datos empleado en este ejemplo es suficientemente pequeño para cargarlo todo en memoria y emplear directamente las magníficas posibilidades de análisis exploratorio que ofrece R, para representar mejor la forma de proceder con grandes volúmenes de datos, se emplean funciones propias de H2O.
## [1] 45222 15## [1] "age" "workclass" "final_weight" "education"
## [5] "education_number" "marital_status" "occupation" "relationship"
## [9] "race" "sex" "capital_gain" "capital_loss"
## [13] "hours_per_week" "native_country" "salario"La función h2o.describe() es muy útil para obtener un análisis rápido que muestre el tipo de datos, la cantidad de valores ausentes, el valor mínimo, máximo, media, desviación típica y el número de categorías (Cardinality) de cada una de las variables. H2O emplea el nombre enum para los datos de tipo factor o character.
Para conocer el índice o nombre de las columnas que son de un determinado tipo, por ejemplo, numérico, se emplea la función h2o.columns_by_type().
## [1] 1 3 5 11 12 13## [1] "age" "final_weight" "education_number" "capital_gain"
## [5] "capital_loss" "hours_per_week"Con la función h2o.cor() se calcula la correlación entre dos o más columnas numéricas.
indices <- h2o.columns_by_type(object = datos_h2o, coltype = "numeric")
h2o.cor(x = datos_h2o[, indices], y = NULL, method = "Pearson", na.rm = TRUE)Para contar el número de observaciones de cada clase en una variable categórica, como es en este caso la variable respuesta salario, se emplea la función h2o.table().
# Se crea una tabla con el número de observaciones de cada tipo.
tabla_muestra <- as.data.frame(h2o.table(datos_h2o$salario))
tabla_muestra# Una vez creada la tabla, se carga en el entorno de R para poder graficar.
ggplot(
data = tabla_muestra,
aes(x = salario, y = Count, fill = salario)) +
geom_col() +
scale_fill_manual(values = c("gray50", "orangered2")) +
theme_bw() +
labs(
x = "Salario", y = "Número de observaciones",
title = "Distribución de la variable Salario") +
theme(legend.position = "none")
Separación de training, validación y test
Evaluar la capacidad predictiva de un modelo consiste en comprobar cómo de próximas son sus predicciones a los verdaderos valores de la variable respuesta. Para poder cuantificar de forma correcta este error, se necesita disponer de un conjunto de observaciones, de las que se conozca la variable respuesta, pero que el modelo no haya “visto”, es decir, que no hayan participado en su ajuste. Con esta finalidad, se separan los datos disponibles en un conjunto de entrenamiento, uno de validación y uno de test. El conjunto de entrenamiento se emplea para ajustar el modelo, el de validación para encontrar los mejores hiperparámetros (model tuning) y el de test para estimar el error que comete el modelo al predecir nuevos datos. El tamaño adecuado de las particiones depende en gran medida de la cantidad de datos disponibles y la seguridad que se necesite en la estimación del error, 70%-15%-15% suele dar buenos resultados.
Si el número de observaciones es limitado, crear 3 particiones puede generar grupos demasiado pequeños y variables. En estos casos, es preferible dividir los datos únicamente en un conjunto de entrenamiento y otro de test, y utilizar validación cruzada (cross validation) sobre el conjunto de entrenamiento durante la optimización de los hiperparámetros del modelo. Para mostrar ambas posibilidades, en el ejemplo de modelo GLM se empleará validación cruzada, mientras que en el resto de modelos (GBM, RF y Deep Learning), se emplea un solo conjunto de validación. Es importante tener en cuenta que, la validación cruzada, requiere ajustar el modelo repetidas veces, lo que supone un coste computacional muy alto. Cuando se trabaja con big data no suele ser factible aplicar validación cruzada, además, al disponer de tantos datos, un conjunto de validación único es suficiente para obtener buenas estimaciones del comportamiento del modelo. Para una descripción más detallada de los diferentes métodos de validación, consultar el Anexo4 del documento Machine learning con R y Caret.
La función h2o.splitFrame() realiza particiones aleatorias, pero no permite hacerlas de forma estratificada, por lo que no asegura que la distribución de clases de variable respuesta sea igual en todas particiones. Esto puede ser problemático con datos muy desbalanceados (alguno de los grupos es muy minoritario).
# Separación de las observaciones en conjunto de entrenamiento y test.
# En los ejemplos de GBM y deep learning se repetirá la separación, pero en
# tres conjuntos en lugar de dos.
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")Número de observaciones de cada clase por partición.
## salario Count
## 1 <=50K 27229
## 2 >50K 9009
##
## [2 rows x 2 columns]## salario Count
## 1 <=50K 6785
## 2 >50K 2199
##
## [2 rows x 2 columns]## salario Count
## 1 NaN 0.7513936
## 2 NaN 0.2486064
##
## [2 rows x 2 columns]## salario Count
## 1 NaN 0.7552315
## 2 NaN 0.2447685
##
## [2 rows x 2 columns]La proporción de los dos grupos de la variable respuesta es muy similar en ambas particiones.
Preprocesado de datos
H2O incorpora y automatiza muchas de las transformaciones necesarias para que los datos puedan ser ingeridos por los algoritmos de machine learning. Esto es distinto a la mayoría de librerías de machine learning de otros lenguajes como Python y R, donde las etapas de preprocesado suelen ser independientes del entrenamiento del modelo. En concreto, H2O automatiza las siguientes transformaciones:
Variables categóricas en H2O
En el momento de ajustar un modelo (GLM, GBM, DRF, Deep Learning, K-Means, Aggregator, XGBoost), H2O identifica automáticamente que variables son categóricas y crea internamente las variables dummy correspondientes. Es altamente recomendable permitir que H2O realice este proceso en lugar de hacerlo externamente, ya que su implementación está muy optimizada. El comportamiento de la codificación de las variables categóricas puede controlarse con el argumento categorical_encoding, por defecto su valor es "AUTO".
Estandarización
Por defecto, H2O estandariza los predictores numéricos antes de ajustar los modelos para que todos tengan media cero y varianza uno. Este comportamiento se puede controlar con el argumento standardize. Es importante tener en cuenta que, para muchos modelos (lasso, ridge, Deep Learning…), es necesario realizar la estandarización de predictores.
Eliminación de variables con varianza cero
No se deben de incluir en un modelo predictores que contengan un único valor (varianza cero), ya que no aportan información. Los algoritmos de H2O excluyen directamente las columnas con valor constante. Este comportamiento se puede controlar mediante el argumento ignore_const_cols.
Balance de clases
En problemas de clasificación, puede ocurrir que los datos estén desbalanceados, es decir, que el número de observaciones pertenecientes a cada grupo sea muy dispar. Por ejemplo, que de 1000 observaciones, 800 pertenezcan a un grupo y solo 200 al otro. En estos casos, los modelos pueden tener dificultades para aprender a identificar las observaciones del grupo minoritario. Con el argumento balance_classes se puede indicar que antes de ajustar el modelo se equilibren las clases mediante undersampling u oversampling. Por defecto, este comportamiento está desactivado.
Todas las trasformaciones descritas anteriormente se aprenden a partir de los datos de entrenamiento en el momento del ajuste, y se aplican automáticamente cuando el modelo se emplea para predecir nuevas observaciones. De esta forma, se garantiza que no se viola la condición de que ninguna información procedente de las observaciones de test participe o influya en el ajuste del modelo.
H2O dispone de un listado de funciones para manipular y transformar los datos, sin embargo, con frecuencia no son suficientes para resolver todas las necesidades de preprocesado. En estos casos es necesario recurrir a R, si el volumen de datos es limitado, o a Spark si la información no se puede cargar en memoria.
Nota: en mi experiencia propia, aunque las tecnologías distribuidas como Spark son muy útiles, antes de recurrir a ellas, recomiendo reflexionar sobre la manipulación que se tiene que realizar y analizar si es necesario que se haga sobre todos los datos a la vez. Por ejemplo, el proceso de imputación se debe aprender únicamente de los datos de entrenamiento y luego aplicarlo al resto de datos. Puede que todo el dataset no se pueda cargar en memeoria, pero sí el conjunto de entrenamiento. Esto permitiría aprender la imputación y después aplicarla por separado a cada conjunto.
Modelos con H2O
Modelos disponibles
H2O incorpora los siguientes algoritmos de machine learning. Todos ellos están implementados de forma que puedan trabajar, en la medida de lo posible, de forma distribuida y/o en paralelo.
Modelos supervisados
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
Modelos no supervisados
Aggregator
Generalized Low Rank Models (GLRM)
K-Means Clustering
Isolation Forest
Principal Component Analysis (PCA)
Otros
- Word2vec
En este documento se muestran ejemplos con Generalized Linear Model (GLM), Gradient Boosting Machine (GBM), Deep Learning (Neural Networks) y Stacked Ensembles.
Parada temprana
Los algoritmos de H2O están pensados para poder ser utilizados de forma eficiente con grandes volúmenes de datos. Emplear miles o millones de datos suele implicar tiempos largos de entrenamiento. H2O incorpora criterios de parada para que no se continúe mejorando un modelo si ya se ha alcanzado una solución aceptable. Dependiendo del algoritmo, el criterio de parada es distinto. En los siguientes apartados se describen con detalle.
Optimización de hiperparámetros
Muchos modelos contienen parámetros que no pueden aprenderse a partir de los datos de entrenamiento y que, por lo tanto, deben de ser establecidos por el analista. A estos se les conoce como hiperparámetros. Los resultados de un modelo pueden depender en gran medida del valor que tomen sus hiperparámetros, sin embargo, no se puede conocer de antemano cuál es el adecuado. Aunque con la práctica, los especialistas en machine learning ganan intuición sobre qué valores pueden funcionar mejor en cada problema, no hay reglas fijas. La forma más común de encontrar los valores óptimos es probando diferentes posibilidades, lo que comúnmente se conoce como tunning.
H2O emplea la función h2o.grid() para realizar la búsqueda de los mejores hiperparámetros, sus argumentos principales son: el nombre del algoritmo, los parámetros del algoritmo, una lista con los valores de los hiperparámetros que se quieren comparar, el tipo de búsqueda (“Cartesian” o “RandomDiscrete”) y, si es de tipo random, un criterio de parada. Una vez que la búsqueda ha finalizado, el objeto grid creado contiene todos los modelos, para acceder a ellos es necesario extraerlos.
Generalized linear models (GLMs)
Introducción
El modelo de regresión lineal (Legendre, Gauss, Galton y Pearson) considera que, dado un conjunto de observaciones, la media \(\mu\) de la variable respuesta \(Y\) se relaciona de forma lineal con la o las variables regresoras \(X\).
\[\mu = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}\]
o en notación matricial (incorporando \(\beta_0\) en el vector \(\mathbf{\beta}\):
\[\mu = \mathbf{X}^T \mathbf{\beta}\]
Ajustar el modelo consiste en estimar los valores de los coeficientes de regresión \(\hat{\beta}\) y la varianza \(\hat{\sigma}^2\) que maximizan la verosimilitud (likelihood) de los datos, es decir, los que dan lugar al modelo que con mayor probabilidad puede haber generado los datos observados. El método empleado con más frecuencia es el ajuste por mínimos cuadrados (OLS):
\[\hat{\beta} = \underset{\beta_0,\beta}{\arg\min} (Y - \mathbf{X}^T \mathbf{\beta})^2\]
\[\hat{\mu} = \mathbf{X}^T \mathbf{\hat{\beta}}\]
\[\hat{\sigma}^2 = \frac{\sum^n_{i=1} \hat{\epsilon}_i^2}{n-p} = \frac{\sum^n_{i=1} (y_i - \hat{\mu})^2}{n-p}\]
Para que esta aproximación sea válida, se necesita que el error se distribuya de forma normal y que la varianza sea constante. El requisito de que la variable respuesta se distribuya de forma normal hace que el modelo de regresión lineal no sea aplicable a muchos problemas reales. Los modelos GLM (John Nelder y Robert Wedderburn) permiten superar esta limitación haciendo una abstracción a un modelo más genérico, en el que la variable respuesta puede seguir cualquier distribución de la familia exponencial (Gaussian, Poisson, binomial…). Para conseguirlo, la media de la variable respuesta ya no está relacionada directamente con los predictores, sino que se hace a través de una función link \(g(.)\). Por ejemplo, el modelo lineal generalizado es equivale al modelo de regresión lineal por mínimos cuadrados cuando la función link es la identidad, o a la regresión logística cuando es binomial (logit).
\[g(\mu) = \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + ... + \beta_p x_{p}\]
\[g(\mu) = \mathbf{X}^T \mathbf{\beta}\]
Además, el modelo lineal generalizado puede incluir regularización durante su ajuste, por lo que también incluye los modelos ridge regression, lasso y elastic net.
Cuando se emplea una función link gausiana (con o sin regularización L2), la máxima verosimilitud equivale a minimizar la suma de errores cuadrados, lo que tiene una solución analítica. Para el resto de familias o cuando se emplea la regularización L1, no existe una solución analítica con la que encontrar la máxima verosimilitud, por lo que se requiere un método de optimización iterativo como el iteratively reweighted least squares (IRLSM), L-BFGS, método de Newton o descenso de gradiente. H2O selecciona automáticamente el método apropiado en dependiendo de la función link especificada por el usuario.
Regularización
H2O incorpora 3 tipos de regularización para los modelos GLM con el objetivo de evitar overfitting, reducir varianza y atenuar el efecto de la correlación entre predictores. Por lo general, aplicando regularización se consigue modelos con mayor poder predictivo (generalización).
El modelo lasso es un modelo lineal por mínimos cuadrados que incorpora una regularización que penaliza la suma del valor absolutos de los coeficientes de regresión \((||\beta||_1 = \sum_{k=1} ^p |\beta_k|)\). A esta penalización se le conoce como l1 y tiene el efecto de forzar a que los coeficientes de los predictores tiendan a cero. Dado que un predictor con coeficiente de regresión cero no influye en el modelo, lasso consigue seleccionar los predictores más influyentes. El grado de penalización está controlado por el hiperparámetro \(\lambda\). Cuando \(\lambda = 0\), el resultado es equivalente al de un modelo lineal por mínimos cuadrados ordinarios. A medida que \(\lambda\) aumenta, mayor es la penalización y más predictores quedan excluidos.
El modelo ridge es un modelo lineal por mínimos cuadrados que incorpora una regularización que penaliza la suma de los coeficientes elevados al cuadrado \((||\beta||^2_2 = \sum_{k=1} ^p \beta^2_k)\). A esta penalización se le conoce como l2 y tiene el efecto de reducir de forma proporcional el valor de todos los coeficientes del modelo pero sin que estos lleguen a cero. Al igual que lasso, el grado de penalización está controlado por el hiperparámetro \(\lambda\).
La principal diferencia práctica entre lasso y ridge es que el primero consigue que algunos coeficientes sean exactamente cero, por lo que realiza selección de predictores, mientras que el segundo no llega a excluir ninguno. Esto supone una ventaja notable de lasso en escenarios donde no todos los predictores son importantes para el modelo y se desea que los menos influyentes queden excluidos. Por otro lado, cuando existen predictores altamente correlacionados (linealmente), ridge reduce la influencia de todos ellos a la vez y de forma proporcional, mientras que lasso tiende a seleccionar uno de ellos, dándole todo el peso y excluyendo al resto. En presencia de correlaciones, esta selección varía mucho con pequeñas perturbaciones (cambios en los datos de entrenamiento), por lo que, las soluciones de lasso, son muy inestables si los predictores están altamente correlacionados.
Para conseguir un equilibrio óptimo entre estas dos propiedades, se puede emplear lo que se conoce como penalización elastic net, que combina ambas estrategias.
- El modelo elastic net incluye una regularización que combina la penalización l1 y l2 \((\alpha \lambda ||\beta||_1 + \frac{1}{2}(1- \alpha)||\beta||^2_2)\). El grado en que influye cada una de las penalizaciones está controlado por el hiperparámetro \(\alpha\). Su valor debe estar comprendido en el intervalo [0,1], cuando \(\alpha = 0\), H2O aplica ridge regression y cuando \(\alpha = 1\) aplica lasso. La combinación de ambas penalizaciones suele dar lugar a buenos resultados. Una estrategia frecuentemente utilizada es asignarle casi todo el peso a la penalización l1 (\(\alpha\) muy próximo a 1) para conseguir seleccionar predictores y un poco a la l2 para dar cierta estabilidad en el caso de que algunos predictores estén correlacionados.
Encontrar el mejor modelo GLM implica identificar los valores óptimos de los hiperparámetros de regularización \(\alpha\) y \(\lambda\). Por defecto, H2O selecciona \(\alpha = 0.5\) y, con cada uno de ellos, analiza un rango de valores \(\lambda\). También es posible hacer una búsqueda del valor óptimo de \(\alpha\) mediante un grid seach.
Lambda Search
La búsqueda del valor óptimo de lambda se habilita indicando lambda_search=TRUE. Además, esta búsqueda puede ser configurada con los siguientes argumentos:
alpha: el valor de \(\alpha\) que se emplea en el modelo para distribuir la penalización l1 y l2. Por defecto,alpha=0.5.validation_frame: método de validación empleado para identificar el mejor modelo, puede ser mediante un único conjunto de validación o mediante validación cruzada (cross-validation). En el segundo caso, no es necesario crear un conjunto de validación, las particiones se generan a partir del conjunto de entrenamiento.n_folds: número de particiones en la validación cruzada.lambda_min_ratioynlambdas: determinan la secuencia de valores lambda que se generan automáticamente y sobre los que se realiza la búsqueda. El rango de valores va desde el valor mínimo con el que se alcanza una penalización total (lambda_max), es decir, el valor de lambda a partir del cual el coeficiente de todos los predictores es cero (valor identificado automáticamente por H2O), hasta un valor de lambda igual a lambda_min_ratio x lambda_max. Dentro de ese rango, se generan un número de lambdas igual al valor especificado connlambdas. Por defecto, se empleannlambda=100ylambda_min_ratio=1e−4.max_active_predictors: establece el número máximo de predictores incluidos en el modelo (predictores con coeficientes de regresión distintos de cero). Esta opción es muy útil para evitar tiempos de computación excesivos cuando se trabaja con sets de datos que contienen cientos de predictores pero se sabe que solo unos pocos son importantes.
Método de ajuste
Todos los modelos GLM se ajustan encontrando el valor de los coeficientes \(\beta\) que maximizan el likelihood. La estimación de estos coeficientes se realiza mediante métodos iterativos, ya que, a excepción de unos pocos casos, no existe una solución analítica. H2O incorpora dos métodos (solvers) distintos, cuya correcta elección, que depende del tipo de datos, puede influir notablemente en el rendimiento computacional.
Iteratively Reweighted Least Squares Method (IRLSM): funciona especialmente bien cuando el set de datos empleado tiene muchas observaciones (largo) pero un número limitado de predictores (estrecho). Los autores recomiendan utilizarlo cuando el número de predictores no es mayor de 500.
Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS): es el método recomendado cuando se dispone de muchos predictores (>500).
Las anteriores son solo recomendaciones generales, en caso de duda, la mejor forma de elegir es comparando el resultado de ambos solvers.
Criterio de parada
Los algoritmos de H2O están pensados para poder ser aplicados a grandes volúmenes de datos, por lo que, entre sus características, incorporan criterios de parada para que no se continúe mejorando un modelo si ya se ha alcanzado una solución aceptable. En el caso de los GLM, el criterio de parada es el número máximo de predictores permitidos, que se especifica mediante el argumento max_active_predictors.
Interacciones
Una de las limitaciones de los modelos GLM es que asumen que la relación entre cada uno de los predictores y la variable respuesta es independiente del resto de predictores, es decir, no contemplan automáticamente interacciones. Una forma de incorporar esta información es creando nuevas variables multiplicando pares de predictores. Este comportamiento puede conseguirse con el argumento interactions e interaction_pairs.
##Modelo
El problema que se describió en la introducción sobre predecir el salario es un problema de clasificación binaria. El principal modelo lineal empleado para clasificaciones binarias es la regresión logística, que es el resultado de emplear, en un modelo lineal generalizado, la función logit como link. Aunque la selección se puede hacer manualmente, H2O es capaz de identificar el tipo de variable respuesta (numérica o categórica) y en función de eso seleccionar automáticamente el tipo de modelo, regresión o clasificación.
# Se comprueba que la variable respuesta es de tipo factor.
datos_train_h2o$salario <- h2o.asfactor(datos_train_h2o$salario)
datos_test_h2o$salario <- h2o.asfactor(datos_test_h2o$salario)
h2o.isfactor(datos_train_h2o$salario)## [1] TRUE# Se define la variable respuesta y los predictores.
var_respuesta <- "salario"
# Para este modelo se emplean todos los predictores disponibles.
predictores <- setdiff(h2o.colnames(datos_h2o), var_respuesta)# Ajuste del modelo y validación mediente 5-CV para estimar su error.
modelo_glm <- h2o.glm(
y = var_respuesta,
x = predictores,
training_frame = datos_train_h2o,
family = "binomial",
link = "logit",
standardize = TRUE,
balance_classes = FALSE,
ignore_const_cols = TRUE,
# Se especifica que hacer con observaciones incompletas
missing_values_handling = "Skip",
# Se hace una búsqueda del hiperparámetro lamba
lambda_search = TRUE,
# Selección automática del solver adecuado
solver = "AUTO",
alpha = 0.95,
# Validación cruzada de 5 folds para estimar el error
# del modelo.
seed = 123,
nfolds = 5,
# Reparto estratificado de las observaciones en la creación
# de las particiones.
fold_assignment = "Stratified",
keep_cross_validation_predictions = FALSE,
model_id = "modelo_glm"
)## Model Details:
## ==============
##
## H2OBinomialModel: glm
## Model ID: modelo_glm
## GLM Model: summary
## family link regularization
## 1 binomial logit Elastic Net (alpha = 0.95, lambda = 5.09E-5 )
## lambda_search
## 1 nlambda = 100, lambda.max = 0.1519, lambda.min = 5.09E-5, lambda.1se = 5.21E-4
## number_of_predictors_total number_of_active_predictors number_of_iterations
## 1 53 45 136
## training_frame
## 1 RTMP_sid_aeaa_11
##
## Coefficients: glm coefficients
## names coefficients standardized_coefficients
## 1 Intercept -7.429849 -2.660849
## 2 occupation.Adm-clerical 0.000000 0.000000
## 3 occupation.Craft-repair 0.017865 0.017865
## 4 occupation.Exec-managerial 0.695254 0.695254
## 5 occupation.Farming-fishing -1.067297 -1.067297
##
## ---
## names coefficients standardized_coefficients
## 49 age 0.025490 0.337455
## 50 final_weight 0.000001 0.073707
## 51 education_number 0.203311 0.519740
## 52 capital_gain 0.000313 2.373427
## 53 capital_loss 0.000630 0.255137
## 54 hours_per_week 0.029095 0.348779
##
## H2OBinomialMetrics: glm
## ** Reported on training data. **
##
## MSE: 0.1046908
## RMSE: 0.3235596
## LogLoss: 0.3269983
## Mean Per-Class Error: 0.1890437
## AUC: 0.9045265
## AUCPR: 0.7689901
## Gini: 0.8090531
## R^2: 0.4395605
## Residual Deviance: 23699.53
## AIC: 23791.53
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 23127 4102 0.150648 =4102/27229
## >50K 2049 6960 0.227439 =2049/9009
## Totals 25176 11062 0.169739 =6151/36238
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.322194 0.693538 226
## 2 max f2 0.145909 0.792705 304
## 3 max f0point5 0.585860 0.709997 128
## 4 max accuracy 0.503406 0.848474 158
## 5 max precision 0.996336 0.992322 2
## 6 max recall 0.000953 1.000000 399
## 7 max specificity 0.999858 0.999890 0
## 8 max absolute_mcc 0.401537 0.587917 195
## 9 max min_per_class_accuracy 0.278949 0.816961 246
## 10 max mean_per_class_accuracy 0.244960 0.821457 261
## 11 max tns 0.999858 27226.000000 0
## 12 max fns 0.999858 8750.000000 0
## 13 max fps 0.000953 27229.000000 399
## 14 max tps 0.000953 9009.000000 399
## 15 max tnr 0.999858 0.999890 0
## 16 max fnr 0.999858 0.971251 0
## 17 max fpr 0.000953 1.000000 399
## 18 max tpr 0.000953 1.000000 399
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
##
## H2OBinomialMetrics: glm
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
##
## MSE: 0.1051063
## RMSE: 0.324201
## LogLoss: 0.3282526
## Mean Per-Class Error: 0.1914324
## AUC: 0.9036885
## AUCPR: 0.7671405
## Gini: 0.8073771
## R^2: 0.4373362
## Residual Deviance: 23790.44
## AIC: 23880.44
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 23278 3951 0.145103 =3951/27229
## >50K 2142 6867 0.237762 =2142/9009
## Totals 25420 10818 0.168138 =6093/36238
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.331430 0.692692 224
## 2 max f2 0.143851 0.791580 306
## 3 max f0point5 0.578472 0.709330 135
## 4 max accuracy 0.480840 0.847922 169
## 5 max precision 0.988455 0.989614 6
## 6 max recall 0.000886 1.000000 399
## 7 max specificity 0.999845 0.999853 0
## 8 max absolute_mcc 0.383540 0.585427 203
## 9 max min_per_class_accuracy 0.278878 0.816628 245
## 10 max mean_per_class_accuracy 0.239945 0.820311 263
## 11 max tns 0.999845 27225.000000 0
## 12 max fns 0.999845 8750.000000 0
## 13 max fps 0.000886 27229.000000 399
## 14 max tps 0.000886 9009.000000 399
## 15 max tnr 0.999845 0.999853 0
## 16 max fnr 0.999845 0.971251 0
## 17 max fpr 0.000886 1.000000 399
## 18 max tpr 0.000886 1.000000 399
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## Cross-Validation Metrics Summary:
## mean sd cv_1_valid cv_2_valid cv_3_valid cv_4_valid
## accuracy 0.83662486 0.0077441493 0.84442 0.8424983 0.8255325 0.8383588
## auc 0.90377724 0.002580588 0.906217 0.9066376 0.90287024 0.90057284
## aucpr 0.7671483 0.008452881 0.7596798 0.77820635 0.7729425 0.7663389
## err 0.16337517 0.0077441493 0.15557998 0.1575017 0.17446749 0.1616412
## err_count 1183.6 43.672646 1132.0 1160.0 1245.0 1174.0
## cv_5_valid
## accuracy 0.83231455
## auc 0.9025886
## aucpr 0.758574
## err 0.16768546
## err_count 1207.0
##
## ---
## mean sd cv_1_valid cv_2_valid cv_3_valid
## precision 0.64973193 0.023413904 0.65447354 0.67998004 0.6202975
## r2 0.4371859 0.0066841524 0.43763807 0.44788808 0.43729803
## recall 0.7489738 0.031996913 0.7535572 0.72408295 0.7899721
## residual_deviance 4757.921 52.38938 4699.8945 4820.0493 4707.5337
## rmse 0.32420433 0.0019110688 0.32094553 0.32402903 0.3254827
## specificity 0.86562604 0.019410979 0.8733466 0.8831145 0.83748364
## cv_4_valid cv_5_valid
## precision 0.66062176 0.63328695
## r2 0.43105683 0.43204856
## recall 0.71070236 0.7665544
## residual_deviance 4789.476 4772.651
## rmse 0.32529965 0.3252647
## specificity 0.88023406 0.8539513Al imprimir el modelo por pantalla, se muestra toda la información disponible: tipo de modelo, coeficientes de regresión obtenidos, métricas, resultados de la validación cruzada… Esta suele ser más información de la que se necesita en un análisis convencional. Para poder acceder directamente a la información de interés, H2O incorpora una serie de funciones que extraen información concreta del modelo.
Coeficientes de regresión
Si se especifica en la creación del modelo que se estandaricen los predictores (standardize=TRUE), se generan dos tipos de coeficientes de regresión:
Coeficientes estandarizados: son los obtenidos directamente por el modelo al emplear los predictores estandarizados. Es importante tener en cuenta que, al estar estandarizados, su valor no está en las mismas unidades que los predictores originales, lo que dificulta su interpretación. Sin embargo, sí que son útiles (su valor absoluto) para ordenar los predictores de mayor a menor influencia en el modelo.
Coeficientes: son los coeficientes “normales” obtenidos a partir de los coeficientes estandarizados pero revirtiendo la estandarización. Sus unidades son las mismas que las de los predictores originales y pueden ser interpretados de forma normal.
Para obtener más información sobre la correcta interpretación de los coeficientes de regresión de un modelo de regresión logística consultar Regresión logística simple y múltiple en R.
# Coeficientes de regresión de cada uno de los predictores.
as.data.frame(modelo_glm@model$coefficients_table) %>% head()Predictores incluidos en el modelo
Si el modelo GLM incorpora penalización L1 (lasso o elastic net), entonces, el propio ajuste del modelo realiza selección de predictores. Los predictores incluidos en el modelo son aquellos cuyo coeficiente de regresión es distinto de cero.
## [1] "Intercept" "occupation.Craft-repair"
## [3] "occupation.Exec-managerial" "occupation.Farming-fishing"
## [5] "occupation.Handlers-cleaners" "occupation.Machine-op-inspct"
## [7] "occupation.Other-service" "occupation.Priv-house-serv"
## [9] "occupation.Prof-specialty" "occupation.Protective-serv"
## [11] "occupation.Sales" "occupation.Tech-support"
## [13] "occupation.Transport-moving" "workclass.Federal-gov"
## [15] "workclass.Private" "workclass.Self-emp-inc"
## [17] "workclass.Self-emp-not-inc" "workclass.State-gov"
## [19] "education.bachillerato" "education.doctorado"
## [21] "education.master" "education.primaria"
## [23] "education.secundaria" "education.universidad"
## [25] "relationship.Husband" "relationship.Other-relative"
## [27] "relationship.Own-child" "relationship.Unmarried"
## [29] "relationship.Wife" "race.Amer-Indian-Eskimo"
## [31] "race.Asian-Pac-Islander" "race.Other"
## [33] "race.White" "native_country.Europa"
## [35] "native_country.Norte America" "native_country.Otros"
## [37] "native_country.Sud America" "marital_status.casado"
## [39] "marital_status.soltero" "sex.Female"
## [41] "age" "final_weight"
## [43] "education_number" "capital_gain"
## [45] "capital_loss" "hours_per_week"Importancia de los predictores
En los modelos lineales generalizados, la importancia de los predictores puede estudiarse a partir del valor absoluto de los coeficientes de regresión estandarizados.
coeficientes <- as.data.frame(modelo_glm@model$coefficients_table)
# Se excluye el intercept.
coeficientes <- coeficientes %>% filter(names != "Intercept")
# Se calcula el valor absoluto.
coeficientes <- coeficientes %>%
mutate(abs_stand_coef = abs(standardized_coefficients))
# Se añade una variable con el signo del coeficiente.
coeficientes <- coeficientes %>%
mutate(signo = if_else(standardized_coefficients > 0,
"Positivo",
"Negativo"))
ggplot(data = coeficientes,
aes(x = reorder(names, abs_stand_coef),
y = abs_stand_coef,
fill = signo)) +
geom_col() +
coord_flip() +
labs(title = "Importancia de los predictores en el modelo GLM",
x = "Predictor",
y = "Valor absoluto coeficiente estandarizado") +
theme_bw() +
theme(legend.position = "bottom")
Métricas de entrenamiento
Se pueden obtener toda una serie de métricas calculadas a partir de los datos de entrenamiento. Al tratarse de métricas de entrenamiento, no son válidas como estimación del error que comete el modelo al predecir nuevos valores (son demasiado optimistas) ni tampoco para comparar modelos, para esto hay que emplear métricas de validación o de test.
## [1] 0.9045265# Mean Squared Error
h2o.mse(modelo_glm, train = TRUE)
# R2
h2o.r2(modelo_glm, train = TRUE)
# LogLoss
h2o.logloss(modelo_glm, train = TRUE)
# Coeficiente de Gini
h2o.giniCoef(modelo_glm, train = TRUE)
# Desviance del modelo nulo
h2o.null_deviance(modelo_glm, train = TRUE)
# Desviance del modelo final
h2o.residual_deviance(modelo_glm, train = TRUE)
# AIC
h2o.aic(modelo_glm, train = TRUE)Se pueden obtener todas a la vez, junto con la matriz de confusión, con la función h2o.performance().
## H2OBinomialMetrics: glm
## ** Reported on training data. **
##
## MSE: 0.1046908
## RMSE: 0.3235596
## LogLoss: 0.3269983
## Mean Per-Class Error: 0.1890437
## AUC: 0.9045265
## AUCPR: 0.7689901
## Gini: 0.8090531
## R^2: 0.4395605
## Residual Deviance: 23699.53
## AIC: 23791.53
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 23127 4102 0.150648 =4102/27229
## >50K 2049 6960 0.227439 =2049/9009
## Totals 25176 11062 0.169739 =6151/36238
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.322194 0.693538 226
## 2 max f2 0.145909 0.792705 304
## 3 max f0point5 0.585860 0.709997 128
## 4 max accuracy 0.503406 0.848474 158
## 5 max precision 0.996336 0.992322 2
## 6 max recall 0.000953 1.000000 399
## 7 max specificity 0.999858 0.999890 0
## 8 max absolute_mcc 0.401537 0.587917 195
## 9 max min_per_class_accuracy 0.278949 0.816961 246
## 10 max mean_per_class_accuracy 0.244960 0.821457 261
## 11 max tns 0.999858 27226.000000 0
## 12 max fns 0.999858 8750.000000 0
## 13 max fps 0.000953 27229.000000 399
## 14 max tps 0.000953 9009.000000 399
## 15 max tnr 0.999858 0.999890 0
## 16 max fnr 0.999858 0.971251 0
## 17 max fpr 0.000953 1.000000 399
## 18 max tpr 0.000953 1.000000 399
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`Métricas de validación y validación cruzada
Al igual que para el entrenamiento, se pueden mostrar las métricas obtenidas en el conjunto de validación (val=TRUE) o el promedio obtenido mediante validación cruzada (xval=TRUE).
## [1] 0.9036885## H2OBinomialMetrics: glm
## ** Reported on cross-validation data. **
## ** 5-fold cross-validation on training data (Metrics computed for combined holdout predictions) **
##
## MSE: 0.1051063
## RMSE: 0.324201
## LogLoss: 0.3282526
## Mean Per-Class Error: 0.1914324
## AUC: 0.9036885
## AUCPR: 0.7671405
## Gini: 0.8073771
## R^2: 0.4373362
## Residual Deviance: 23790.44
## AIC: 23880.44
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 23278 3951 0.145103 =3951/27229
## >50K 2142 6867 0.237762 =2142/9009
## Totals 25420 10818 0.168138 =6093/36238
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.331430 0.692692 224
## 2 max f2 0.143851 0.791580 306
## 3 max f0point5 0.578472 0.709330 135
## 4 max accuracy 0.480840 0.847922 169
## 5 max precision 0.988455 0.989614 6
## 6 max recall 0.000886 1.000000 399
## 7 max specificity 0.999845 0.999853 0
## 8 max absolute_mcc 0.383540 0.585427 203
## 9 max min_per_class_accuracy 0.278878 0.816628 245
## 10 max mean_per_class_accuracy 0.239945 0.820311 263
## 11 max tns 0.999845 27225.000000 0
## 12 max fns 0.999845 8750.000000 0
## 13 max fps 0.000886 27229.000000 399
## 14 max tps 0.000886 9009.000000 399
## 15 max tnr 0.999845 0.999853 0
## 16 max fnr 0.999845 0.971251 0
## 17 max fpr 0.000886 1.000000 399
## 18 max tpr 0.000886 1.000000 399
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`Predicciones
Una vez que el modelo ha sido entrenado, puede emplearse para predecir nuevas observaciones con la función h2o.predict(), que recibe como argumentos: un modelo y un nuevo set de datos.
El resultado devuelto por la función h2o.predict() para un modelo GLM logístico es una tabla con 3 columnas, una con la clase predicha y otras dos con la probabilidad de pertenecer a cada una de las clases. El orden de las filas es el mismo que el de los datos pasados al argumento newdata.
Si el nuevo set de datos incluye la variable respuesta, se pueden calcular métricas que cuantifican el grado de acierto.
## H2OBinomialMetrics: glm
##
## MSE: 0.1027842
## RMSE: 0.3205997
## LogLoss: 0.3209669
## Mean Per-Class Error: 0.1935457
## AUC: 0.9058415
## AUCPR: 0.7746252
## Gini: 0.8116829
## R^2: 0.4439796
## Residual Deviance: 5767.134
## AIC: 5859.134
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 5942 843 0.124245 =843/6785
## >50K 578 1621 0.262847 =578/2199
## Totals 6520 2464 0.158170 =1421/8984
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.365986 0.695261 212
## 2 max f2 0.129705 0.790623 309
## 3 max f0point5 0.670778 0.714700 109
## 4 max accuracy 0.511600 0.851514 162
## 5 max precision 0.999907 1.000000 0
## 6 max recall 0.004048 1.000000 394
## 7 max specificity 0.999907 1.000000 0
## 8 max absolute_mcc 0.365986 0.590662 212
## 9 max min_per_class_accuracy 0.278373 0.817539 246
## 10 max mean_per_class_accuracy 0.259330 0.820227 254
## 11 max tns 0.999907 6785.000000 0
## 12 max fns 0.999907 2148.000000 0
## 13 max fps 0.000687 6785.000000 399
## 14 max tps 0.004048 2199.000000 394
## 15 max tnr 0.999907 1.000000 0
## 16 max fnr 0.999907 0.976808 0
## 17 max fpr 0.000687 1.000000 399
## 18 max tpr 0.004048 1.000000 394
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`# Cálculo manual de accuracy
mean(as.vector(predicciones$predict) == as.vector(datos_test_h2o$salario))## [1] 0.8285841Grid search
En el modelo anterior, se ha hecho una búsqueda del valor óptimo de lambda pero no de alpha (se ha empleado un valor fijo de alpha = 0.95). A continuación, se repite el ajuste del modelo, pero esta vez, comparando también distintos valores de alpha. Para estimar la capacidad predictiva de cada modelo se emplea validación cruzada con 10 particiones.
Nota: con el número de observaciones disponible, emplear validación simple sería suficiente y los resultados se obtendrían mucho más rápido.
# Valores de alpha que se van a comparar.
hiperparametros <- list(alpha = c(0, 0.1, 0.5, 0.95, 1))
grid_glm <- h2o.grid(
# Algoritmo y parámetros
algorithm = "glm",
family = "binomial",
link = "logit",
# 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 = 10,
# Reparto estratificado de las observaciones en la creación
# de las particiones
fold_assignment = "Stratified",
keep_cross_validation_predictions = FALSE,
grid_id = "grid_glm"
)# Se muestran los modelos ordenados de mayor a menor AUC.
resultados_grid <- h2o.getGrid(
grid_id = "grid_glm",
sort_by = "auc",
decreasing = TRUE
)
print(resultados_grid)## H2O Grid Details
## ================
##
## Grid ID: grid_glm
## Used hyper parameters:
## - alpha
## Number of models: 5
## Number of failed models: 0
##
## Hyper-Parameter Search Summary: ordered by decreasing auc
## alpha model_ids auc
## 1 [0.0] grid_glm_model_1 0.9036504646332404
## 2 [0.1] grid_glm_model_2 0.9035758863699662
## 3 [0.5] grid_glm_model_3 0.9035081566940981
## 4 [1.0] grid_glm_model_5 0.903469539629516
## 5 [0.95] grid_glm_model_4 0.9034615496108757De entre los valores comparados, el mayor auc promedio de validación cruzada se obtiene con alpha = 0.0 (modelo ridge), aunque para este problema, la diferencia es mínima.
Para conocer con más detalle la distribución de auc obtenido por cada modelo en cada partición de la validación cruzada, se tiene que extraer la información de cada uno de los modelos.
# Identificador de los modelos creados por validación cruzada.
id_modelos <- unlist(resultados_grid@model_ids)
# Se crea una lista donde se almacenarán los resultados.
auc_xvalidacion <- vector(mode = "list", length = length(id_modelos))
# Se recorre cada modelo almacenado en el grid y se extraen la métrica (auc)
# obtenida en cada partición.
for (i in seq_along(id_modelos)) {
modelo <- h2o.getModel(resultados_grid@model_ids[[i]])
metricas_xvalidacion_modelo <- modelo@model$cross_validation_metrics_summary
names(auc_xvalidacion)[i] <- modelo@model$model_summary$regularization
auc_xvalidacion[[i]] <- as.numeric(metricas_xvalidacion_modelo["auc", -c(1,2)])
}
# Se eliminan los espacios en blanco del nombre de los modelos.
library(stringr)
names(auc_xvalidacion) <- str_remove_all(string = names(auc_xvalidacion),
pattern = "[ )=]")
names(auc_xvalidacion) <- str_replace_all(string = names(auc_xvalidacion),
pattern = "[(,]",
replacement = "_")
# Se convierte la lista en dataframe.
auc_xvalidacion_df <- as.data.frame(auc_xvalidacion) %>%
mutate(resample = row_number()) %>%
gather(key = "modelo", value = "auc", -resample) %>%
mutate(modelo = str_replace_all(string = modelo,
pattern = "_" ,
replacement = " \n "))
# Gráfico
ggplot(data = auc_xvalidacion_df, aes(x = modelo, y = auc, color = modelo)) +
geom_boxplot(alpha = 0.6, outlier.shape = NA) +
geom_jitter(width = 0.1, alpha = 0.6) +
stat_summary(fun.y = "mean", colour = "red", size = 2, geom = "point") +
theme_bw() +
labs(title = "Accuracy obtenido por 10-CV") +
coord_flip() +
theme(legend.position = "none")
En este caso, los resultados de los 5 modelos son prácticamente idénticos.
Una vez identificado el mejor modelo, mediante h2o.grid(), se extrae del objeto grid y se almacena por separado.
# Se extrae el mejor modelo, que es el que ocupa la primera posición tras haber
# ordenado los resultados de mayor a menor AUC (métrica utilizada en este caso).
modelo_glm_final <- h2o.getModel(resultados_grid@model_ids[[1]])Gradient Boosting Machine (GBM)
Introducción
Boosting Machine (BM) es un tipo de modelo obtenido al combinar (ensembling) múltiples modelos sencillos, también conocidos como weak learners. Está combinación se realiza de manera secuencial de tal forma que, cada nuevo modelo que se incorpora al conjunto, intenta corregir los errores de los anteriores modelos. Como resultado de la combinación de múltiples modelos, Boosting Machine consigue aprender relaciones no lineales entre la variable respuesta y los predictores. Si bien los modelos combinados pueden ser muy variados, H2O, al igual que la mayoría de librerías, utiliza como weak learners modelos basados en árboles.
Gradient Boosting Machine (GBM) es una generalización del modelo de Boosting Machine que permite aplicar el método de descenso de gradiente para optimizar cualquier función de coste durante el ajuste del modelo.
El valor predicho por un modelo GBM es la agregación (normalmente la moda en problemas de clasificación y la media en problemas de regresión) de las predicciones de todos los modelos individuales que forman el ensemble.
Dado que el ajuste de los weak learners que forman un GBM se hace de forma secuencial (el árbol i se ajusta a partir de los resultados del árbol i-1), las posibilidades de paralelizar el algoritmo son limitadas. H2O consigue acelerar el proceso en cierta medida paralelizando el ajuste de los árboles individuales.
Para información más detallada sobre métodos boosting consultar Árboles de predicción.
Detención temprana del entrenamiento (Early Stopping)
Una de las características de GBM es que, con el número suficiente de weak learners, el modelo final tiende a ajustarse perfectamente a los datos de entrenamiento causando overfitting. Este comportamiento implica que el analista tiene que encontrar el número adecuado de árboles y, para ello, suele tener que entrenar el modelo con cientos o miles de árboles hasta identificar el momento en el que empieza el overfitting. Esto suele ser poco eficiente en términos de tiempo, ya que, posiblemente, se estén ajustando muchos árboles innecesarios.
De nuevo, la experiencia de los creadores de H2O queda reflejada en su herramienta, esta vez incluyendo toda una serie de estrategias para detener el proceso de ajuste de un modelo GBM a partir del momento en el que este deja de mejorar. Por defecto, la métrica empleada se calcula con los datos de validación. Cuatro argumentos controlan la estrategia de parada:
stopping_metric: métrica empleada para cuantificar cuánto mejora el modelo.score_tree_interval: número de árboles tras los que se evalúa el modelo. Por ejemplo, si el valor es 10, entonces, tras cada 10 nuevos árboles que se añaden al modelo, se calcula lastopping_metric. Si bien la evaluación del modelo puede hacerse tras cada nuevo árbol que se añade, esto puede ralentizar el entrenamiento.stopping_tolerance: porcentaje mínimo de mejora entre dos mediciones consecutivas (score_tree_interval) por debajo del cual se considera que el modelo no ha mejorado.stopping_rounds: número de mediciones consecutivas en las que no se debe superar elstopping_tolerancepara que el algoritmo se detenga.
Algunos ejemplos de criterios de parada son:
Detener el entrenamiento de un modelo cuando el error de clasificación no se reduce más de un 1% durante dos mediciones consecutivas:
stopping_rounds=1,stopping_tolerance=0.01ystopping_metric="misclassification".Detener el entrenamiento de un modelo cuando el logloss no se reduce nada durante 3 mediciones consecutivas:
stopping_rounds=3,stopping_tolerance=0ystopping_metric="logloss".Detener el entrenamiento de un modelo cuando la media de AUC no aumenta más de un 0.1% durante cinco mediciones consecutivas:
stopping_rounds=5,stopping_tolerance=0.001ystopping_metric="AUC".Para detener el modelo aunque la métrica converja, se tiene que deshabilitar este comportamiento de parada automática indicando
stopping_rounds=0. En este caso, el algoritmo se detendrá cuando se alcance otra condición de parada, por ejemplo, el número máximo de árboles que pueden formar el modelo.
También es posible detener el entrenamiento del modelo cuando se supera un tiempo máximo especificado con el argumento max_runtimesecs > 0.
En la práctica, es conveniente establecer un número máximo de árboles muy elevado (potencialmente por encima de la cantidad necesaria) y activar la parada temprana.
Comportamiento estocástico
Stochastic Gradient Boosting es el resultado de emplear, en el entrenamiento de cada árbol, únicamente una selección aleatoria de observaciones y predictores del conjunto de entrenamiento. Esta modificación suele reducir el overfitting del modelo y mejorar su generalización. Se puede especificar el % de datos aleatorios seleccionados en cada ajuste con los argumentos:
sample_rate: porcentaje de observaciones aleatorias empleadas para el ajuste de cada árbol del ensemble.sample_rate_per_class: cuando el set de datos con el que se trabaja está desbalanceado (el % de observaciones de cada grupo es muy distinto), al hacer una selección aleatoria de observaciones para el entrenamiento de los árboles, puede ocurrir que, en algunos árboles, la selección excluya a alguno de los grupos. Para evitarlo, con el argumentosample_rate_per_classse indica que la selección debe hacerse de forma estratificada, asegurando que todos los grupos queden representados.col_sample_rate: porcentaje de predictores (columnas) seleccionados de forma aleatoria para el ajuste de cada árbol del ensemble.
Los tres parámetros deben tener un valor entre 0 y 1. Por defecto, el valor empleado es 1, es decir, se emplean todos los datos y no hay comportamiento estocástico.
Distribución y función de coste
Dependiendo el tipo de modelo (regresión, clasificación binaria, clasificación múltiple…) se emplea una distribución y una función de coste distinta. Indicando la distribución con el argumento distribution, H2O selecciona automáticamente la función de coste correspondiente.
Las distribuciones disponibles son: bernoulli, quasibinomial, multinomial, poisson, laplace, tweedie, gaussian, huber, gamma y quantile. De entre todas ellas, las más empleadas suelen ser: gaussian para regresión, bernoulli para clasificación binaria y multinomial para clasificación multiclase. Aunque menos frecuentes, también pueden resultar muy útiles la distribución huber, cuando se quiere que los errores penalicen más que su valor absoluto pero menos que el cuadrado, y la distribución quantile cuando en lugar de la media se quiere predecir uno de los cuantiles, por ejemplo la mediana.
Otros hiperparámetros
Además de los ya descritos (criterio de parada, comportamiento estocástico, distribución), el modelo GBM de H2O incorpora muchos otros parámetros e hiperparámetros con los que se puede controlar el comportamiento del algoritmo. A continuación se describen algunos de los que tienen más influencia en el modelo final. Se puede encontrar un listado completo en la documentación.
Complejidad de los árboles
ntrees: número de árboles que forman el ensemble.max_depth: profundidad máxima que pueden tener los árboles (número máximo de nodos).min_rows: número mínimo de observaciones que debe tener cada nodo.min_split_improvement: reducción mínima en el error cuadrático que debe de conseguir la división de un nodo para que se lleve a cabo.
Velocidad de aprendizaje
learn_rate: determina cuanto influye cada árbol individual en el conjunto del ensemble. Cuanto más pequeño es el valor de learning rate, más lentamente aprende el modelo y más árboles se necesitan, sin embargo, al aprender poco a poco, se reduce el riesgo de overfitting. Suelen emplearse valores inferiores a 0.1.learn_rate_annealing: emplear un valor delearn_ratemuy pequeño puede conllevar que el modelo necesite varios miles de árboles y, por lo tanto, el entrenamiento tarde mucho tiempo. En lugar de utilizar un valor muy bajo desde el principio, se puede iniciar el proceso con un valor moderadamente alto e ir reduciéndolo tras cada nuevo árbol que se incorpora al modelo. Ellearn_rate_annealingdetermina el factor por el que se reduce ellearn_rateen cada iteración. Para N árboles, GBM empieza conlearn_ratey termina conlearn_rate*learn_rate_annealing^N. Por ejemplo, en lugar de emplearearn_rate=0.01, se puede emplear unlearn_rate=0.05y unlearn_rate_annealing=0.99. Con la segunda opción, el algoritmo converge mucho más rápido con aproximadamente la misma capacidad predictiva final.
Comportamiento monotónico
Los modelos GBM son capaces de aprender interacciones complejas entre los predictores y la variable respuesta. Esto puede convertirse en una desventaja cuando, debido al ruido de los datos, el modelo aprende comportamientos que no reflejan la realidad. Con el argumento monotone_constraints se puede forzar a que la relación de cada predictor (numérico) con la variable respuesta tenga siempre la misma dirección (positiva o negativa). Puede encontrarse una descripción más detallada en GBM monotónicos.
Para conseguir maximizar la capacidad predictiva de un modelo GBM, es imprescindible entender bien cómo influye cada uno de los parámetros e hiperparámetros. Algunos de ellos pueden encontrarse descritos con más detalle en el documento Boosting.
Modelo
A continuación, se ajusta un modelo GBM que prediga la variable respuesta salario. En esta ocasión, se emplea validación simple en lugar de validación cruzada, por lo que dividen de nuevo los datos, esta vez, en un conjunto de entrenamiento (60%), validación (20%) y test (20%).
particiones <- h2o.splitFrame(data = datos_h2o, ratios = c(0.6, 0.20),
seed = 123)
datos_train_h2o <- h2o.assign(data = particiones[[1]], key = "datos_train_H2O")
datos_val_h2o <- h2o.assign(data = particiones[[2]], key = "datos_val_H2O")
datos_test_h2o <- h2o.assign(data = particiones[[3]], key = "datos_test_H2O")Se entrena un modelo inicial en el que únicamente se especifican el número máximo de árboles, la complejidad de los mismos y la estrategia de parada temprana, el resto de valores se dejan por defecto.
modelo_gbm <- h2o.gbm(
# Tipo de distribución (clasificación binaria)
distribution = "bernoulli",
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Datos de validación para estimar el error
validation_frame = datos_val_h2o,
# Número de árboles
ntrees = 500,
# Complejidad de los árboles
max_depth = 3,
min_rows = 10,
# Aprendizaje
learn_rate = 0.01,
# Detención temprana
score_tree_interval = 5,
stopping_rounds = 3,
stopping_metric = "AUC",
stopping_tolerance = 0.001,
model_id = "modelo_gbm",
seed = 123
)## Model Details:
## ==============
##
## H2OBinomialModel: gbm
## Model ID: modelo_gbm
## Model Summary:
## number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1 95 95 15088 3
## max_depth mean_depth min_leaves max_leaves mean_leaves
## 1 3 3.00000 8 8 8.00000
##
##
## H2OBinomialMetrics: gbm
## ** Reported on training data. **
##
## MSE: 0.1244416
## RMSE: 0.3527628
## LogLoss: 0.4020413
## Mean Per-Class Error: 0.2194013
## AUC: 0.8896691
## AUCPR: 0.7635744
## Gini: 0.7793381
## R^2: 0.3359533
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 18247 2157 0.105715 =2157/20404
## >50K 2263 4531 0.333088 =2263/6794
## Totals 20510 6688 0.162512 =4420/27198
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.340967 0.672155 117
## 2 max f2 0.151212 0.773090 175
## 3 max f0point5 0.469619 0.711687 60
## 4 max accuracy 0.380419 0.846717 85
## 5 max precision 0.683568 0.995885 1
## 6 max recall 0.118008 1.000000 181
## 7 max specificity 0.700388 0.999951 0
## 8 max absolute_mcc 0.367128 0.566731 93
## 9 max min_per_class_accuracy 0.263324 0.784074 156
## 10 max mean_per_class_accuracy 0.278211 0.798695 149
## 11 max tns 0.700388 20403.000000 0
## 12 max fns 0.700388 6575.000000 0
## 13 max fps 0.118008 20404.000000 181
## 14 max tps 0.118008 6794.000000 181
## 15 max tnr 0.700388 0.999951 0
## 16 max fnr 0.700388 0.967766 0
## 17 max fpr 0.118008 1.000000 181
## 18 max tpr 0.118008 1.000000 181
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`
## H2OBinomialMetrics: gbm
## ** Reported on validation data. **
##
## MSE: 0.1206073
## RMSE: 0.3472855
## LogLoss: 0.3929575
## Mean Per-Class Error: 0.2051175
## AUC: 0.898275
## AUCPR: 0.7743462
## Gini: 0.7965499
## R^2: 0.3480206
##
## Confusion Matrix (vertical: actual; across: predicted) for F1-optimal threshold:
## <=50K >50K Error Rate
## <=50K 6105 720 0.105495 =720/6825
## >50K 675 1540 0.304740 =675/2215
## Totals 6780 2260 0.154314 =1395/9040
##
## Maximum Metrics: Maximum metrics at their respective thresholds
## metric threshold value idx
## 1 max f1 0.339734 0.688268 89
## 2 max f2 0.157774 0.773203 136
## 3 max f0point5 0.379802 0.727403 64
## 4 max accuracy 0.379802 0.857190 64
## 5 max precision 0.700388 1.000000 0
## 6 max recall 0.118008 1.000000 143
## 7 max specificity 0.700388 1.000000 0
## 8 max absolute_mcc 0.379802 0.589602 64
## 9 max min_per_class_accuracy 0.265986 0.801806 119
## 10 max mean_per_class_accuracy 0.281448 0.809153 114
## 11 max tns 0.700388 6825.000000 0
## 12 max fns 0.700388 2144.000000 0
## 13 max fps 0.118008 6825.000000 143
## 14 max tps 0.118008 2215.000000 143
## 15 max tnr 0.700388 1.000000 0
## 16 max fnr 0.700388 0.967946 0
## 17 max fpr 0.118008 1.000000 143
## 18 max tpr 0.118008 1.000000 143
##
## Gains/Lift Table: Extract with `h2o.gainsLift(<model>, <data>)` or `h2o.gainsLift(<model>, valid=<T/F>, xval=<T/F>)`Evolución del modelo
Tras ajustar un modelo GBM, es fundamental evaluar la como aprende el modelo a medida que se añaden nuevos árboles al ensemble. Una forma de hacerlo es representar la evolución de una métrica (en este caso el área bajo la curva roc, AUC) de entrenamiento y validación.
# H2O almacena las métricas de entrenamiento y test bajo el nombre de scoring.
# Los valores se encuentran almacenados dentro del modelo.
scoring <- as.data.frame(modelo_gbm@model$scoring_history)
head(scoring)Tal y como se ha indicado en el argumento (score_tree_interval=5), cada 5 árboles se han recalculado las métricas de rendimiento en el conjunto de entrenamiento y validación. El argumento stopping_metric="AUC" indica que se emplee el área bajo la curva como criterio para decidir si el algoritmo debe detenerse de forma temprana, pero se calculan siempre todas (AUC, logloss, lift, rmse, classification_error).
scoring <- scoring %>%
select(-timestamp, -duration) %>%
gather(key = "metrica", value = "valor", -number_of_trees) %>%
separate(col = metrica, into = c("conjunto", "metrica"), sep = "_")
scoring %>% filter(metrica == "auc") %>%
ggplot(aes(x = number_of_trees, y = valor, color = conjunto)) +
geom_line() +
labs(x = "AUC", y = "número de árboles (ntrees)",
title = "Evolución del AUC vs número de árboles") +
theme_bw()
En el gráfico puede observarse que no se han alcanzado los 500 árboles indicados en el modelo, esto es así porque, alcanzados los 95 árboles, el valor de AUC se ha estabilizado lo suficiente como para que se ejecute la detención temprana.
Importancia de los predictores
En los modelos GBM, se puede estudiar la influencia de los predictores cuantificando la reducción total de error cuadrático que ha conseguido cada predictor en el conjunto de todos los árboles que forman el modelo.
# Se extraen los valores de importancia
importancia <- as.data.frame(modelo_gbm@model$variable_importances)
importanciaggplot(data = importancia,
aes(x = reorder(variable, scaled_importance), y = scaled_importance)) +
geom_col() +
coord_flip() +
labs(title = "Importancia de los predictores en el modelo GBM",
subtitle = "Importancia en base a la reducción del error cuadrático medio",
x = "Predictor",
y = "Importancia relativa") +
theme_bw()
Predicciones y error
Una vez ajustado el modelo, se puede predecir nuevas observaciones y estimar su error promedio.
## [1] 0.8928113Grid search
Empleando los valores por defecto de los hiperparámetros, se ha conseguido un AUC de test de 0.893. A continuación, se intenta encontrar los valores óptimos de los hiperparámetros con los que se consiga mejorar aún más el AUC.
Dada la gran cantidad de hiperparámetros que tiene el algoritmo GBM y, aunque todos interaccionan entre sí, no es posible explorar todo el espacio de posibles valores a la vez. Una estrategia que suele funcionar bien es establecer primero los hiperparámetros más influyentes (número de árboles, learning rate, complejidad de los árboles, % de observaciones empleadas en cada ajuste).
Número de árboles (
ntree): incorporar tantos árboles como sea necesario hasta que el error de validación se estabilice o empiece a aumentar. Dado que H2O incorpora estrategias de parada temprana, no hay problema en indicar un número elevado porque, de no ser necesario, no se ajustarán.Tasa de aprendizaje (
learning_rate): por lo general, valores pequeños (<0.01) dan lugar a en mejores modelos, pero requieren más árboles para aprender. Se recomienda reducir el valor tanto como permitan las limitaciones de tiempo/computación.Profundidad de los árboles (
max_depth): es la principal forma de controlar la complejidad de los árboles que forman el ensemble. Al tratarse de weak learners, la complejidad de los árboles debe ser baja, aunque depende mucho de los datos, el valor óptimo raramente está fuera del rango [1, 20].Comportamiento estocástico (
sample_rate,col_sample_rate): la finalidad de emplear solo un subconjunto de observaciones y/o columnas en el entrenamiento de cada árbol es conseguir mejor generalización, sin embargo, puede tener el resultado contrario cuando se dispone de muy pocos datos. El valor óptimo suele estar dentro del rango de [0.7-0.9]. En situaciones en las que la variable respuesta es cualitativa y los grupos están desbalanceados, es conveniente realizar una selección aleatoria estratificada, este comportamiento puede activarse consample_rate_per_class.
# Hiperparámetros que se quieren comparar.
hiperparametros <- list(
learn_rate = c(0.001, 0.05, 0.1),
max_depth = c(5, 10, 15, 20),
sample_rate = c(0.8, 1)
)
grid_gbm <- h2o.grid(
# Algoritmo
algorithm = "gbm",
distribution = "bernoulli",
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Datos de validación
validation_frame = datos_val_h2o,
# Preprocesado
ignore_const_cols = TRUE,
# Detención temprana
score_tree_interval = 100,
stopping_rounds = 3,
stopping_metric = "AUC",
stopping_tolerance = 0.001,
# Hiperparámetros fijados
ntrees = 5000,
min_rows = 10,
# Hiperparámetros optimizados
hyper_params = hiperparametros,
# Tipo de búsqueda
search_criteria = list(strategy = "Cartesian"),
seed = 123,
grid_id = "grid_gbm"
)Se muestran los modelos ordenados de mayor a menor AUC.
resultados_grid <- h2o.getGrid(
grid_id = "grid_gbm",
sort_by = "auc",
decreasing = TRUE
)
data.frame(resultados_grid@summary_table) %>% select(-model_ids)Los 3 modelos con mayor AUC de validación tienen todos un learning rate de 0.1 y un max depth de 15 y 20, el sample rate parece no influir de forma sustancial. Teniendo en cuenta estos resultados, pude concluirse que los valores óptimos son: learn_rate=0.1, max_depth=20 y sample_rate=0.8.
Se extrae del grid el mejor modelo y se calcula su métrica de test.
modelo_gbm_final <- h2o.getModel(resultados_grid@model_ids[[1]])
# AUC de test
h2o.performance(model = modelo_gbm_final, newdata = datos_test_h2o)@metrics$AUC## [1] 0.9283551Optimizando los hiperparámetros más influyentes, se ha conseguido mejorar el AUC de test de 0.893 a 0.914. Por último, se intenta encontrar el valor óptimo del resto de hiperparámetros, esta vez mediante una búsqueda aleatoria.
# Hiperparámetros que se quieren optimizar mediante búsqueda aleatoria.
# Para realizar esta búsqueda se tiene que pasar un vector de posibles valores
# de cada hiperparámetro, entre los que se escoge aleatoriamente.
hiperparametros <- list(
min_rows = seq(from = 5, to = 50, by = 10),
nbins = 2^seq(4, 10, 1),
nbins_cats = 2^seq(4, 12, 1),
min_split_improvement = c(0, 1e-8, 1e-6, 1e-4),
histogram_type = c("UniformAdaptive","QuantilesGlobal","RoundRobin")
)
# Al ser una búsqueda aleatoria, hay que indicar criterios de parada.
search_criteria <- list(
strategy = "RandomDiscrete",
# Tiempo máximo de búsqueda (5 minutos)
max_runtime_secs = 5*60,
# Número máximo de modelos
max_models = 100,
# Reproducible
seed = 1234
)
grid_gbm2 <- h2o.grid(
# Algoritmo
algorithm = "gbm",
distribution = "bernoulli",
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Datos de validación
validation_frame = datos_val_h2o,
# Preprocesado
ignore_const_cols = TRUE,
# Detención temprana
score_tree_interval = 100,
stopping_rounds = 3,
stopping_metric = "AUC",
stopping_tolerance = 0.001,
# Hiperparámetros fijados
ntrees = 5000,
learn_rate = 0.1,
max_depth = 20,
sample_rate = 0.8,
# Hiperparámetros optimizados
hyper_params = hiperparametros,
# Tipo de búsqueda
search_criteria = search_criteria,
seed = 123,
grid_id = "grid_gbm2"
)# Se muestran los modelos ordenados de mayor a menor AUC.
resultados_grid <- h2o.getGrid(
grid_id = "grid_gbm2",
sort_by = "auc",
decreasing = TRUE
)
data.frame(resultados_grid@summary_table) %>% select(-model_ids)modelo_gbm_final <- h2o.getModel(resultados_grid@model_ids[[1]])
# AUC de test
h2o.performance(model = modelo_gbm_final, newdata = datos_test_h2o)@metrics$AUC## [1] 0.8970854El AUC de test obtenido es prácticamente idéntico al del anterior modelo.
Deep Learning (Neural Networks)
Introducción
El término deep learning engloba a todo un conjunto de modelos basados en redes neuronales artificiales (artificial neural networks) que contienen múltiples capas intermedias (ocultas). En concreto, H2O incorpora redes neuronales de tipo Multi-layer - feedforward - neural networks. Estas redes se caracterizan por:
Pueden tener una o múltiples capas intermedias, también conocidas como capas ocultas. Son las capas de neuronas que hay entre la capa de entrada y la de salida.
La estructura es full conected, lo que significa que cada neurona está conectada con todas las neuronas de la capa siguiente.
El peso (participación de cada neurona en la red) se aprende mediante el proceso feed-forward. A grandes rasgos, la estrategia de aprendizaje consiste en un proceso iterativo. En cada iteración (época), la red predice las observaciones de entrenamiento, se calcula el error de estas predicciones comparando los resultados frente al valor real y se reajustan los pesos de las neuronas con el objetivo de minimizar el error (backpropagation y descenso de gradiente).
Estas redes son unas de las más empleadas, pero existen otros tipos, por ejemplo, las Convolutional Neural Networks (CNNs) son las que mejores resultados dan cuando se trabaja con imágenes.
Comprender los fundamentos de las redes neuronales y deep learning requiere de un estudio notable, puede encontrarse mucha de la información necesaria en el curso gratuito de Stanford Convolutional Neural Networks for Visual Recognition y en el libro Deep Learning An MIT Press book.
Flexibilidad de las redes neuronales
Los modelos basados en redes neuronales tienen la característica de ser potencialmente capaces de aprender cualquier función. La velocidad con la que aprenden está influenciada principalmente (aunque hay muchos otros hiperparámetros) por el número de observaciones, el número de neuronas y el número de iteraciones de aprendizaje (comúnmente llamadas epoch). A continuación, se muestra un ejemplo sencillo que ilustra este concepto.
Supóngase las siguientes observaciones generadas a partir de un modelo no lineal de espiral.
# Número de observaciones por clase
N <- 300
# Número de dimensiones
D <- 2
# Número clases
K <- 3
# Matriz para almacenar las observaciones
x_1 <- vector(mode = "numeric")
x_2 <- vector(mode = "numeric")
y <- vector(mode = "numeric")
# Simulación de los datos
for(i in 1:K){
set.seed(123)
r <- seq(from = 0, to = 1, length.out = N)
t <- seq(from = i*4, to = (i+1)*4, length.out = N) + rnorm(n = N) * 0.35
x_1 <- c(x_1, r * sin(t))
x_2 <- c(x_2, r*cos(t))
y <- c(y, rep(letters[i], N))
}
datos_espiral <- data.frame(y = as.factor(y), x_1, x_2)
ggplot(data = datos_espiral, aes(x = x_1, y = x_2, color = y)) +
geom_point() +
theme_bw() +
theme(legend.position = "none",
text = element_blank(),
axis.ticks = element_blank())
Se comparan 3 redes neuronales: 10 neuronas y 50 epochs, 20 neuronas y 1000 epochs, 2 capas de 200 neuronas neuronas y 1000 epochs.
datos_espiral_h2o <- as.h2o(datos_espiral)
modelo_dl_10 <- h2o.deeplearning(
x = c("x_1", "x_2"),
y = "y",
distribution = "multinomial",
training_frame = datos_espiral_h2o,
standardize = TRUE,
activation = "Rectifier",
hidden = 10,
stopping_rounds = 0,
epochs = 50,
seed = 123,
model_id = "modelo_dl_10"
)
modelo_dl_20 <- h2o.deeplearning(
x = c("x_1", "x_2"),
y = "y",
distribution = "multinomial",
training_frame = datos_espiral_h2o,
standardize = TRUE,
activation = "Rectifier",
hidden = 20,
stopping_rounds = 0,
epochs = 1000,
seed = 123,
model_id = "modelo_dl_20"
)
modelo_dl_200_200 <- h2o.deeplearning(
x = c("x_1", "x_2"),
y = "y",
distribution = "multinomial",
training_frame = datos_espiral_h2o,
standardize = TRUE,
activation = "Rectifier",
hidden = c(200, 200),
stopping_rounds = 0,
epochs = 1000,
seed = 123,
model_id = "modelo_dl_200_200"
)Una vez entrenados los modelos, se representan sus fronteras de clasificación.
grid_predicciones <- expand.grid(x_1 = seq(from = -1, to = 1, length = 75),
x_2 = seq(from = -1, to = 1, length = 75))
grid_predicciones_h2o <- as.h2o(grid_predicciones)
predicciones_10 <- h2o.predict(object = modelo_dl_10,
newdata = grid_predicciones_h2o)
grid_predicciones$y_10 <- as.vector(predicciones_10$predict)
predicciones_20 <- h2o.predict(object = modelo_dl_20,
newdata = grid_predicciones_h2o)
grid_predicciones$y_20 <- as.vector(predicciones_20$predict)
predicciones_200_200 <- h2o.predict(object = modelo_dl_200_200,
newdata = grid_predicciones_h2o)
grid_predicciones$y_200_200 <- as.vector(predicciones_200_200$predict)p1 <- ggplot(data = grid_predicciones, aes(x = x_1, y = x_2, color = y_10)) +
geom_point(size = 0.5) +
theme_bw() +
labs(title = "1 capa oculta, 10 neuronas, 50 epochs") +
theme(legend.position = "none",
plot.title = element_text(size=11),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
p2 <- ggplot(data = grid_predicciones, aes(x = x_1, y = x_2, color = y_20)) +
geom_point(size = 0.5) +
labs(title = "1 capa oculta, 20 neuronas, 1000 epochs") +
theme_bw() +
theme(legend.position = "none",
plot.title = element_text(size=11),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
p3 <- ggplot(data = grid_predicciones, aes(x = x_1, y = x_2, color = y_200_200)) +
geom_point(size = 0.5) +
labs(title = "2 capas ocultas, 200 neuronas, 1000 epochs") +
theme_bw() +
theme(legend.position = "none",
plot.title = element_text(size=11),
axis.text = element_blank(),
axis.title = element_blank(),
axis.ticks = element_blank())
ggarrange(p1, p2, p3, nrow = 1)
Este ejemplo debe servir como reflexión sobre la alta capacidad que tienen los modelos basados en redes para adaptarse a los datos de entrenamiento, por lo tanto, cuando la complejidad de los datos no es muy grande o se dispone de pocas observaciones, se debe tener muy presente el riesgo de overfitting.
Parámetros e hiperparámetros
Los modelos de Deep Learning ofrecidos por H2O tienen un número muy elevado de parámetros configurables (listado). Para la gran mayoría de casos, los valores por defecto dan buenos resultados, sin embargo, es conveniente conocer, al menos, los más influyentes.
Arquitectura
hidden: número de capas ocultas y número de neuronas por capa. La estructura se define mediante un vector de números enteros. Por ejemplo, una red con una única capa de 100 neuronas se crea conhidden = c(100), y una red con dos capas ocultas de 10 neuronas cada una conhidden = c(10, 10). Las capas de entrada y salida se crean automáticamente. La capa de entrada tiene tantas neuronas como predictores (tras codificar las variables categóricas). La capa de salida tiene una única neurona en problemas de regresión, y tantas neuronas como clases en los problemas de clasificación.
Preprocesado
standardize: por defecto, se estandarizan los valores para que tengan media cero y varianza uno. Los modelos basados en redes no son independientes de la escala de los predictores, por lo que es fundamental estandarizarlos.missing_values_handling: las redes neuronales no aceptan observaciones con valores ausentes. H2O ofrece dos opciones: excluirlosSkipo imputarlos con la mediaMeanImputation.shuffle_training_data: mezclar aleatoriamente las observaciones antes de entrenar el modelo.
Aprendizaje
activation: función de activación de las neuronas de las capas intermedias (Tahn, Tahn with dropout, Rectifier, Rectifier with dropout, Maxout, Maxout with dropout). Las funciones Rectifier y Rectifier with dropout suelen dar buenos resultados para un abanico amplio de situaciones. La función de activación de la capa de salida se selecciona automáticamente dependiendo de si es un problema de regresión o clasificación.loss: función de coste que se intenta minimizar durante el entrenamiento de la red. Esta es la forma en que se cuantifica el error que comete la red y en función de ello aprende, por lo tanto, es importante escoger la función de coste más adecuada para el problema en cuestión (esto es así para todos los algoritmos, no solo para redes). Por ejemplo, en problemas de regresión, las dos funciones más comunes son el error absoluto y el error cuadrático. Con el primero, el modelo resultante prioriza predecir bien la mayoría de observaciones, aunque para unas pocas se equivoque por mucho. Con el segundo, como los errores se elevan al cuadrado, el modelo resultante tiende a evitar grandes errores a cambio de aumentar un poco el error en la mayoría de observaciones.epochs: número de iteraciones de aprendizaje durante el entrenamiento de la red. Encontrar el valor óptimo de épocas es muy importante ya que, si son muy pocas, la red puede no aprender lo suficiente, pero si son demasiadas, se produce overfiting. Normalmente, se puede tener una idea aproximada del valor correcto en función de la complejidad de la red (cuanto más compleja más pesos hay que aprender) y el número de observaciones de entrenamiento. Gracias al sistema de early stopping implementado en H2O, se puede indicar un número muy elevado de épocas sin riesgo de overfitting, ya que el entrenamiento se detendrá cuando la métrica de validación elegida no mejore. Aun así, siempre es recomendable representar la evolución del error de entrenamiento y validación en función del número de épocas.adaptive_rate: si se indica esta opción (activada por defecto), se emplea el algoritmo (ADADELTA) como algoritmo de aprendizaje. Con él, se consigue que el learning rate se adapte automáticamente (reduciéndose) a medida que el entrenamiento avanza. Es recomendable probar en primer lugar esta opción, ya que consigue un balance entre velocidad de aprendizaje y resultados bastante bueno. Este algoritmo está controlado por dos parámetros:rho: controla el ratio en el que se va reduciendo el learning rate en cada iteración de aprendizaje.
epsilon: corrector para evitar divisiones entre cero.
Si se desactiva el adaptive_rate, entonces, el usuario debe especificar el valor de learning rate empleado por la red.
rate: Cuanto menor sea este valor, más estable es el modelo final pero más tarda en aprender (llegar a convergencia).nesterov_accelerated_gradient: activa el método de descenso de gradiente Nesterov Accelerated Gradient.
Regularización
input_dropout_ratio: porcentaje de dropout en la capa de entrada. Este tipo de regularización controla que ningún predictor influya en exceso en la red. Se recomiendan valores de 0.1 o 0.2.hidden_dropout_ratios: porcentaje de dropout en las capas de intermedias. Solo es aplicable si se selecciona como función de activación TanhWithDropout, RectifierWithDropout, o MaxoutWithDropout. Por defecto, el valor es 0.5.l1: penalización l1. Ver apartado de regularización GLM para más información.l2: penalización l2. Ver apartado de regularización GLM para más información.
El siguiente listado contiene el valor por defecto de cada uno de los argumentos de la función h2o.deeplearning().
nfolds= 0keep_cross_validation_predictions= FALSEkeep_cross_validation_fold_assignment= FALSEfold_assignment= c(“AUTO”, “Random”, “Modulo”, “Stratified”)fold_column= NULLignore_const_cols= TRUEscore_each_iteration= FALSEweights_column= NULLoffset_column= NULLbalance_classes= FALSEclass_sampling_factors= NULLmax_after_balance_size= 5max_hit_ratio_k= 0checkpoint= NULLpretrained_autoencoder= NULLoverwrite_with_best_model= TRUEuse_all_factor_levels= TRUEstandardize= TRUEactivation= c(“Tanh”, “TanhWithDropout”, “Rectifier”, “RectifierWithDropout”, “Maxout”, “MaxoutWithDropout”)hidden= c(200, 200)epochs= 10train_samples_per_iteration= -2target_ratio_comm_to_comp= 0.05seed= -1adaptive_rate= TRUErho= 0.99epsilon= 1e-08rate= 0.005rate_annealing= 1e-06rate_decay= 1momentum_start= 0momentum_ramp= 1e+06momentum_stable= 0nesterov_accelerated_gradient= TRUEinput_dropout_ratio= 0hidden_dropout_ratios= NULLl1= 0l2= 0max_w2= 3.4028235e+38initial_weight_distribution= c(“UniformAdaptive”, “Uniform”, “Normal”)initial_weight_scale= 1initial_weights= NULLinitial_biases= NULLloss= c(“Automatic”, “CrossEntropy”, “Quadratic”, “Huber”, “Absolute”, “Quantile”)distribution= c(“AUTO”, “bernoulli”, “multinomial”, “gaussian”, “poisson”, “gamma”, “tweedie”, “laplace”, “quantile”, “huber”)quantile_alpha= 0.5tweedie_power= 1.5huber_alpha= 0.9score_interval= 5score_training_samples= 10000score_validation_samples= 0score_duty_cycle= 0.1classification_stop= 0regression_stop= 1e-06stopping_rounds= 5stopping_metric= c(“AUTO”, “deviance”, “logloss”, “MSE”, “RMSE”, “MAE”, “RMSLE”, “AUC”, “lift_top_group”, “misclassification”, “mean_per_class_error”)stopping_tolerance= 0max_runtime_secs= 0score_validation_sampling= c(“Uniform”, “Stratified”)diagnostics= TRUEfast_mode= TRUEforce_load_balance= TRUEvariable_importances= TRUEreplicate_training_data= TRUEsingle_node_mode= FALSEshuffle_training_data= FALSEmissing_values_handling= c(“MeanImputation”, “Skip”)quiet_mode= FALSEautoencoder= FALSEsparse= FALSEcol_major= FALSEaverage_activation= 0sparsity_beta= 0max_categorical_features= 2147483647reproducible= FALSEexport_weights_and_biases= FALSEmini_batch_size= 1categorical_encoding= c(“AUTO”, “Enum”, “OneHotInternal”, “OneHotExplicit”, “Binary”, “Eigen”, “LabelEncoder”, “SortByResponse”, “EnumLimited”)elastic_averaging= FALSEelastic_averaging_moving_rate= 0.9elastic_averaging_regularization= 0.001verbose= FALSE
Importancia de los predictores
Para modelos de redes neuronales, una forma de cuantificar la importancia de los predictores es mediante el método Gedeon. Por defecto, esta opción está desactivada debido a que puede ser un proceso muy lento.
Modelo
# Hiperparámetros que se quieren comparar.
hiperparametros <- list(hidden = list(c(64), c(128), c(256), c(512), c(1024),
c(64,64), c(128,128), c(256,256),
c(512, 512)))
grid_dl <- h2o.grid(
# Algoritmo
algorithm = "deeplearning",
activation = "RectifierWithDropout",
epochs = 500,
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento.
training_frame = datos_train_h2o,
shuffle_training_data = FALSE,
# Datos de validación
validation_frame = datos_val_h2o,
# Preprocesado
standardize = TRUE,
missing_values_handling = "Skip",
# Detención temprana
stopping_rounds = 3,
stopping_metric = "AUC",
stopping_tolerance = 0.01,
# Hiperparámetros optimizados
hyper_params = hiperparametros,
# Regularización
l1 = 1e-5,
l2 = 1e-5,
# Tipo de búsqueda
search_criteria = list(strategy = "Cartesian"),
seed = 123,
grid_id = "grid_dl"
)# Se muestran los modelos ordenados de mayor a menor AUC.
resultados_grid <- h2o.getGrid(
grid_id = "grid_dl",
sort_by = "auc",
decreasing = TRUE
)
data.frame(resultados_grid@summary_table) %>% select(-model_ids)El modelo que consigue mayor AUC de validación es el que tiene una arquitectura de una única capa con 128 neuronas.
modelo_dl_final <- h2o.getModel(resultados_grid@model_ids[[1]])
plot(modelo_dl_final, timestep = "epochs", metric = "classification_error")
## [1] 0.9128137Stacked Ensembles
Introducción
A modo general, el término model ensembling hace referencia a la combinación de las predicciones de dos o más modelos distintos con el objetivo de mejorar las predicciones finales. Esta estrategia se basa en la asunción de que, distintos modelos entrenados independientemente, emplean distintos aspectos de los datos para realizar las predicciones, es decir, cada uno es capaz de identificar parte de la “verdad” pero no toda ella. Combinando la perspectiva de cada uno de ellos, se obtiene una descripción más detallada de la verdadera estructura subyacente en los datos. A modo de analogía, imagínese un grupo de estudiantes que se enfrentan a un examen multidisciplinal. Aunque todos obtengan aproximadamente la misma nota, cada uno de ellos habrá conseguido más puntos con las preguntas que tratan sobre las disciplinas en las que destacan. Si en lugar de hacer el examen de forma individual, lo hacen en grupo, cada uno podrá contribuir en los aspectos que más domina, y el resultado final será probablemente superior a cualquiera de los resultados individuales.
La clave para que el ensembling consiga mejorar los resultados es la diversidad de los modelos. Si todos los modelos combinados son similares entre ellos, no podrán compensarse unos a otros. Por esta razón, se tiene que intentar combinar modelos que sean lo mejor posible a nivel individual y lo más diferentes entre ellos.
Las formas más simples de combinar las predicciones de varios modelos es emplear la media para problemas de regresión y la moda para problemas de clasificación. Sin embargo existen otras aproximaciones más complejas capaces de conseguir mejores resultados:
Ponderar las agregaciones dando distinto peso a cada modelo, por ejemplo, en proporción al accuracy que han obtenido de forma individual.
Super Learner (stacked regression): emplear un nuevo modelo para que decida la mejor forma de combinar las predicciones de los otros modelos.
Algoritmo Super Learner
La implementación de Super Learner disponible en H2O sigue el siguiente algoritmo:
Definición del ensemble
Definir un listado con los algoritmos base (cada uno con los hiperparámetros pertinentes).
Definir el algoritmo de metalearning que defina como se entrena en modelo superior. Entre los metalearners disponibles están (GLM, GBM, Random Forest y Deep Learning).
Entrenamiento del ensemble
Entrenar cada uno de los algoritmos base con el conjunto de entrenamiento.
Realizar k-fold cross-validation con cada uno de los algoritmos base y almacenar las predicciones hechas en cada una de las k particiones.
Combinar las predicciones del paso 2 en una única matriz NxL (N = número de observaciones en el conjunto de entrenamiento, L = número de modelos base).
Entrenar el metalearning con la variable respuesta y la matriz NxL como predictores.
El Super learner final está formado por los modelos base y el modelo metalearning.
Predicción
Predecir la nueva observación con cada uno de los modelos base.
Emplear las predicciones de los modelos base como input del metalearner para obtener la predicción final.
Modelo
Se entrena un super learner empleando como algoritmos base un GLM y un GBM. Cada uno de los modelos base se entrena por separado empleando como hiperparámetros los valores óptimos identificados en los grid search de los apartados anteriores.
Es importante tener en cuenta que, en el entrenamiento de los modelos base, es necesario:
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).
# Entrenamiento del modelo GLM
modelo_glm <- h2o.glm(
y = var_respuesta,
x = predictores,
training_frame = datos_train_h2o,
family = "binomial",
link = "logit",
standardize = TRUE,
balance_classes = FALSE,
ignore_const_cols = TRUE,
# Se especifica que hacer con observaciones incompletas
missing_values_handling = "Skip",
# Se hace una búsqueda del hiperparámetro lamba
lambda_search = TRUE,
# Selección automática del solver adecuado
solver = "AUTO",
alpha = 0,
# Validación cruzada de 5 folds para estimar el error
# del modelo.
seed = 123,
nfolds = 5,
# Reparto estratificado de las observaciones en la creación
# de las particiones.
fold_assignment = "Stratified",
keep_cross_validation_predictions = TRUE,
model_id = "modelo_glm"
)
# Entrenamiento del modelo GBM
modelo_gbm <- h2o.gbm(
# Tipo de distribución (clasificación binaria)
distribution = "bernoulli",
# Variable respuesta y predictores
y = var_respuesta,
x = predictores,
# Datos de entrenamiento
training_frame = datos_train_h2o,
# Número de árboles
ntrees = 5000,
# Complejidad de los árboles
max_depth = 20,
min_rows = 10,
# Aprendizaje
learn_rate = 0.1,
sample_rate = 0.8,
# Detención temprana
score_tree_interval = 100,
stopping_rounds = 3,
stopping_metric = "AUC",
stopping_tolerance = 0.01,
# Validación cruzada de 5 folds para estimar el error
# del modelo.
seed = 123,
nfolds = 5,
# Reparto estratificado de las observaciones en la creación
# de las particiones.
fold_assignment = "Stratified",
keep_cross_validation_predictions = TRUE,
model_id = "modelo_gbm"
)# Entrenar el ensemble con los dos modelos anteriores.
modelo_ensemble <- h2o.stackedEnsemble(
y = var_respuesta,
x = predictores,
training_frame = datos_train_h2o,
base_models = list(modelo_glm, modelo_gbm),
model_id = "modelo_ensemble"
)## [1] 0.9071069Interpretación de modelos
En muchas aplicaciones reales, conocer la importancia que tiene cada predictor dentro de un modelo no es suficiente, además, se necesita entender de qué forma influye cada uno en el valor de las predicciones.
En los modelos lineales, gracias a los coeficientes de regresión, es fácil conocer en qué dirección y magnitud influyen variaciones de un predictor en el valor predicho. Sin embargo, en modelos más complejos como GBM, Random Forest o redes neuronales, no es tan inmediato. Una estrategia para estos últimos son los gráficos PDP e ICE.
Gráficos PDP
Los gráficos PDP muestran cómo varían en promedio las predicciones a medida que se modifica el valor de un predictor, manteniendo constantes el resto. Pueden obtenerse con la función h2o.partialPlot() pasándole como argumentos un modelo (object), las observaciones con las que se ha entrenado el modelo (data) y el nombre de los predictores para los que se quiere obtener el modelo (cols).
Véase el resultado de calcular las gráficas PDP del predictor age en el modelo GBM ajustado anteriormente.
h2o.partialPlot(
object = modelo_gbm_final,
data = datos_h2o,
cols = c("age"),
plot = TRUE,
plot_stddev = TRUE
)
El gráfico PDP muestra que, la relación entre la edad y la probabilidad de ganar más de 50000 dólares anuales no es constante ni monotónica. En el intervalo de los 20-40 años, a medida que aumenta la edad, mayor la probabilidad promedio de que una persona supere ese salario. En el intervalo 40-60 se estabiliza y a partir de los 60 decae.
También es posible representar la influencia de dos variables a la vez con gráficos PDP-2D. Para crear estos gŕaficos se necesita tener instalados los paquetes plot3Drgl y rgl.
p <- h2o.partialPlot(
object = modelo_glm_final,
data = datos_h2o,
cols = c("age", "education_number", "sex"),
col_pairs_2dpdp = list(c("age", "education_number"), c("age", "sex")),
plot = TRUE
)Gráficos ICE
Los gráficos ICE pueden considerarse una extensión de los gráficos de dependencia parcial Partial Dependence Plots (PDP). La diferencia entre ambos reside en que, los PDP, muestran, con una única curva, cómo varía en promedio la predicción de la variable respuesta a medida que se modifica uno de los predictores. Los ICE, muestran cómo varía la predicción para cada una de las observaciones (una curva distinta por cada observación). Para información más detallada consultar Gráficos ICE para interpretar modelos predictivos.
Con la función partial del paquete iml se pueden obtener, entre otras cosas, los gráficos ICE de modelos predictivos. En el caso particular de los modelos de H2O, como no se crean en el entorno de R, para obtener los gráficos ICE son necesarios dos pasos adicionales.
Definir una función
predict()que envíe los nuevos datos al cluster H2O en el formato adecuado, y que devuelva un vector con un valor numérico de predicción.Pasar los datos en formato
h2ode nuevo adata.frame.
Al tratarse de clasificación binaria, la función h2o.predict() devuelve por cada predicción 3 valores: la clase predicha y la probabilidad de cada una. Se necesita una función que recupere la probabilidad de >50K, que en este caso está en la tercera columna.
## predict <=50K >50K
## 1 <=50K 0.9999998 2.045886e-07
##
## [1 row x 3 columns]# Función predict especial para un modelo H2O de clasificación binaria.
predict_custom <- function(object, newdata){
h2o.no_progress()
as.vector(predict(object, newdata = as.h2o(newdata))[[3]])
}
# Ejemplo
predict_custom(object = modelo_gbm_final, newdata = datos_h2o[1,])## [1] 2.045886e-07Para poder obtener los gráficos ICE de varios predictores, se crea una función que engloba a la función partial() y que después se aplica a cada predictor con map().
library(pdp)
calcular_graficar_ice <- function(predictor, modelo_h2o, funcion_predict,
grid.resolution, datos_train_h2o, fraccion){
# Función para calcular y graficar curvas ICE de uno o varios predictores
datos_train <- as.data.frame(datos_h2o)
datos_train <- dplyr::sample_frac(tbl = datos_train, size = fraccion)
plot_ice <- partial(
pred.var = predictor,
object = modelo_h2o,
pred.fun = funcion_predict,
grid.resolution = grid.resolution,
train = datos_train,
ice = TRUE,
center = TRUE,
quantiles = FALSE,
plot = TRUE,
contour = FALSE,
rug = TRUE,
alpha = 0.1,
parallel = FALSE,
plot.engine = "ggplot2"
)
plot_ice <- plot_ice + geom_rug(data = datos_train,
aes(x = !!sym(predictor)),
inherit.aes = FALSE)
plot(plot_ice)
}
purrr::walk(
.x = c("age", "education_number"),
.f = calcular_graficar_ice,
modelo_h2o = modelo_gbm_final,
funcion_predict = predict_custom,
grid.resolution = 30,
datos_train_h2o = datos_train_h2o,
# Se emplea una fracción muy pequeña dato el elevado número de observaciones
fraccion = 0.3
)
Guardar un modelo H2O
Para guardar un modelo H2O ya entrenado se emplean la función h2o.saveModel().
# Se guarda el modelo en el directorio actual
h2o.saveModel(object = modelo_glm_final, path = getwd(), force = TRUE)# Se guarda el modelo en el directorio actual
h2o.saveModel(object = modelo_glm_final, path = getwd(), force = TRUE)
file.rename(file.path(getwd(), modelo_glm_final@model_id), "./modelo_glm")Para cargar un modelo previamente guardado se emplea la función h2o.loadModel().
Información sesión
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)Bibliografía
Generalized Linear Modeling with H2O Nykodym, T., Kraljevic, T., Hussami, N., Rao, A., and Wang, A. (Nov 2017).
Gradient Boosting Machine with H2O by Michal Malohlava & Arno Candel with assitance from Cliff Click, Hank Roark, & Viraj Parmar
Machine Learning with R and H2O by Mark Landry with assistance from Spencer Aiello, Eric Eckstrand, Anqi Fu, & Patrick Aboyoun Edited by: Angela Bartz
Dua, D. and Karra Taniskidou, E. (2017). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science
Practical Machine Learning with H2O by Darren Cook
Top 10 Deep Learning Tips and Tricks - Arno Candel
Linear Models with R by Julian J.Faraway
An Introduction to Statistical Learning: with Applications in R (Springer Texts in Statistics)
Extending the Linear Model with R: Generalized Linear, Mixed Effects and Nonparametric Regression Models by Julian J.Faraway
¿Cómo citar este documento?
Machine Learning con H2O y R por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/44_machine_learning_con_h2o_y_r
¿Te ha gustado el artículo? Tu ayuda es importante
Mantener un sitio web tiene unos costes elevados, tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊
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.