Redes neuronales con R


Más sobre ciencia de datos en cienciadedatos.net

Introducción


El ámbito de las redes neuronales y su hermano mayor, el deep learning, es complejo y amplio. Durante los últimos años, el interés y la aplicación de este tipo de modelos han experimentado tal expansión que se ha convertido en una disciplina por sí misma. Si bien entender bien sus fundamentos requiere de una cantidad notable de tiempo y práctica, esto no significa que se necesiten adquirir todos ellos para empezar a sacarles partido. En este documento se presenta una introducción, más intuitiva que rigurosa, sobre los modelos de redes neuronales y de cómo crearlos con R.

En el ecosistema de R, existen múltiples paquetes que permiten crear modelos basados en redes neuronales. Conviene diferenciar entre dos casos de uso ya que, dependiendo de estos, son más adecuadas unas librerías u otras:

  • Modelos de redes simples (perceptrón simple): estos modelos se caracterizan por tener arquitecturas relativamente sencillas por lo que los requerimientos computacionales no son elevados y no es necesario el uso el uso de GPUs. Dentro de este grupo destacan las implementaciones de nnet y neuralnet.

  • Deep learning: son modelos más complejos (perceptrón multicapa, redes convolucionales, redes recurrentes, LSTM…) cuyos requerimientos computacionales hacen necesario el uso de muy optimizado de CPUs o GPUs. Para este tipo de modelos se tiene que recurrir a frameworks especializados como H2O, Tensorflow-Keras o Torch.

En este documento se muestran ejemplos de modelos relativamente sencillos (perceptrón simple y perceptrón multicapa) haciendo uso del paquete H2O.

Redes neuronales


Estructura de la red


Las redes neuronales son modelos creados al ordenar operaciones matemáticas siguiendo una determinada estructura. La forma más común de representar la estructura de una red neuronal es mediante el uso de capas (layers), formadas a su vez por neuronas (unidades, units o neurons). Cada neurona, realiza una operación sencilla y está conectada a las neuronas de la capa anterior y de la capa siguiente mediante pesos, cuya función es regular la información que se propaga de una neurona a otra.

Representación de una red neuronal perceptrón simple. Fuente: Computer Age Statistical Inference 2016.



La primera capa de la red neuronal (color verde) se conoce como capa de entrada o input layer y recibe los datos en bruto, es decir, el valor de los predictores. La capa intermedia (color azul), conocida como capa oculta o hidden layer, recibe los valores de la capa de entrada, ponderados por los pesos (flechas grises). La última capa, llamada output layer, combina los valores que salen de la capa intermedia para generar la predicción.

Para facilitar la comprensión de la estructura de las redes, es útil representar una red equivalente a un modelo de regresión lineal.

\[y = w_1 x_1 + ... + w_d x_d + b\]

Representación de una red neuronal equivalente a un modelo lineal con 4 predictores. Fuente: COMS W4995 Applied Machine Learning.



Cada neurona de la capa de entrada representa el valor de uno de los predictores. Las flechas representan los coeficientes de regresión, que en términos de redes se llaman pesos, y la neurona de salida representa el valor predicho. Para que esta representación equivalga a la ecuación de un modelo lineal, faltan dos cosas:

  • El bias del modelo.

  • Las operaciones de multiplicación y suma que combinan el valor de los predictores con los pesos del modelo.

Cada neurona de la capa intermedia tiene un valor de bias, pero suele omitirse en las representaciones gráficas. En cuanto a las operaciones matemáticas, es el elemento clave que ocurre dentro de las neuronas y conviene verlo con detalle.

La neurona (unidad)


La neurona es la unidad funcional de los modelos de redes. Dentro de cada neurona ocurren simplemente dos operaciones: la suma ponderada de sus entradas y la aplicación de una función de activación.

En la primera parte, se multiplica cada valor de entrada \(x_i\) por su peso asociado \(w_i\) y se suman junto con el bias. Este es el valor neto de entrada a la neurona. A continuación, este valor se pasa por una función, conocida como función de activación, que transforma el valor neto de entrada en un valor de salida.

Si bien el valor que llega a la neurona, siempre es una combinación lineal, gracias a la función de activación, se pueden generar salidas muy diversas. Es en la función de activación donde reside el potencial de los modelos de redes para aprender relaciones no lineales.

Representación de una neurona. Fuente: Deep Learning A Practitioner’s Approach by Josh Patterson and Adam Gibson.



La anterior ha sido una explicación intuitiva del funcionamiento de una neurona. Véase ahora una definición más matemática.

El valor neto de entrada a una neurona es la suma de los valores que le llegan, ponderados por el peso de las conexiones, más el bias.

\[entrada = \sum^n_{i=1} x_i w_i + b\]

En lugar de utilizar el sumatorio, esta operación suele representarse como el producto matricial, donde \(\textbf{X}\) representa el vector de los valores de entrada y \(\textbf{W}\) el vector de pesos.

\[entrada = \textbf{X}\textbf{W} + b\]

A este valor se le aplica una función de activación (\(g\)) que lo transforma en lo que se conoce como valor de activación (\(a\)), que es lo que finalmente sale de la neurona.

\[a = g(entrada) = g(\textbf{X}\textbf{W} + b)\]

Para la capa de entrada, donde únicamente se quiere incorporar el valor de los predictores, la función de activación es la unidad, es decir, sale lo mismo que entra. En la capa de salida, la función de activación utilizada suele ser la identidad para problemas de regresión y softmax para clasificación.

Función de activación


Las funciones de activación controlan en gran medida qué información se propaga desde una capa a la siguiente (forward propagation). Estas funciones convierten el valor neto de entrada a la neuronal (combinación de los input, pesos y bias) en un nuevo valor. Gracias a combinar funciones de activación no lineales con múltiples capas (ver más adelante), los modelos de redes son capaces de aprender relaciones no lineales.

La gran mayoría de funciones de activación convierten el valor de entrada neto de la neurona en un valor dentro del rango (0, 1) o (-1, 1). Cuando el valor de activación de una neurona (salida de su función de activación) es cero, se dice que la neurona está inactiva, ya que no pasa ningún tipo de información a las siguientes neuronas. A continuación, se describen las funciones de activación más empleadas.

Rectified linear unit (ReLU)

La función de activación ReLu aplica una transformación no lineal muy simple, activa la neurona solo si el input está por encima de cero. Mientras el valor de entrada está por debajo de cero, el valor de salida es cero, pero cuando es superior, el valor de salida aumenta de forma lineal con el de entrada.

\[\operatorname{ReLU}(x) = \max(x, 0)\]

De esta forma, la función de activación retiene únicamente los valores positivos y descarta los negativos dándoles una activación de cero.

Representación función activación ReLU.



ReLU es con diferencia la función de activación más empleada por sus buenos resultados en aplicaciones diversas. La razón de esto reside en el comportamiento de su derivada (gradiente), que es cero o constante. Gracias a esto se evita el problema de vanishing gradients que limita la capacidad de aprendizaje de los modelos de redes.

Sigmoide

La función sigmoide transforma valores en el rango de (-inf, +inf) a valores en el rango (0, 1).

\[\operatorname{sigmoid}(x) = \frac{1}{1 + \exp(-x)}\]

Representación función activación sigmoide.



Aunque la función de activación sigmoide se utilizó mucho en los inicios de los modelos de redes, en la actualidad, suele preferirse la función ReLU.

Un caso en el que la función de activación sigmoide sigue siendo la función utilizada por defecto es en las neuronas de la capa de salida de los modelos de clasificación binaria, ya que su salida puede interpretarse como probabilidad.

Tangente hiperbólica (Tanh)

La función de activación Tanh, se comporta de forma similar a la función sigmoide, pero su salida está acotada en el rango (-1, 1).

\[\operatorname{tanh}(x) = \frac{1 - \exp(-2x)}{1 + \exp(-2x)}\]



Representación función activación Tanh.



Sin las funciones de activación, las redes neuronales solo pueden aprender relaciones lineales.

Función de coste (loss function)


La función de coste (\(l\)), también llamada función de pérdida, loss function o cost function, es la encargada de cuantificar la distancia entre el valor real y el valor predicho por la red, en otras palabras, mide cuánto se equivoca la red al realizar predicciones. En la mayoría de casos, la función de coste devuelve valores positivos. Cuanto más próximo a cero es el valor de coste, mejor son las predicciones de la red (menor error), siendo cero cuando las predicciones se corresponden exactamente con el valor real.

La función de coste puede calcularse para una única observación o para un conjunto de datos (normalmente promediando el valor de todas las observaciones). Es el segundo caso el que se utiliza para dirigir el entrenamiento de los modelos.

Dependiendo del tipo de problema, regresión o clasificación, es necesario utilizar una función de coste u otra. En problemas de regresión, las más utilizadas son el error cuadrático medio y el error absoluto medio. En problemas de clasificación suele emplearse la función log loss, también llamada logistic loss o cross-entropy loss.

Error cuadrático medio

El error cuadrático medio (mean squared error, MSE) es con diferencia la función de coste más utilizada en problemas de regresión. Para una determinada observación \(i\), el error cuadrático se calcula como la diferencia al cuadrado entre el valor predicho \(\hat{y}\) y el valor real \(y\).

\[l^{(i)}(\mathbf{w}, b) = \left(\hat{y}^{(i)} - y^{(i)}\right)^2\]

Las funciones de coste suelen escribirse con la notación \(l(\mathbf{w}, b)\) para hacer referencia a que su valor depende de los pesos y bias del modelo, ya que son estos los que determinan el valor de las predicciones \(y^{(i)}\).

Con frecuencia, esta función de coste se encuentra multiplicada por \(\frac{1}{2}\), esto es simplemente por conveniencia matemática para simplificar el cálculo de su derivada.

\[l^{(i)}(\mathbf{w}, b) = \frac{1}{2} \left(\hat{y}^{(i)} - y^{(i)}\right)^2\]

Para cuantificar el error que comete el modelo todo un conjunto de datos, por ejemplo los de entrenamiento, se promedia el error de todas las \(N\) observaciones.

\[L(\mathbf{w}, b) =\frac{1}{n}\sum_{i=1}^n l^{(i)}(\mathbf{w}, b) = \frac{1}{n}\sum_{i=1}^n \left(\hat{y}^{(i)} - y^{(i)}\right)^2\]

Cuando un modelo se entrena utilizando el error cuadrático medio como función de coste, está aprendiendo a predecir la media de la variable respuesta.

Error medio absoluto

El error medio absoluto (mean absolute error, MAE) consiste en promediar el error absoluto de las predicciones.

\[L(\mathbf{w}, b) =\frac{1}{n}\sum_{i=1}^n |\hat{y}^{(i)} - y^{(i)}|\]

El error medio absoluto es más robusto frente a outliers que el error cuadrático medio. Esto significa que, el entrenamiento del modelo, se ve menos influenciado por datos anómalos que pueda haber en el conjunto de entrenamiento. Cuando un modelo se entrena utilizando el error absoluto medio como función de coste, está aprendiendo a predecir la mediana de la variable respuesta.

Log loss, logistic loss o cross-entropy loss

En problemas de clasificación, la capa de salida utiliza como función de activación la función softmax. Gracias a esta función, la red devuelve una serie de valores que pueden interpretarse como la probabilidad de que la observación predicha pertenezca a cada una de las posibles clases.

Cuando la clasificación es de tipo binaria, donde la variable respuesta es 1 o 0, y \(p = \operatorname{Pr}(y = 1)\), la función de coste log-loss se define como:

\[L_{\log}(y, p) = -\log \operatorname{Pr}(y|p) = -(y \log (p) + (1 - y) \log (1 - p))\]

Para problemas de clasificación con más de dos clases, esta fórmula se generaliza a:

\[L_{\log}(Y, P) = -\log \operatorname{Pr}(Y|P) = - \frac{1}{N} \sum_{i=0}^{N-1} \sum_{k=0}^{K-1} y_{i,k} \log p_{i,k}\]

En ambos casos, minimizar esta la función equivale a que la probabilidad predicha para la clase correcta tienda a 1 y a 0 en las demás clases.

Dado que esta función se ha utilizado en campos diversos, se le conoce por nombres distintos: Log loss, logistic loss o cross-entropy loss, pero todos hacen referencia a lo mismo. Puede encontrarse una explicación más detallada de esta función de coste aquí.

Múltiples capas


El modelo de red neuronal con una única capa (single-layer perceptron), aunque supuso un gran avance en el campo del machine learning, solo es capaz de aprender patrones sencillos. Para superar esta limitación, los investigadores descubrieron que, combinando múltiples capas ocultas, la red puede aprender relaciones mucho más complejas entre los predictores y la variable respuesta. A esta estructura se le conoce como perceptrón multicapa o multilayer perceptron (MLP), y puede considerarse como el primer modelo de deep learning.

La estructura de un perceptrón multicapa consta de varias capas ocultas. Cada neurona está conectada a todas las neuronas de la capa anterior y a las de la capa posterior. Aunque no es estrictamente necesario, todas las neuronas que forman parte de una misma capa suelen emplear la misma función de activación.

Combinando múltiples capas ocultas y funciones de activación no lineales, los modelos de redes pueden aprender prácticamente cualquier patrón. De hecho, está demostrado que, con suficientes neuronas, un MLP es un aproximador universal para cualquier función.

Representación de una red neuronal perceptrón multicapa. Fuente: Computer Age Statistical Inference 2016.



Entrenamiento


El proceso de entrenamiento de una red neuronal consiste en ajustar el valor de los pesos y bias de tal forma que, las predicciones que se generen, tengan el menor error posible. Gracias a esto, el modelo es capaz de identificar qué predictores tienen mayor influencia y de qué forma están relacionados entre ellos y con la variable respuesta.

La idea intuitiva de cómo entrenar una red neuronal es la siguiente:

  1. Iniciar la red con valores aleatorios de los pesos y bias.

  2. Para cada observación de entrenamiento, calcular el error que comete la red al hacer su predicción. Promediar los errores de todas las observaciones.

  3. Identificar la responsabilidad que ha tenido cada peso y bias en el error de las predicciones.

  4. Modificar ligeramente los pesos y bias de la red (de forma proporcional a su responsabilidad en el error) en la dirección correcta para que se reduzca el error.

  5. Repetir los pasos 2, 3, 4 y 5 hasta que la red sea suficientemente buena.

Si bien la idea parece sencilla, alcanzar una forma de implementarla ha requerido la combinación de múltiples métodos matemáticos, en concreto, el algoritmo de retropropagación (backpropagation) y la optimización por descenso de gradiente (gradient descent).

Backpropagation


Backpropagation es el algoritmo que permite cuantificar la influencia que tiene cada peso y bias en las predicciones de la red. Para conseguirlo, hace uso de la regla de la cadena (chain rule) para calcular el gradiente, que no es más que el vector formado por las derivadas parciales de una función.

En el caso de las redes, la derivada parcial del error respecto a un parámetro (peso o bias) mide cuánta “responsabilidad” ha tenido ese parámetro en el error cometido. Gracias a esto, se puede identificar qué pesos de la red hay que modificar para mejorarla. El siguiente paso necesario, es determinar cuánto y cómo modificarlos (optimización).

Descenso de gradiente


Descenso de gradiente (gradient descent) es un algoritmo de optimización que permite minimizar una función haciendo actualizaciones de sus parámetros en la dirección del valor negativo de su gradiente. Aplicado a las redes neuronales, el descenso de gradiente permite ir actualizando los pesos y bias del modelo para reducir su error.

Dado que, calcular el error del modelo para todas las observaciones de entrenamiento en cada iteración puede ser computacionalmente muy costoso, suele utilizarse una alternativa al método de descenso de gradiente llamada gradiente estocástico (stochastic gradient descent, SGD). Este método consiste en dividir el conjunto de entrenamiento en lotes (minibatch o batch) y actualizar los parámetros de la red con cada uno. De esta forma, en lugar de esperar a evaluar todas las observaciones para actualizar los parámetros, se pueden ir actualizando de forma progresiva. Una ronda completa sobre todos los batch se llama época. El número de épocas con las que se entrena una red equivale al número de veces que la red ve cada ejemplo de entrenamiento.

Preprocesado


A la hora de entrenar modelos basados en redes neuronales, es necesario aplicar a los datos, al menos, dos tipos de transformaciones.

Binarización (one hot encoding) de las variables categóricas

La binarización (one-hot-encoding) consiste en crear nuevas variables dummy con cada uno de los niveles de las variables cualitativas. Por ejemplo, una variable llamada color que contenga los niveles rojo, verde y azul, se convertirá en tres nuevas variables (color_rojo, color_verde, color_azul), todas con el valor 0 excepto la que coincide con la observación, que toma el valor 1.

Estandarización y escalado de variables numéricas

Cuando los predictores son numéricos, la escala en la que se miden, así como la magnitud de su varianza pueden influir en gran medida en el modelo. Si no se igualan de alguna forma los predictores, aquellos que se midan en una escala mayor o que tengan más varianza dominarán el modelo aunque no sean los que más relación tienen con la variable respuesta. Existen principalmente 2 estrategias para evitarlo:

  • Centrado: consiste en restarle a cada valor la media del predictor al que pertenece. Si los datos están almacenados en un dataframe, el centrado se consigue restándole a cada valor la media de la columna en la que se encuentra. Como resultado de esta transformación, todos los predictores pasan a tener una media de cero, es decir, los valores se centran en torno al origen.

  • Normalización (estandarización): consiste en transformar los datos de forma que todos los predictores estén aproximadamente en la misma escala.

    • Normalización Z-score: dividir cada predictor entre su desviación típica después de haber sido centrado, de esta forma, los datos pasan a tener un rango de [\(+\infty\), \(-\infty\)].

    • Estandarización max-min: transformar los datos de forma que estén dentro del rango [0, 1].

Nunca se deben estandarizar las variables después de ser binarizadas.

Hiperparámetros


La gran “flexibilidad” que tienen las redes neuronales es un arma de doble filo. Por un lado, son capaces de generar modelos que aprenden relaciones muy complejas, sin embargo, sufren fácilmente el problema de sobreajuste (overfitting) lo que los incapacita al tratar de predecir nuevas observaciones. La forma de minimizar este problema y conseguir modelos útiles pasa por configurar de forma adecuada sus hiperparámetros. Son muchos los hiperparámetros de un modelo basado en redes y su nomenclatura varía de unas implementaciones a otras, sin embargo, los de mayor impacto siempre están presentes:

  • Número y tamaño de capas

  • Learning rate (tasa de aprendizaje)

  • Algoritmo de optimización

  • Regularización

Número y tamaño de capas


La arquitectura de una red, el número de capas y el número de neuronas que forman parte de cada capa, determinan en gran medida la complejidad del modelo y con ello su potencial capacidad de aprendizaje.

La capa de entrada y salida son sencillas de establecer. La capa de entrada tiene tantas neuronas como predictores y la capa de salida tiene una neurona en problemas de regresión y tantas como clases en problemas de clasificación. En la mayoría de implementaciones, estos valores se establecen automáticamente en función del conjunto de entrenamiento. El usuario suele especificar únicamente el número de capas intermedias (ocultas) y el tamaño de las mismas.

Cuantas más neuronas y capas, mayor la complejidad de las relaciones que puede aprender el modelo. Sin embargo, dado que cada neurona está conectada por pesos al resto de neuronas de las capas adyacentes, el número de parámetros a aprender aumenta y con ello el tiempo de entrenamiento.

Learning rate


El learning rate o ratio de aprendizaje establece cómo de rápido pueden cambiar los parámetros de un modelo a medida que se optimiza (aprende). Este hiperparámetro es uno de los más complicados de establecer, ya que depende mucho de los datos e interacciona con el resto de hiperparámetros. Si el learning rate es muy grande, el proceso de optimización puede ir saltando de una región a otra sin que el modelo sea capaz de aprender. Si por el contrario, el learning rate es muy pequeño, el proceso de entrenamiento puede tardar demasiado y no llegar a completarse. Algunas de las recomendaciones heurísticas basadas en prueba y error son:

  • Utilizar un learning rate lo más pequeño posible siempre y cuando el tiempo de entrenamiento no supere las limitaciones temporales disponibles.

  • No utilizar un valor constante de learning rate durante todo el proceso de entrenamiento. Por lo general, utilizar valores mayores al inicio y pequeños al final.

Algoritmo de optimización


El descenso de gradiente y el descenso de gradiente estocástico fueron de los primeros métodos de optimización utilizados para entrenar modelos de redes neuronales. Ambos utilizan directamente el gradiente para dirigir la optimización. Pronto se vio que esto genera problemas a medida que las redes aumentan de tamaño (neuronas y capas). En muchas regiones del espacio de búsqueda, el gradiente es muy próximo a cero, lo que hace que la optimización quede estancada. Para evitar este problema, se han desarrollado modificaciones del descenso de gradiente capaces de adaptar el learning rate en función del gradiente y subgradiente. De esta forma, el proceso de aprendizaje se ralentiza o acelera dependiendo de las características de la región del espacio de búsqueda en el que se encuentren. Aunque existen multitud de adaptaciones, suele recomendarse:

  • Para conjuntos de datos pequeños: l-bfgs

  • Para conjuntos de datos grandes: adam o rmsprop

La elección del algoritmo de optimización puede tener un impacto notable en el aprendizaje de los modelos, sobre todo en deep learning. Puede encontrarse una excelente descripción más detallada en el libro gratuito Dive into Deep Learning.

Regularización


Los métodos de regularización tienen el objetivo de reducir el sobreajuste (overfitting) de los modelos. Un modelo con sobreajuste memoriza los datos de entrenamiento pero es incapaz de predecir correctamente nuevas observaciones.

Los modelos de redes neuronales pueden considerarse como modelos sobre parametrizados, por lo tanto, las estrategias de regularización son fundamentales. De entre las muchas que existen, destacan la regularización L1 y L2 (weight decay) y el dropout.

Regularización L1 y L2

El objetivo de la regularización L1 y L2, esta última también conocida como weight decay, es evitar que los pesos tomen valores excesivamente elevados. De esta forma, se evita que unas pocas neuronas dominen el comportamiento de la red y se fuerza a que las características poco informativas (ruido) tengan pesos próximos o iguales a cero.

Dropout

Este proceso consiste en desactivar aleatoriamente una serie de neuronas durante el proceso de entrenamiento. En concreto, durante cada iteración del entrenamiento, se ponen a cero los pesos de una fracción aleatoria de neuronas por capa. El método de dropout, descrito por Srivastava et al. en 2014, se ha convertido en un estándar para entrenar redes neuronales. El porcentaje de neuronas que suele desactivarse por capa (dropout rate) suele ser un valor entre 0.2 y 0.5.

Modelos de redes neuronales con H2O


¿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.

Modelos


Para crear modelos basados en redes neuronales con H2O, se utiliza la función h2o.deeplearning(). Son muchos los argumentos que controlan el comportamiento de este tipo de modelos. Afortunadamente, los responsables de su implementación han establecido valores por defecto que suelen funcionar adecuadamente en muchos escenarios. A continuación, se muestran 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 con hidden = c(100) o hidden = 100, y una red con dos capas ocultas de 10 neuronas cada una con hidden = 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 la opción de excluirlos Skip o imputarlos con la media MeanImputation.

  • shuffle_training_data: mezclar aleatoriamente las observaciones antes de entrenar el modelo.

  • categorical_encoding: codificado de las variables categóricas.

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 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 overfitting. Normalmente, se puede tener una idea aproximada del valor correcto en función de la complejidad de la red 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: learning 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 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.

  • l2: penalización l2.



En este link puede encontrarse el valor por defecto de cada uno de los argumentos de la función h2o.deeplearning().



Ejemplo clasificación


En este primer ejemplo se muestra cómo, dependiendo de la arquitectura (capas ocultas y tamaño de las mismas), un modelo basado en redes neuronales puede aprender funciones no lineales con gran facilidad. Para evitar problemas de overfitting, es importante identificar la combinación de hiperparámetros que consigue un equilibrio adecuado en el aprendizaje. Se trata de un ejemplo muy sencillo, cuyo objetivo es que el lector se familiarice con la flexibilidad que tienen este tipo de modelos y en cómo realizar una búsqueda de hiperparámetros mediante grid search y validación cruzada.

Packages


Los paquetes utilizadas en este ejemplo son:

# Tratamiento de datos y gráficos
# ==============================================================================
library(tidyverse)
library(ggthemes)
library(ggpubr)

# Modelado
# ==============================================================================
# sudo apt install default-jre
library(h2o)



Iniciación del cluster


Una vez que la librería H2O ha sido cargada, hay que iniciar el cluster. 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"
)
##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         4 hours 34 minutes 
##     H2O cluster timezone:       Europe/Madrid 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.32.1.3 
##     H2O cluster version age:    1 month and 21 days  
##     H2O cluster name:           H2O_started_from_R_ximo_tcy429 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   5.93 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.1.0 (2021-05-18)
# 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()



Datos


Se simulan observaciones en dos dimensiones, pertenecientes a tres grupos, cuya separación no es lineal.

# Datos simulados
# ==============================================================================
datos = read_csv('https://raw.githubusercontent.com/JoaquinAmatRodrigo/Estadistica-machine-learning-python/master/data/blobs.csv')

datos <- datos %>% mutate(y = as.factor(y))

ggplot(data = datos, aes(x = x_1, y = x_2, fill = y)) + 
  geom_point(shape = 21, size = 2) +
  theme_fivethirtyeight() + 
  theme(
    legend.position = "none",
    text = element_blank(),
    axis.ticks =  element_blank()
  )

datos_h2o   <- as.h2o(datos)
particiones <- h2o.splitFrame(data = datos_h2o, ratios = c(0.6, 0.2), seed = 123)
datos_train <- h2o.assign(data = particiones[[1]], key = "datos_train")
datos_validation  <- h2o.assign(data = particiones[[2]], key = "datos_validacion")
datos_test  <- h2o.assign(data = particiones[[3]], key = "datos_test")



Arquitectura de la red


Se procede a crear 4 modelos en orden creciente de complejidad (número de neuronas y capas), para comprobar cómo la arquitectura de la red afecta a su capacidad de aprendizaje. 

# Modelos
# ==============================================================================
modelo_1 <- h2o.deeplearning(
                  x               = c("x_1", "x_2"),
                  y               = "y",
                  distribution    = "multinomial",
                  training_frame  = datos_train,
                  standardize     = TRUE,
                  activation      = "Rectifier",
                  adaptive_rate   = FALSE,
                  hidden          = 1,
                  stopping_rounds = 0,
                  epochs          = 1000,
                  seed            = 123,
                  model_id        = "modelo_1"
           )

modelo_2 <- h2o.deeplearning(
                x               = c("x_1", "x_2"),
                y               = "y",
                distribution    = "multinomial",
                training_frame  = datos_train,
                standardize     = TRUE,
                activation      = "Rectifier",
                adaptive_rate   = FALSE,
                hidden          = 10,
                stopping_rounds = 0,
                epochs          = 1000,
                seed            = 123,
                model_id        = "modelo_2"
            )

modelo_3 <- h2o.deeplearning(
                  x               = c("x_1", "x_2"),
                  y               = "y",
                  distribution    = "multinomial",
                  training_frame  = datos_train,
                  standardize     = TRUE,
                  activation      = "Rectifier",
                  adaptive_rate   = FALSE,
                  hidden          = c(10, 10),
                  stopping_rounds = 0,
                  epochs          = 1000,
                  seed            = 123,
                  model_id        = "modelo_3"
            )

modelo_4 <- h2o.deeplearning(
                  x               = c("x_1", "x_2"),
                  y               = "y",
                  distribution    = "multinomial",
                  training_frame  = datos_train,
                  standardize     = TRUE,
                  activation      = "Rectifier",
                  adaptive_rate   = FALSE,
                  hidden          = c(50, 50, 50),
                  stopping_rounds = 0,
                  epochs          = 1000,
                  seed            = 123,
                  model_id        = "modelo_4"
            )
# Predicciones de cada modelo
# ==============================================================================

grid_predicciones <- expand.grid(
                        x_1 = seq(from = min(datos$x_1), to = max(datos$x_1), length = 75),
                        x_2 = seq(from = min(datos$x_2), to = max(datos$x_2), length = 75)
                     )

grid_predicciones_h2o <- as.h2o(grid_predicciones)


predicciones_1 <- h2o.predict(
                    object  = modelo_1,
                    newdata = grid_predicciones_h2o
                  )

predicciones_2 <- h2o.predict(
                    object  = modelo_2,
                    newdata = grid_predicciones_h2o
                  )

predicciones_3 <- h2o.predict(
                    object  = modelo_3,
                    newdata = grid_predicciones_h2o
                  )

predicciones_4 <- h2o.predict(
                    object  = modelo_4,
                    newdata = grid_predicciones_h2o
                  )


grid_predicciones$modelo_1 <- as.vector(predicciones_1$predict)
grid_predicciones$modelo_2 <- as.vector(predicciones_2$predict)
grid_predicciones$modelo_3 <- as.vector(predicciones_3$predict)
grid_predicciones$modelo_4 <- as.vector(predicciones_4$predict)
# Gráfico de predicciones
# ==============================================================================
p1 <- ggplot(data = grid_predicciones, aes(x = x_1, y = x_2, color = modelo_1)) + 
      geom_point(size = 0.5) +
      theme_fivethirtyeight() +  
      labs(title = "Arquitectura: (5)") +
      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 = modelo_2)) + 
      geom_point(size = 0.5) +
      labs(title = "Arquitectura: (10)") +
      theme_fivethirtyeight() + 
      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 = modelo_3)) + 
      geom_point(size = 0.5) +
      labs(title = "Arquitectura: (20, 20)") +
      theme_fivethirtyeight() + 
      theme(legend.position = "none",
            plot.title = element_text(size=11),
            axis.text  = element_blank(),
            axis.title = element_blank(),
            axis.ticks =  element_blank())

p4 <- ggplot(data = grid_predicciones, aes(x = x_1, y = x_2, color = modelo_4)) + 
      geom_point(size = 0.5) +
      labs(title = "Arquitectura: (50, 50, 50)") +
      theme_fivethirtyeight() + 
      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, p4, nrow = 2, ncol = 2)

Puede observarse cómo, a medida que aumenta la complejidad de la red (más neuronas y más capas), las fronteras de decisión se adaptan más y más a los datos de entrenamiento.

Optimización de hiperparámetros


En este apartado, se muestra cómo afectan al aprendizaje algunos de los hiperparámetros más influyentes. Para ello, se combinan los métodos de grid search, random search y early stopping.

Detención temprana del entrenamiento (Early Stopping)


Una de las características de los modelos de deep learning es que, con el número suficiente de épocas, 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 épocas y, para ello, suele tener que entrenar el modelo con cientos o miles de épocas 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 realizando iteraciones innecesarias.

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 Deeplearning 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. Tres argumentos controlan la estrategia de parada:

  • stopping_metric: métrica empleada para cuantificar cuánto mejora el modelo.

  • stopping_tolerance: porcentaje mínimo de mejora entre dos mediciones consecutivas 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 el stopping_tolerance para 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=2, stopping_tolerance=0.01 y stopping_metric="misclassification".

  • Detener el entrenamiento de un modelo cuando el logloss no se reduce nada durante 3 mediciones consecutivas: stopping_rounds=3, stopping_tolerance=0 y stopping_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.001 y stopping_metric="AUC".

También es posible detener el entrenamiento del modelo cuando se supera un tiempo máximo especificado con el argumento max_runtimesecs > 0.

modelo <- h2o.deeplearning(
                      x = c("x_1", "x_2"),
                      y = "y",
                      distribution    = "multinomial",
                      training_frame  = datos_train,
                      validation_frame = datos_validation,
                      standardize     = TRUE,
                      activation      = "Rectifier",
                      adaptive_rate   = FALSE,
                      rate            = 0.1,
                      hidden          = 10,
                      stopping_rounds = 5,
                      stopping_metric = "misclassification",
                      stopping_tolerance = 0.001,
                      score_validation_samples = 1000,
                      epochs          = 100000000, # Múchas epocas
                      seed            = 123
          )
modelo@model$scoring_history 
plot(modelo)



Caso de uso: predicción de precios vivienda


El set de datos SaratogaHouses del paquete mosaicData contiene información sobre la precio de 1728 viviendas situadas en Saratoga County, New York, USA en el año 2006. Además del precio, incluye 15 variables adicionales:

  • price: precio de la vivienda.
  • lotSize: metros cuadrados de la vivienda.
  • age: antigüedad de la vivienda.
  • landValue: valor del terreno.
  • livingArea: metros cuadrados habitables.
  • pctCollege: porcentaje del vecindario con título universitario.
  • bedrooms: número de dormitorios.
  • firplaces: número de chimeneas.
  • bathrooms: número de cuartos de baño (el valor 0.5 hace referencia a cuartos de baño sin ducha).
  • rooms: número de habitaciones.
  • heating: tipo de calefacción.
  • fuel: tipo de alimentación de la calefacción (gas, electricidad o diesel).
  • sewer: tipo de desagüe.
  • waterfront: si la vivienda tiene vistas al lago.
  • newConstruction: si la vivienda es de nueva construcción.
  • centralAir: si la vivienda tiene aire acondicionado.

El objetivo es obtener un modelo de red neuronal capaz de predecir el precio de la vivienda.

Paquetes


# Tratamiento de datos y gráficos
# ==============================================================================
library(tidymodels)
library(tidyverse)
library(skimr)
library(DataExplorer)
library(ggpubr)
library(mosaicData)

# Modelado
# ==============================================================================
# sudo apt install default-jre
library(h2o)



Datos


data("SaratogaHouses", package = "mosaicData")
datos <- SaratogaHouses

# Se renombran las columnas para que sean más descriptivas
colnames(datos) <- c("precio", "metros_totales", "antiguedad", "precio_terreno",
                     "metros_habitables", "universitarios",
                     "dormitorios", "chimenea", "banyos", "habitaciones",
                     "calefaccion", "consumo_calefacion", "desague",
                     "vistas_lago","nueva_construccion", "aire_acondicionado")



Análisis exploratorio


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

Tabla resumen

# Tabla resumen
# ==============================================================================
skim(datos)
Data summary
Name datos
Number of rows 1728
Number of columns 16
_______________________
Column type frequency:
factor 6
numeric 10
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
calefaccion 0 1 FALSE 3 hot: 1121, ele: 305, hot: 302
consumo_calefacion 0 1 FALSE 3 gas: 1197, ele: 315, oil: 216
desague 0 1 FALSE 3 pub: 1213, sep: 503, non: 12
vistas_lago 0 1 FALSE 2 No: 1713, Yes: 15
nueva_construccion 0 1 FALSE 2 No: 1647, Yes: 81
aire_acondicionado 0 1 FALSE 2 No: 1093, Yes: 635

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
precio 0 1 211966.71 98441.39 5000 1.45e+05 189900.00 259000.00 775000.0 ▅▇▂▁▁
metros_totales 0 1 0.50 0.70 0 1.70e-01 0.37 0.54 12.2 ▇▁▁▁▁
antiguedad 0 1 27.92 29.21 0 1.30e+01 19.00 34.00 225.0 ▇▁▁▁▁
precio_terreno 0 1 34557.19 35021.17 200 1.51e+04 25000.00 40200.00 412600.0 ▇▁▁▁▁
metros_habitables 0 1 1754.98 619.94 616 1.30e+03 1634.50 2137.75 5228.0 ▇▇▂▁▁
universitarios 0 1 55.57 10.33 20 5.20e+01 57.00 64.00 82.0 ▁▃▆▇▁
dormitorios 0 1 3.15 0.82 1 3.00e+00 3.00 4.00 7.0 ▃▇▅▁▁
chimenea 0 1 0.60 0.56 0 0.00e+00 1.00 1.00 4.0 ▆▇▁▁▁
banyos 0 1 1.90 0.66 0 1.50e+00 2.00 2.50 4.5 ▁▇▇▁▁
habitaciones 0 1 7.04 2.32 2 5.00e+00 7.00 8.25 12.0 ▃▇▇▅▂

Todas las columnas tienen el tipo adecuado.

# Número de datos ausentes por variable
# ==============================================================================
datos %>% map_dbl(.f = function(x){sum(is.na(x))})
##             precio     metros_totales         antiguedad     precio_terreno 
##                  0                  0                  0                  0 
##  metros_habitables     universitarios        dormitorios           chimenea 
##                  0                  0                  0                  0 
##             banyos       habitaciones        calefaccion consumo_calefacion 
##                  0                  0                  0                  0 
##            desague        vistas_lago nueva_construccion aire_acondicionado 
##                  0                  0                  0                  0
plot_missing(
  data    = datos, 
  title   = "Porcentaje de valores ausentes",
  ggtheme = theme_bw(),
  theme_config = list(legend.position = "none")
)

Todas las columnas están completas, no hay valores ausentes.

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

# Tabla de estadísticos principales 
summary(datos$precio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  145000  189900  211967  259000  775000

Los modelos de redes neuronales son de tipo no paramétrico, no asumen ningún tipo de distribución de la variable respuesta, por lo tanto, no es necesario que esta siga ninguna distribución concreta (normal, gamma…). Aun así, siempre es recomendable hacer un estudio mínimo, ya que, a fin de cuentas, es lo que interesa predecir. En este caso, la variable precio tiene una distribución asimétrica con una cola positiva debido a que, unas pocas viviendas, tienen un precio muy superior a la media.

# Gráfico de distribución para cada variable numérica
# ==============================================================================
plot_density(
  data    = datos %>% select(-precio),
  ncol    = 3,
  title   = "Distribución variables continuas",
  ggtheme = theme_bw(),
  theme_config = list(
                  plot.title = element_text(size = 14, face = "bold"),
                  strip.text = element_text(colour = "black", size = 12, face = 2)
                 )
  )

La variable chimenea, aunque es de tipo numérico, apenas toma unos pocos valores y la gran mayoría de observaciones pertenecen a solo dos de ellos. En casos como este, suele ser conveniente tratar la variable como cualitativa.

# Valores observados de chimenea
# ==============================================================================
table(datos$chimenea)
## 
##   0   1   2   3   4 
## 740 942  42   2   2
# Se convierte la variable chimenea a factor.
datos <- datos %>%
         mutate(chimenea = as.factor(chimenea))
# Gráfico para cada variable cualitativa
# ==============================================================================
plot_bar(
  datos,
  ncol    = 3,
  title   = "Número de observaciones por grupo",
  ggtheme = theme_bw(),
  theme_config = list(
                   plot.title = element_text(size = 14, face = "bold"),
                   strip.text = element_text(colour = "black", size = 8, face = 2),
                   legend.position = "none"
                  )
)

Si alguno de los niveles de una variable cualitativa tiene muy pocas observaciones en comparación a los otros niveles, puede ocurrir que, durante la validación cruzada o bootstrapping, algunas particiones no contengan ninguna observación de dicha clase (varianza cero), lo que puede dar lugar a errores. Para este caso, hay que tener precaución con la variable chimenea. Se unifican los niveles de 2, 3 y 4 en un nuevo nivel llamado “2_mas”.

table(datos$chimenea)
## 
##   0   1   2   3   4 
## 740 942  42   2   2
datos <- datos %>%
         mutate(
           chimenea = recode_factor(
                        chimenea,
                        `2` = "2_mas",
                        `3` = "2_mas",
                        `4` = "2_mas"
                      )
         )

table(datos$chimenea)
## 
## 2_mas     0     1 
##    46   740   942



División train y test


Con el objetivo de poder estimar el error que comete el modelo al predecir nuevas observaciones, se dividen los datos en dos grupos, uno de entrenamiento y otro de test (80%, 20%).

# Reparto de datos en train y test
# ==============================================================================
set.seed(123)
split_inicial <- initial_split(
                    data   = datos,
                    prop   = 0.8,
                    strata = precio
                 )

datos_train <- training(split_inicial)
datos_test  <- testing(split_inicial)

Tras realizar el reparto, se verifica que los dos grupos son similares.

summary(datos_train$precio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5000  145000  189050  212423  259000  775000
summary(datos_test$precio)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   20000  145000  189950  210157  257750  620000



Preprocesado


Los modelos de redes neuronales requieren como mínimo de dos tipos de preprocesado: binarización (one hot encoding) de las variables categóricas y estandarización de las variables continuas.

# Se almacenan en un objeto `recipe` todos los pasos de preprocesado y, finalmente,
# se aplican a los datos.
transformer <- recipe(
                  formula = precio ~ .,
                  data =  datos_train
               ) %>%
               step_naomit(all_predictors()) %>%
               step_nzv(all_predictors()) %>%
               step_center(all_numeric(), -all_outcomes()) %>%
               step_scale(all_numeric(), -all_outcomes()) %>%
               step_dummy(all_nominal(), -all_outcomes())

transformer
## Data Recipe
## 
## Inputs:
## 
##       role #variables
##    outcome          1
##  predictor         15
## 
## Operations:
## 
## Removing rows with NA values in all_predictors()
## Sparse, unbalanced variable filter on all_predictors()
## Centering for all_numeric(), -all_outcomes()
## Scaling for all_numeric(), -all_outcomes()
## Dummy variables from all_nominal(), -all_outcomes()

Una vez que se ha definido el objeto recipe, con la función prep() se aprenden las transformaciones con los datos de entrenamiento y se aplican a los dos conjuntos con bake().

# Se entrena el objeto recipe
transformer_fit <- prep(transformer)

# Se aplican las transformaciones al conjunto de entrenamiento y de test
datos_train_prep <- bake(transformer_fit, new_data = datos_train)
datos_test_prep  <- bake(transformer_fit, new_data = datos_test)

glimpse(datos_train_prep)
## Rows: 1,380
## Columns: 19
## $ metros_totales              <dbl> -0.43894231, 0.26607619, 0.48189818, 0.611…
## $ antiguedad                  <dbl> 3.6226807, 0.1147142, 0.2866733, -0.917040…
## $ precio_terreno              <dbl> -0.7674967, -0.5829431, -0.3570715, -0.354…
## $ metros_habitables           <dbl> 0.31203819, -0.96981699, -0.19293506, -0.2…
## $ universitarios              <dbl> -0.4558859, -3.2928560, -0.4558859, -0.455…
## $ dormitorios                 <dbl> 1.0465989, 1.0465989, -0.1681248, -0.16812…
## $ banyos                      <dbl> -1.3602353, -1.3602353, -0.6045490, 0.1511…
## $ habitaciones                <dbl> 0.425229547, 0.425229547, 0.425229547, -0.…
## $ precio                      <int> 109000, 120000, 90000, 120000, 85860, 1270…
## $ chimenea_X0                 <dbl> 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, …
## $ chimenea_X1                 <dbl> 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, …
## $ calefaccion_hot.water.steam <dbl> 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, …
## $ calefaccion_electric        <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, …
## $ consumo_calefacion_electric <dbl> 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, …
## $ consumo_calefacion_oil      <dbl> 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 1, …
## $ desague_public.commercial   <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, …
## $ desague_none                <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
## $ nueva_construccion_No       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
## $ aire_acondicionado_No       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, …

Tras el preprocesado de los datos, se han generado un total de 19 variables (18 predictores y la variable respuesta).

Modelado


# Inicialización del cluster
# ==============================================================================
h2o.init(
  nthreads = -1,
  max_mem_size = "4g"
)
##  Connection successful!
## 
## R is connected to the H2O cluster: 
##     H2O cluster uptime:         4 hours 42 minutes 
##     H2O cluster timezone:       Europe/Madrid 
##     H2O data parsing timezone:  UTC 
##     H2O cluster version:        3.32.1.3 
##     H2O cluster version age:    1 month and 22 days  
##     H2O cluster name:           H2O_started_from_R_ximo_tcy429 
##     H2O cluster total nodes:    1 
##     H2O cluster total memory:   5.91 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.1.0 (2021-05-18)
# Se eliminan los datos del cluster por si ya había sido iniciado.
h2o.removeAll()
h2o.no_progress()

Se transfieren los datos al cluster de H2O..

datos_train  <- as.h2o(datos_train_prep, key = "datos_train")
datos_test   <- as.h2o(datos_test_prep, key = "datos_test")
# Espacio de búsqueda de cada hiperparámetro
# ==============================================================================
hiperparametros <- list(
                      epochs = c(50, 100, 500),
                      hidden = list(5, 10, 25, 50, c(10, 10))
                    )


# Búsqueda por validación cruzada
# ==============================================================================
variable_respuesta <- 'precio'
predictores <- setdiff(colnames(datos_train), variable_respuesta)

grid <- h2o.grid(
              algorithm    = "deeplearning",
              activation   = "Rectifier",
              x            = predictores,
              y            = variable_respuesta,
              training_frame  = datos_train,
              nfolds       = 3, #validacion cruzada
              standardize  = FALSE,
              hyper_params = hiperparametros,
              search_criteria = list(strategy = "Cartesian"),
              seed         = 123,
              grid_id      = "grid"
          )
# Resultados del grid
# ==============================================================================
resultados_grid <- h2o.getGrid(
                     sort_by = 'rmse',
                     grid_id = "grid",
                     decreasing = FALSE
                   )
data.frame(resultados_grid@summary_table)
# Mejor modelo encontrado
# ==============================================================================
modelo_final <- h2o.getModel(resultados_grid@model_ids[[1]])



Error de test


Aunque mediante los métodos de validación se consiguen buenas estimaciones del error que tiene un modelo al predecir nuevas observaciones, la mejor forma de evaluar un modelo final es prediciendo un conjunto test, es decir, un conjunto de observaciones que se ha mantenido al margen del proceso de entrenamiento y optimización.

predicciones <- h2o.predict(
                  object  = modelo_final,
                  newdata = datos_test
                )

predicciones <- predicciones %>%
                as_tibble() %>%
                mutate(valor_real = as.vector(datos_test$precio))

predicciones %>% head(5)
rmse(predicciones, truth = valor_real, estimate = predict, na_rm = TRUE)



Conclusión


La combinación de hiperparámetros con la que se obtienen mejores resultados acorde a las métricas de validación cruzada es:

modelo_final@allparameters
## $model_id
## [1] "grid_model_12"
## 
## $training_frame
## [1] "datos_train_prep_sid_a36b_13"
## 
## $nfolds
## [1] 3
## 
## $keep_cross_validation_models
## [1] TRUE
## 
## $keep_cross_validation_predictions
## [1] FALSE
## 
## $keep_cross_validation_fold_assignment
## [1] FALSE
## 
## $fold_assignment
## [1] "Random"
## 
## $ignore_const_cols
## [1] TRUE
## 
## $score_each_iteration
## [1] FALSE
## 
## $balance_classes
## [1] FALSE
## 
## $max_after_balance_size
## [1] 5
## 
## $max_confusion_matrix_size
## [1] 20
## 
## $overwrite_with_best_model
## [1] FALSE
## 
## $use_all_factor_levels
## [1] TRUE
## 
## $standardize
## [1] FALSE
## 
## $activation
## [1] "Rectifier"
## 
## $hidden
## [1] 50
## 
## $epochs
## [1] 504.5542
## 
## $train_samples_per_iteration
## [1] -2
## 
## $target_ratio_comm_to_comp
## [1] 0.05
## 
## $seed
## [1] 123
## 
## $adaptive_rate
## [1] TRUE
## 
## $rho
## [1] 0.99
## 
## $epsilon
## [1] 1e-08
## 
## $rate
## [1] 0.005
## 
## $rate_annealing
## [1] 1e-06
## 
## $rate_decay
## [1] 1
## 
## $momentum_start
## [1] 0
## 
## $momentum_ramp
## [1] 1e+06
## 
## $momentum_stable
## [1] 0
## 
## $nesterov_accelerated_gradient
## [1] TRUE
## 
## $input_dropout_ratio
## [1] 0
## 
## $l1
## [1] 0
## 
## $l2
## [1] 0
## 
## $max_w2
## [1] 3.402823e+38
## 
## $initial_weight_distribution
## [1] "UniformAdaptive"
## 
## $initial_weight_scale
## [1] 1
## 
## $loss
## [1] "Automatic"
## 
## $distribution
## [1] "gaussian"
## 
## $quantile_alpha
## [1] 0.5
## 
## $tweedie_power
## [1] 1.5
## 
## $huber_alpha
## [1] 0.9
## 
## $score_interval
## [1] 5
## 
## $score_training_samples
## [1] 10000
## 
## $score_validation_samples
## [1] 0
## 
## $score_duty_cycle
## [1] 0.1
## 
## $classification_stop
## [1] 0
## 
## $regression_stop
## [1] 1e-06
## 
## $stopping_rounds
## [1] 0
## 
## $stopping_tolerance
## [1] 0
## 
## $max_runtime_secs
## [1] 0
## 
## $score_validation_sampling
## [1] "Uniform"
## 
## $diagnostics
## [1] TRUE
## 
## $fast_mode
## [1] TRUE
## 
## $force_load_balance
## [1] TRUE
## 
## $variable_importances
## [1] TRUE
## 
## $replicate_training_data
## [1] TRUE
## 
## $single_node_mode
## [1] FALSE
## 
## $shuffle_training_data
## [1] FALSE
## 
## $missing_values_handling
## [1] "MeanImputation"
## 
## $quiet_mode
## [1] FALSE
## 
## $autoencoder
## [1] FALSE
## 
## $sparse
## [1] FALSE
## 
## $col_major
## [1] FALSE
## 
## $average_activation
## [1] 0
## 
## $sparsity_beta
## [1] 0
## 
## $max_categorical_features
## [1] 2147483647
## 
## $reproducible
## [1] FALSE
## 
## $export_weights_and_biases
## [1] FALSE
## 
## $mini_batch_size
## [1] 1
## 
## $categorical_encoding
## [1] "OneHotInternal"
## 
## $elastic_averaging
## [1] FALSE
## 
## $elastic_averaging_moving_rate
## [1] 0.9
## 
## $elastic_averaging_regularization
## [1] 0.001
## 
## $auc_type
## [1] "AUTO"
## 
## $x
##  [1] "metros_totales"              "antiguedad"                 
##  [3] "precio_terreno"              "metros_habitables"          
##  [5] "universitarios"              "dormitorios"                
##  [7] "banyos"                      "habitaciones"               
##  [9] "chimenea_X0"                 "chimenea_X1"                
## [11] "calefaccion_hot.water.steam" "calefaccion_electric"       
## [13] "consumo_calefacion_electric" "consumo_calefacion_oil"     
## [15] "desague_public.commercial"   "desague_none"               
## [17] "nueva_construccion_No"       "aire_acondicionado_No"      
## 
## $y
## [1] "precio"
# Se apaga el cluster H2O
#h2o.shutdown(prompt = FALSE)



Información sesión


sesion_info <- devtools::session_info()
dplyr::select(
  tibble::as_tibble(sesion_info$packages),
  c(package, loadedversion, source)
)



Bibliografía


The Elements of Statistical Learning (2nd edition)

Hastie, Tibshirani and Friedman (2009). Springer-Verlag.

Applied Predictive Modeling by Max Kuhn and Kjell Johnson

Deep Learning by Josh Patterson, Adam Gibson



¿Cómo citar este documento?

Redes neuronales con R por Joaquín Amat Rodrigo, disponible con licencia CC BY-NC-SA 4.0 en https://www.cienciadedatos.net/documentos/68-redes-neuronales-r.html


¿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! 😊


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