Capítulo 8 Árboles de decisión

Un árbol de decisiones es un algoritmo del aprendizaje supervisado que se puede utilizar tanto para problemas de clasificación como de regresión. Es un clasificador estructurado en árbol, donde los nodos internos representan las características de un conjunto de datos, las ramas representan las reglas de decisión y cada nodo hoja representa el resultado. La idea básica de los árboles es buscar puntos de cortes en las variables de entrada para hacer predicciones, ir dividiendo la muestra, y encontrar cortes sucesivos para refinar las predicciones.

En un árbol de decisión, hay dos tipos nodos, el nodo de decisión o nodos internos (Decision Node) y el nodo hoja o nodo terminal (Leaf node). Los nodos de decisión se utilizan para tomar cualquier decisión y tienen múltiples ramas, mientras que los nodos hoja son el resultado de esas decisiones y no contienen más ramas.

  • Regresión: En el caso de la regresión de árboles de decisión, en los nodos finales se calcula el promedio de la variable de respuesta. El promedio será la estimación del modelo.

  • Clasificación: Por otro lado, en los árboles de clasificación se calcula la proporción de elementos de cada categoría en los nodos finales. De esta manera se calcula la probabilidad de pertenencia a la categoría.

8.1 Ajuste del modelo

En un árbol de decisión, para predecir la clase del conjunto de datos, el algoritmo comienza desde el nodo raíz del árbol. Este algoritmo compara los valores de la variable raíz con la variable de registro y, según la comparación, sigue una rama y salta al siguiente nodo.

Para el siguiente nodo, el algoritmo vuelve a comparar el valor de la siguiente variable con los otros sub-nodos y avanza. Continúa el proceso hasta que se llega a un nodo hoja. El proceso completo se puede comprender mejor con los siguientes pasos:

  1. Comenzamos el árbol con el nodo raíz (llamado S), que contiene el conjunto de entrenamiento completo.

  2. Encuentre la mejor variable en el conjunto de datos usando Attribute Selective Measure (ASM).

  3. Divida la S en subconjuntos que contengan valores posibles para la mejor variable.

  4. Genere el nodo del árbol de decisión, que contiene la mejor variable.

  5. Cree de forma recursiva nuevos árboles de decisión utilizando los subconjuntos del conjunto de datos creado en el paso 3. Continúe este proceso hasta que se alcance una etapa en la que no pueda particionar más los nodos y este nodo final sera un nodo hoja.

  6. Para clasificación nos quedaremos la moda de la variable respuesta del nodo hoja y para regresión usaremos la media de la variable respuesta.

8.1.1 Attribute Selective Measure (ASM)

Al implementar un árbol de decisión, surge el problema principal de cómo seleccionar la mejor variable para el nodo raíz y para los sub-nodos. Para resolver este problemas existe una técnica que se llama medida de selección de atributos o ASM. Mediante esta medición, podemos seleccionar fácilmente la mejor variable para los nodos del árbol. Una de las técnicas más populares para ASM es:

  • Índice de Gini

La medida del grado de probabilidad de que una variable en particular se clasifique incorrectamente cuando se elige al azar se llama índice de Gini o impureza de Gini. Los datos se distribuyen por igual según el índice de Gini.

\[Gini = \sum_{i=1}^{n}\hat{p_i}(1-\hat{p}_i)\]

Con \(p_i\) como la probabilidad de que un objeto se clasifique en una clase particular.

Esta métrica puede analizarse como una métrica de impureza. Cuando todos o la mayoría de elementos dentro de un nodo pertenecen a una misma clase, el índice de Gini toma valores cercanos a cero.

Cuando se utiliza el índice de Gini como criterio seleccionar la variable para el nodo raíz, seleccionaremos la variable con el índice de Gini menor.

8.2 Regularización de árboles

Para asegurarse de que no exista sobre-ajuste en el modelo, es importante considerar algunas regularizaciones a los hiper-parámetros implementados. Posteriormente, se determinará cuál de las posibles combinaciones produce mejores resultados.

8.2.1 Nivel de profundidad de árbol

Podríamos preguntarnos cuándo dejar de crecer un árbol. Pueden existir problemas que tengan un gran conjunto de variables y esto da como resultado una gran cantidad de divisiones, lo que a su vez genera un árbol de decisión muy grande. Estos árboles son complejos y pueden provocar un sobre-ajuste. Entonces, necesitamos saber cuándo parar.

  1. Una forma de hacer esto, es establecer un número mínimo de entradas de entrenamiento para dividir un nodo (min_n).

  2. Otra forma, es establecer la profundidad máxima del modelo. La profundidad máxima se refiere a la longitud del camino más largo desde el nodo raíz hasta un nodo hoja (max_depth).

8.3 Aprendizaje conjunto

El aprendizaje conjunto da crédito a la idea de la “sabiduría de las multitudes”, lo que sugiere que la toma de decisiones de un grupo más grande de individuos (modelos) suele ser mejor que la de un individuo.

El aprendizaje en conjunto es un grupo (o conjunto) de individuos o modelos, que trabajan colectivamente para lograr una mejor predicción final. Un solo modelo, también conocido como aprendiz básico puede no funcionar bien individualmente debido a una gran variación o un alto sesgo, sin embargo, cuando se agregan individuos débiles, pueden formar un individuo fuerte, ya que su combinación reduce el sesgo o la varianza, lo que produce un mejor rendimiento del modelo.

Los métodos de conjunto se ilustran con frecuencia utilizando árboles de decisión, ya que este algoritmo puede ser propenso a sobre ajustar (alta varianza y bajo sesgo) y también puede prestarse a desajuste (baja varianza y alto sesgo) cuando es muy pequeño, como un árbol de decisión con un nivel.

Nota: Cuando un algoritmo se adapta o no se adapta a su conjunto de entrenamiento, no se puede generalizar bien a nuevos conjuntos de datos, por lo que se utilizan métodos de conjunto para contrarrestar este comportamiento y permitir la generalización del modelo a nuevos conjuntos de datos.

8.4 Bagging

Primero tenemos que definir qué es la Agregación de Bootstrap o Bagging. Este es un algoritmo de aprendizaje automático diseñado para mejorar la estabilidad y precisión de algoritmos de ML usados en clasificación estadística y regresión. Además reduce la varianza y ayuda a evitar el sobre-ajuste. Aunque es usualmente aplicado a métodos de árboles de decisión, puede ser usado con cualquier tipo de método. Bagging es un caso especial del promediado de modelos.

Los métodos de bagging son métodos donde los algoritmos simples son usados en paralelo. El principal objetivo de los métodos en paralelo es el de aprovecharse de la independencia que hay entre los algoritmos simples, ya que el error se puede reducir bastante al promediar las salidas de los modelos simples. Es como si, queriendo resolver un problema entre varias personas independientes unas de otras, damos por bueno lo que eligiese la mayoría de las personas.

Para obtener la agregación de las salidas de cada modelo simple e independiente, bagging puede usar la votación para los métodos de clasificación y el promedio para los métodos de regresión.

El bagging o agregación bootstrap, es un método de aprendizaje por conjuntos que se usa comúnmente para reducir la varianza dentro de un conjunto de datos ruidoso.

8.5 Random Forest

Un bosque aleatorio es un algoritmo de aprendizaje automático supervisado que se construye a partir de algoritmos de árbol de decisión. Este algoritmo se aplica en diversas industrias, como la banca y el comercio electrónico, para predecir el comportamiento y los resultados.

En esta clase se dará una descripción general del algoritmo de bosque aleatorio, cómo funciona y las características del algoritmo.

También se señalan las ventajas y desventajas de este algoritmo.

8.5.1 ¿Qué es?

Un bosque aleatorio es una técnica de aprendizaje automático que se utiliza para resolver problemas de regresión y clasificación. Utiliza el aprendizaje por conjuntos, que es una técnica que combina muchos clasificadores para proporcionar soluciones a problemas complejos.

Este algoritmo consta de muchos árboles de decisión. El “bosque” generado se entrena mediante agregación de bootstrap (bagging), el cual es es un meta-algoritmo de conjunto que mejora la precisión de los algoritmos de aprendizaje automático.

El algoritmo establece el resultado en función de las predicciones de los árboles de decisión. Predice tomando el promedio o la media de la salida de varios árboles. El aumento del número de árboles aumenta la precisión del resultado.

Un bosque aleatorio erradica las limitaciones de un algoritmo de árbol de decisión. Reduce el sobre-ajuste de conjuntos de datos y aumenta la precisión. Genera predicciones sin requerir muchas configuraciones.

8.5.2 Características de los bosques aleatorios

  • Es más preciso que el algoritmo árbol de decisiones.

  • Proporciona una forma eficaz de gestionar los datos faltantes.

  • Puede producir una predicción razonable sin ajuste de hiperparámetros.

  • Resuelve el problema del sobre-ajuste en los árboles de decisión.

  • En cada árbol forestal aleatorio, se selecciona aleatoriamente un subconjunto de características en el punto de división del nodo.

8.5.3 Aplicar árboles de decisión en un bosque aleatorio

La principal diferencia entre el algoritmo de árbol de decisión y el algoritmo de bosque aleatorio es que el establecimiento de nodos raíz y la desagregación de nodos se realiza de forma aleatoria en este último. El bosque aleatorio emplea el método de bagging para generar la predicción requerida.

El método bagging implica el uso de diferentes muestras de datos (datos de entrenamiento) en lugar de una sola muestra. Los árboles de decisión producen diferentes resultados, dependiendo de los datos de entrenamiento alimentados al algoritmo de bosque aleatorio.

Nuestro primer ejemplo todavía se puede utilizar para explicar cómo funcionan los bosques aleatorios. Supongamos que solo tenemos cuatro árboles de decisión. En este caso, los datos de entrenamiento que comprenden las observaciones y características de estudio se dividirán en cuatro nodos raíz. Supongamos que queremos modelar si un cliente compra o no compra un teléfono.

Los nodos raíz podrían representar cuatro características que podrían influir en la elección de un cliente (precio, almacenamiento interno, cámara y RAM). El bosque aleatorio dividirá los nodos seleccionando características al azar. La predicción final se seleccionará en función del resultado de los cuatro árboles.

El resultado elegido por la mayoría de los árboles de decisión será la elección final.

Si tres árboles predicen la compra y un árbol predice que no comprará, entonces la predicción final será la compra. En este caso, se prevé que el cliente comprará.

El siguiente diagrama muestra un clasificador de bosque aleatorio simple.

8.5.4 Ventajas y desventjas de bosques aleatorios

Ventajas

  • Puede realizar tareas de regresión y clasificación.

  • Un bosque aleatorio produce buenas predicciones que se pueden entender fácilmente.

  • Puede manejar grandes conjuntos de datos de manera eficiente.

  • Proporciona un mayor nivel de precisión en la predicción de resultados sobre el algoritmo del árbol de decisión.

Desventajas

  • Cuando se usa un bosque aleatorio, se requieren bastantes recursos para el cálculo.

  • Consume más tiempo en comparación con un algoritmo de árbol de decisiones.

  • No producen buenos resultados cuando los datos son muy escasos. En este caso, el subconjunto de características y la muestra de arranque producirán un espacio invariante. Esto conducirá a divisiones improductivas, que afectarán el resultado.

8.6 Implementación en Python

Usaremos los pipelines antes implementados para ajustar tanto el modelo de regresión como el de clasificación. Una ventaja que se explorará al usar KFCV es el conjunto de hiperparámetros para elegir el mejor modelo posible.

Los pasos a seguir, son los siguientes:

  1. Carga de librerías
  2. Carga y separación inicial de datos ( test, train ).
  3. Pre-procesamiento e ingeniería de variables.
  4. Selección de tipo de modelo con hiperparámetros iniciales.
  5. Cálculo de métricas de desempeño.
  6. Creación de grid search y métricas de desempeño.
  7. Entrenamiento de modelos con hiperparámetros definidos.
  8. Análisis de métricas de error e hiperparámetros (Vuelve al paso 3, si es necesario).
  9. Selección de modelo a usar.
  10. Ajuste de modelo final con todos los datos.
  11. Validar poder predictivo con datos de prueba.

8.6.1 Regresión

Paso 0: Carga de librerías

from mlxtend.feature_selection import ColumnSelector
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.ensemble import RandomForestRegressor
from sklearn.pipeline import Pipeline
from sklearn.metrics import mean_absolute_error, mean_absolute_percentage_error
from sklearn.metrics import mean_squared_error, r2_score, make_scorer
from sklearn.model_selection import train_test_split, KFold, cross_val_score, cross_validate
from sklearn.model_selection import GridSearchCV
from sklearn.utils import shuffle
from sklearn import set_config

from plydata.one_table_verbs import pull
from plydata.tidy import pivot_longer
from mizani.formatters import comma_format, dollar_format
from plotnine import *
from siuba import *

import pandas as pd
import numpy as np
import pickle

Paso 1: Carga y separación inicial de datos ( test, train )

Comenzamos por cargar los datos completos e identificar la variable de respuesta para separarla de las explicativas

#### CARGA DE DATOS ####
ames = pd.read_csv("data/ames.csv")
ames = (ames >> 
 mutate(
  Second_Flr_SF = _.Second_Flr_SF.astype(float),
  Full_Bath = _.Full_Bath.astype(int),
  Half_Bath = _.Half_Bath.astype(int))
 )

ames_y = ames >> pull("Sale_Price")    # ames[["Sale_Price"]]
ames_x = select(ames, -_.Sale_Price)   # ames.drop('Sale_Price', axis=1)

#### DIVISIÓN DE DATOS ####
ames_x_train, ames_x_test, ames_y_train, ames_y_test = train_test_split(
 ames_x, ames_y, 
 train_size = 0.80, 
 random_state = 195
 )

Contando con datos de entrenamiento, procedemos a realizar el feature engineering para extraer las mejores características que permitirán realizar las estimaciones en el modelo.

Paso 2: Pre-procesamiento e ingeniería de variables

En este paso se pone a prueba la imaginación y conocimiento de transformaciones estadísticas que permitan ofrecer un mejor ajuste a la relación entre datos dependientes y de respuesta.

EL primer pipeline es para facilitar las transformaciones sobre las columnas sugeridas inicialmente:

## SELECCIÓN DE VARIABLES

# Seleccionamos las variales numéricas de interés
num_cols = ["Full_Bath", "Half_Bath", "Second_Flr_SF"]

# Seleccionamos las variables categóricas de interés
cat_cols = ["Overall_Cond"]

# Juntamos todas las variables de interés
columnas_seleccionadas = num_cols + cat_cols

pipe = ColumnSelector(columnas_seleccionadas)
ames_x_train_selected = pipe.fit_transform(ames_x_train)

ames_train_selected = (pd.DataFrame(
  ames_x_train_selected, 
  columns = columnas_seleccionadas) >> 
 mutate(
  Second_Flr_SF = _.Second_Flr_SF.astype(float),
  Full_Bath = _.Full_Bath.astype(int),
  Half_Bath = _.Half_Bath.astype(int))
 )

ames_train_selected.info()
## <class 'pandas.core.frame.DataFrame'>
## RangeIndex: 2344 entries, 0 to 2343
## Data columns (total 4 columns):
##  #   Column         Non-Null Count  Dtype  
## ---  ------         --------------  -----  
##  0   Full_Bath      2344 non-null   int64  
##  1   Half_Bath      2344 non-null   int64  
##  2   Second_Flr_SF  2344 non-null   float64
##  3   Overall_Cond   2344 non-null   object 
## dtypes: float64(1), int64(2), object(1)
## memory usage: 73.4+ KB

Y una vez que se han sugerido algunas columnas, se procede a realizar el pipeline de transformación:

## TRANSFORMACIÓN DE COLUMNAS

# ColumnTransformer para aplicar transformaciones
preprocessor = ColumnTransformer(
    transformers = [
        ('scaler', StandardScaler(), num_cols),
        ('onehotencoding', OneHotEncoder(drop='first', sparse_output=False), cat_cols)
    ],
    verbose_feature_names_out = False,
    remainder = 'passthrough'  # Mantener las columnas restantes sin cambios
)

transformed_data = preprocessor.fit_transform(ames_train_selected)
new_column_names = preprocessor.get_feature_names_out()

transformed_df = pd.DataFrame(
  transformed_data,
  columns=new_column_names
  )

transformed_df
##       Full_Bath  Half_Bath  Second_Flr_SF  Overall_Cond_Average  \
## 0          0.78       1.24           0.68                  1.00   
## 1          0.78      -0.75          -0.78                  1.00   
## 2          0.78      -0.75           0.73                  1.00   
## 3         -1.01      -0.75          -0.78                  1.00   
## 4         -1.01       1.24          -0.78                  0.00   
## ...         ...        ...            ...                   ...   
## 2339      -1.01       1.24           0.50                  0.00   
## 2340       0.78      -0.75          -0.78                  1.00   
## 2341      -1.01      -0.75          -0.78                  1.00   
## 2342       0.78       1.24           1.38                  1.00   
## 2343       0.78      -0.75           1.36                  1.00   
## 
##       Overall_Cond_Below_Average  Overall_Cond_Excellent  Overall_Cond_Fair  \
## 0                           0.00                    0.00               0.00   
## 1                           0.00                    0.00               0.00   
## 2                           0.00                    0.00               0.00   
## 3                           0.00                    0.00               0.00   
## 4                           0.00                    0.00               0.00   
## ...                          ...                     ...                ...   
## 2339                        1.00                    0.00               0.00   
## 2340                        0.00                    0.00               0.00   
## 2341                        0.00                    0.00               0.00   
## 2342                        0.00                    0.00               0.00   
## 2343                        0.00                    0.00               0.00   
## 
##       Overall_Cond_Good  Overall_Cond_Poor  Overall_Cond_Very_Good  \
## 0                  0.00               0.00                    0.00   
## 1                  0.00               0.00                    0.00   
## 2                  0.00               0.00                    0.00   
## 3                  0.00               0.00                    0.00   
## 4                  1.00               0.00                    0.00   
## ...                 ...                ...                     ...   
## 2339               0.00               0.00                    0.00   
## 2340               0.00               0.00                    0.00   
## 2341               0.00               0.00                    0.00   
## 2342               0.00               0.00                    0.00   
## 2343               0.00               0.00                    0.00   
## 
##       Overall_Cond_Very_Poor  
## 0                       0.00  
## 1                       0.00  
## 2                       0.00  
## 3                       0.00  
## 4                       0.00  
## ...                      ...  
## 2339                    0.00  
## 2340                    0.00  
## 2341                    0.00  
## 2342                    0.00  
## 2343                    0.00  
## 
## [2344 rows x 11 columns]
transformed_df.info()
## <class 'pandas.core.frame.DataFrame'>
## RangeIndex: 2344 entries, 0 to 2343
## Data columns (total 11 columns):
##  #   Column                      Non-Null Count  Dtype  
## ---  ------                      --------------  -----  
##  0   Full_Bath                   2344 non-null   float64
##  1   Half_Bath                   2344 non-null   float64
##  2   Second_Flr_SF               2344 non-null   float64
##  3   Overall_Cond_Average        2344 non-null   float64
##  4   Overall_Cond_Below_Average  2344 non-null   float64
##  5   Overall_Cond_Excellent      2344 non-null   float64
##  6   Overall_Cond_Fair           2344 non-null   float64
##  7   Overall_Cond_Good           2344 non-null   float64
##  8   Overall_Cond_Poor           2344 non-null   float64
##  9   Overall_Cond_Very_Good      2344 non-null   float64
##  10  Overall_Cond_Very_Poor      2344 non-null   float64
## dtypes: float64(11)
## memory usage: 201.6 KB

Recordemos que la función ColumnTransformes() solo son los pasos a seguir, necesitamos usar el método fit() que nos devuelve una receta actualizada con las estimaciones y la función transform() que nos devuelve la matriz transformada. Estos pasos pueden resumirse con el método: fit_transform()

Una vez que la receta de transformación de datos está lista, procedemos a implementar el pipeline del modelo de interés.

Paso 3: Selección de tipo de modelo con hiperparámetros iniciales

Ahora se procede a unir en un mismo flujo el proceso de ingeniería de datos y un modelo inicial sugerido. En este primer ejemplo, se muestra un modelo de 5 vecinos más cercanos.

# Crear el pipeline con la regresión por Random Forest
pipeline = Pipeline([
   ('preprocessor', preprocessor),
   ('regressor', RandomForestRegressor(
     n_estimators=10,
     min_samples_split=2,
     min_samples_leaf=2,
     random_state=12345))
])

# Entrenar el pipeline
results = pipeline.fit(ames_train_selected, ames_y_train)

## PREDICCIONES
y_pred = pipeline.predict(ames_x_test)

ames_test = (
  ames_x_test >>
  mutate(Sale_Price_Pred = y_pred, Sale_Price = ames_y_test)
)

(
ames_test >>
  select(_.Sale_Price, _.Sale_Price_Pred)
)
##       Sale_Price  Sale_Price_Pred
## 390       165000        224063.10
## 1235      124000        195821.07
## 2288       75000         91954.24
## 107       206000        132702.51
## 1861      190000        192461.90
## ...          ...              ...
## 116       171000        154051.67
## 398       120500        132702.51
## 1253      146000        224063.10
## 78        125000        130667.37
## 714       110000        125676.00
## 
## [586 rows x 2 columns]

Estas son las predicciones logradas con el modelo inicial.

8.6.1.1 Métricas de desempeño

Se procede en el siguiente paso a cuantificar los errores producidos por la predicción.

Paso 4: Cálculo de métricas de desempeño

pd.options.display.float_format = '{:.2f}'.format

y_obs = ames_test["Sale_Price"]
y_pred = ames_test["Sale_Price_Pred"]

me = np.mean(y_obs - y_pred)
mae = mean_absolute_error(y_obs, y_pred)
mape = mean_absolute_percentage_error(y_obs, y_pred)
mse = mean_squared_error(y_obs, y_pred)
rmse = np.sqrt(mse)
r2 = r2_score(y_obs, y_pred)

n = len(y_obs)  # Número de observaciones
p = 11  # Número de predictores 
r2_adj = 1 - (n - 1) / (n - p - 1) * (1 - r2)

metrics_data = {
    "Metric": ["ME", "MAE", "MAPE", "MSE", "RMSE", "R^2", "R^2 Adj"],
    "Value": [me, mae, mape, mse, rmse, r2, r2_adj]
}

metrics_df = pd.DataFrame(metrics_data)
metrics_df
##     Metric         Value
## 0       ME       -326.52
## 1      MAE      38207.14
## 2     MAPE          0.22
## 3      MSE 3298090728.47
## 4     RMSE      57429.01
## 5      R^2          0.50
## 6  R^2 Adj          0.49

Una manera amigable de dimensionar los errores y el buen desempeño es mediante gráficos

#### Gráficos de desempeño de modelo

(
  ames_test >>
    ggplot(aes(x = "Sale_Price", y = "Sale_Price_Pred")) +
    geom_point() +
    scale_y_continuous(labels = dollar_format(digits=0, big_mark=','), limits = [0, 600000] ) +
    scale_x_continuous(labels = dollar_format(digits=0, big_mark=','), limits = [0, 500000] ) +
    geom_abline(color = "red") +
    coord_equal() +
    labs(
      title = "Comparación entre predicción y observación",
      y = "Predicción",
      x = "Observación")
)
## <Figure Size: (640 x 480)>

(
ames_test >>
  select(_.Sale_Price, _.Sale_Price_Pred) >>
  mutate(error = _.Sale_Price - _.Sale_Price_Pred) >>
  ggplot(aes(x = "error")) +
  geom_histogram(color = "white", fill = "black") +
  geom_vline(xintercept = 0, color = "red") +
  scale_x_continuous(labels=dollar_format(big_mark=',', digits=0)) + 
  ylab("Conteos de clase") + xlab("Errores") +
  ggtitle("Distribución de error")
)
## <Figure Size: (640 x 480)>

(
ames_test >>
  select(_.Sale_Price, _.Sale_Price_Pred) >>
  mutate(error = _.Sale_Price - _.Sale_Price_Pred) >>
  ggplot(aes(sample = "error")) +
  geom_qq(alpha = 0.3) + stat_qq_line(color = "red") +
  scale_y_continuous(labels=dollar_format(big_mark=',', digits = 0)) + 
  xlab("Distribución normal") + ylab("Distribución de errores") +
  ggtitle("QQ-Plot")
)
## <Figure Size: (640 x 480)>

(
ames_test >>
  select(_.Sale_Price, _.Sale_Price_Pred) >>
  mutate(error = _.Sale_Price - _.Sale_Price_Pred) >>
  ggplot(aes(x = "Sale_Price")) +
  geom_linerange(aes(ymin = 0, ymax = "error"), colour = "purple") +
  geom_point(aes(y = "error"), size = 0.05, alpha = 0.5) +
  geom_abline(intercept = 0, slope = 0) +
  scale_x_continuous(labels=dollar_format(big_mark=',', digits=0)) + 
  scale_y_continuous(labels=dollar_format(big_mark=',', digits=0)) +
  xlab("Precio real") + ylab("Error de estimación") +
  ggtitle("Relación entre error y precio de venta")
)
## <Figure Size: (640 x 480)>

8.6.1.2 Validación cruzada

Para determinar cuáles son los hiper-parámetros que funcionan mejor, es necesario realizar experimentos mediante ensayo-error hasta determinar la mejor solución. En cada partición del método de muestreo KFCV se implementan las distintas configuraciones y se calculan predicciones. Con las predicciones hechas en cada fold, se obtienen intervalos de confianza para conocer la variación asociada al modelo a través de los hiper-parámetros implementados.

Usaremos las recetas antes implementadas para ajustar tanto el modelo de regresión como el de clasificación. Exploraremos un conjunto de hiperparámetros para elegir el mejor modelo, sin embargo, para realizar este proceso de forma ágil, se inicializará un flujo de trabajo que se encargue de realizar todos los experimentos deseados y elegir el modelo adecuado.

Paso 5: Creación de grid search y métricas de desempeño

# Definir el objeto K-Fold Cross Validator
k = 5
kf = KFold(n_splits=k, shuffle=True, random_state=42)

param_grid = {
 'max_depth': range(2, 5),
 'min_samples_split': range(2, 8),
 'min_samples_leaf': range(2, 8),
 'max_features': range(1, 5)
}

Una vez definidos los posibles hiperparámetros, procedemos a definir las métricas que serán usadas para cuantificar la bondad del ajuste.

# Definir las métricas de desempeño que deseas calcular como funciones de puntuación

def adjusted_r2_score(y_true, y_pred, n, p):
  r2 = r2_score(y_true, y_pred)
  adjusted_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)
  return adjusted_r2

scoring = {
    'neg_mean_squared_error': make_scorer(mean_squared_error, greater_is_better=False),
    'r2': make_scorer(adjusted_r2_score, greater_is_better=True, 
                      n=np.ceil(len(ames_train_selected)), p=len(ames_train_selected.columns)),
    'neg_mean_absolute_error': make_scorer(mean_absolute_error, greater_is_better=False),
    'mape': make_scorer(mean_absolute_percentage_error, greater_is_better=False)
}

Paso 6: Entrenamiento de modelos con hiperparámetros definidos

Teniendo todos los elementos anteriores listos, se procede con el ajuste de todas las posibles configuraciones del modelo a través de la validación cruzada. Esto permitirá contar con medidas de tendencia central para los resultados de cada uno de las configuraciones y evaluar si estadísticamente hay diferencia entre los mejores resultados.

pipeline = Pipeline([
    ('preprocessor', preprocessor),
    ('regressor', GridSearchCV(
      RandomForestRegressor(), 
      param_grid = param_grid, 
      cv=kf, 
      scoring=scoring, 
      refit='neg_mean_squared_error',
      verbose=2, 
      n_jobs=7,
      error_score='raise')
     )
])

pipeline.fit(ames_train_selected, ames_y_train)

pickle.dump(pipeline, open('models/grid_search_random_forest.pkl', 'wb'))

Paso 7: Análisis de métricas de error e hiperparámetros (Vuelve al paso 3, si es necesario)

Podemos obtener las métricas de los resultados de cada fold:

pipeline = pickle.load(open('models/grid_search_random_forest.pkl', 'rb'))

results_cv = pipeline.named_steps['regressor'].cv_results_

# Convierte los resultados en un DataFrame
pd.set_option('display.max_columns', 500)
results_df = pd.DataFrame(results_cv)
results_df.columns
## Index(['mean_fit_time', 'std_fit_time', 'mean_score_time', 'std_score_time',
##        'param_max_depth', 'param_max_features', 'param_min_samples_leaf',
##        'param_min_samples_split', 'params',
##        'split0_test_neg_mean_squared_error',
##        'split1_test_neg_mean_squared_error',
##        'split2_test_neg_mean_squared_error',
##        'split3_test_neg_mean_squared_error',
##        'split4_test_neg_mean_squared_error',
##        'mean_test_neg_mean_squared_error', 'std_test_neg_mean_squared_error',
##        'rank_test_neg_mean_squared_error', 'split0_test_r2', 'split1_test_r2',
##        'split2_test_r2', 'split3_test_r2', 'split4_test_r2', 'mean_test_r2',
##        'std_test_r2', 'rank_test_r2', 'split0_test_neg_mean_absolute_error',
##        'split1_test_neg_mean_absolute_error',
##        'split2_test_neg_mean_absolute_error',
##        'split3_test_neg_mean_absolute_error',
##        'split4_test_neg_mean_absolute_error',
##        'mean_test_neg_mean_absolute_error', 'std_test_neg_mean_absolute_error',
##        'rank_test_neg_mean_absolute_error', 'split0_test_mape',
##        'split1_test_mape', 'split2_test_mape', 'split3_test_mape',
##        'split4_test_mape', 'mean_test_mape', 'std_test_mape',
##        'rank_test_mape'],
##       dtype='object')
# Puedes seleccionar las columnas de interés, por ejemplo:

summary_df = (
  results_df >>
  select(-_.contains("split._"), -_.contains("time"), -_.params)
)
summary_df
##     param_max_depth param_max_features param_min_samples_leaf  \
## 0                 2                  1                      2   
## 1                 2                  1                      2   
## 2                 2                  1                      2   
## 3                 2                  1                      2   
## 4                 2                  1                      2   
## ..              ...                ...                    ...   
## 427               4                  4                      7   
## 428               4                  4                      7   
## 429               4                  4                      7   
## 430               4                  4                      7   
## 431               4                  4                      7   
## 
##     param_min_samples_split  mean_test_neg_mean_squared_error  \
## 0                         2                    -4970280637.82   
## 1                         3                    -4974587424.48   
## 2                         4                    -5005446361.21   
## 3                         5                    -5019964053.35   
## 4                         6                    -5034494921.64   
## ..                      ...                               ...   
## 427                       3                    -3576638787.58   
## 428                       4                    -3563145452.88   
## 429                       5                    -3584718778.74   
## 430                       6                    -3599654298.71   
## 431                       7                    -3583862726.28   
## 
##      std_test_neg_mean_squared_error  rank_test_neg_mean_squared_error  \
## 0                       502974264.93                               398   
## 1                       520187648.69                               399   
## 2                       596477224.42                               402   
## 3                       475974437.54                               408   
## 4                       462785558.52                               409   
## ..                               ...                               ...   
## 427                     416707863.11                                30   
## 428                     410466613.96                                19   
## 429                     415571196.41                                37   
## 430                     437276292.66                                42   
## 431                     416165232.48                                35   
## 
##      mean_test_r2  std_test_r2  rank_test_r2  \
## 0            0.21         0.02           398   
## 1            0.21         0.02           399   
## 2            0.21         0.03           402   
## 3            0.21         0.02           408   
## 4            0.20         0.01           410   
## ..            ...          ...           ...   
## 427          0.44         0.03            27   
## 428          0.44         0.03            19   
## 429          0.43         0.03            37   
## 430          0.43         0.03            42   
## 431          0.43         0.03            35   
## 
##      mean_test_neg_mean_absolute_error  std_test_neg_mean_absolute_error  \
## 0                            -49246.68                           1474.96   
## 1                            -49247.84                            993.79   
## 2                            -49481.49                           1732.09   
## 3                            -49624.64                            814.17   
## 4                            -49623.45                            724.66   
## ..                                 ...                               ...   
## 427                          -39875.82                           1197.55   
## 428                          -39734.62                           1124.42   
## 429                          -39856.54                           1130.38   
## 430                          -39917.63                           1252.06   
## 431                          -39806.33                           1282.96   
## 
##      rank_test_neg_mean_absolute_error  mean_test_mape  std_test_mape  \
## 0                                  398           -0.31           0.02   
## 1                                  399           -0.31           0.01   
## 2                                  404           -0.32           0.01   
## 3                                  409           -0.32           0.01   
## 4                                  408           -0.32           0.01   
## ..                                 ...             ...            ...   
## 427                                 29           -0.24           0.01   
## 428                                  6           -0.24           0.01   
## 429                                 26           -0.24           0.01   
## 430                                 33           -0.24           0.01   
## 431                                 14           -0.24           0.01   
## 
##      rank_test_mape  
## 0               399  
## 1               398  
## 2               402  
## 3               406  
## 4               409  
## ..              ...  
## 427              27  
## 428               6  
## 429              28  
## 430              36  
## 431              20  
## 
## [432 rows x 16 columns]

En la siguiente gráfica observamos las distintas métricas de error asociados a los hiperparámetros elegidos:

(
  summary_df >>
  ggplot(aes(x = "param_max_features", y = "mean_test_r2", fill = "param_max_depth")) +
  geom_point() +
  ggtitle("Parametrización de Random Forest vs R^2") +
  xlab("Parámetro: Número de features por árbol") +
  ylab("R^2 promedio")
)
## <Figure Size: (640 x 480)>

Seleccionando la mejor parametrización, se puede definir la siguiente gráfica para comparar la diferencia entre hiperparámetros a un nivel más granular.

(
  summary_df >>
  select(_.param_max_depth, _.param_max_features, _.param_min_samples_leaf, 
         _.param_min_samples_split, _.mean_test_r2) >>
  pivot_longer(
    cols = ["param_max_depth", "param_max_features", "param_min_samples_leaf", "param_min_samples_split"],
    names_to="parameter",
    values_to="value") >>
  ggplot(aes(x = "value", y = "mean_test_r2")) +
  geom_point(size = 1, ) +
  facet_wrap("~parameter", scales = "free_x") +
  xlab("Parameter value") +
  ylab("R^2 promedio") +
  ggtitle("Parametrización de Random Forest vs R^2")
)
## <Figure Size: (640 x 480)>

Paso 8: Selección de modelo a usar

Habiendo realizado un análisis de los hiperparámetros, se procede a elegir el mejor modelo. Esto a veces a muy evidente y otras no. En cualquier caso, puede automatizarse la extracción del mejor modelo:

best_params = pipeline.named_steps['regressor'].best_params_
best_params
## {'max_depth': 4, 'max_features': 4, 'min_samples_leaf': 2, 'min_samples_split': 5}

Sabiendo cuáles son los mejores hiperparámetros, se procede a extraer el modelo que DEBEREMOS AJUSTAR A TODOS LOS DATOS DE ENTRENAMIENTO.

best_estimator = pipeline.named_steps['regressor'].best_estimator_
best_estimator
## RandomForestRegressor(max_depth=4, max_features=4, min_samples_leaf=2,
##                       min_samples_split=5)

Paso 9: Ajuste de modelo final con todos los datos (Vuelve al paso 2, si es necesario)

Ahora obtendremos el modelo que mejor desempeño tiene y haremos las predicciones del conjunto de prueba con este modelo.

Es importante volver a hacer el ajuste con el modelo elegido.

final_rf_pipeline = Pipeline([
   ('preprocessor', preprocessor),
   ('regressor', best_estimator)
])

# Entrenar el pipeline
final_rf_pipeline.fit(ames_train_selected, ames_y_train)
## Pipeline(steps=[('preprocessor',
##                  ColumnTransformer(remainder='passthrough',
##                                    transformers=[('scaler', StandardScaler(),
##                                                   ['Full_Bath', 'Half_Bath',
##                                                    'Second_Flr_SF']),
##                                                  ('onehotencoding',
##                                                   OneHotEncoder(drop='first',
##                                                                 sparse_output=False),
##                                                   ['Overall_Cond'])],
##                                    verbose_feature_names_out=False)),
##                 ('regressor',
##                  RandomForestRegressor(max_depth=4, max_features=4,
##                                        min_samples_leaf=2,
##                                        min_samples_split=5))])

Este último objeto creado es el modelo final entrenado, el cual contiene toda la información del pre-procesamiento de datos, por lo que en caso de ponerse en producción, sólo se necesita de los nuevos datos y de este último elemento para poder realizar nuevas predicciones.

Paso 10: Validar poder predictivo con datos de prueba

Imaginemos por un momento que pasa un mes de tiempo desde que hicimos nuestro modelo, es hora de ponerlo a prueba prediciendo valores de nuevos elementos:

## Predicciones finales
y_pred_rf = final_rf_pipeline.predict(ames_x_test)

results_reg = (
  ames_x_test >>
  mutate(final_rf_pred = y_pred_rf, Sale_Price = ames_y_test) >>
  select(_.Sale_Price, _.final_rf_pred)
)
results_reg
##       Sale_Price  final_rf_pred
## 390       165000      217381.45
## 1235      124000      170935.03
## 2288       75000      157630.42
## 107       206000      131537.36
## 1861      190000      256303.03
## ...          ...            ...
## 116       171000      153121.89
## 398       120500      131537.36
## 1253      146000      217381.45
## 78        125000      135611.99
## 714       110000      146506.48
## 
## [586 rows x 2 columns]

Métricas de desempeño

Ahora para calcular las métricas de desempeño usaremos la paquetería MLmetrics. Es posible definir nuestro propio conjunto de métricas que deseamos reportar creando el objeto metric_set:

me = np.mean(y_obs - y_pred_rf)
mae = mean_absolute_error(y_obs, y_pred_rf)
mape = mean_absolute_percentage_error(y_obs, y_pred_rf)
mse = mean_squared_error(y_obs, y_pred_rf)
rmse = np.sqrt(mse)
r2 = r2_score(y_obs, y_pred_rf)
r2_adj = adjusted_r2_score(y_true = y_obs, y_pred = y_pred_rf,
  n=np.ceil(len(ames_train_selected)), p=len(ames_train_selected.columns))

metrics_data = {
    "Metric": ["ME", "MAE", "MAPE", "MSE", "RMSE", "R^2", "R^2 Adj"],
    "Value": [me, mae, mape, mse, rmse, r2, r2_adj]
}

metrics_df = pd.DataFrame(metrics_data)
metrics_df
##     Metric         Value
## 0       ME       -926.42
## 1      MAE      42217.89
## 2     MAPE          0.25
## 3      MSE 3860886848.32
## 4     RMSE      62136.04
## 5      R^2          0.41
## 6  R^2 Adj          0.41
(
  results_reg >>
    ggplot(aes(x = "final_rf_pred", y = "Sale_Price")) +
    geom_point() +
    geom_abline(color = "red") +
    xlab("Prediction") +
    ylab("Observation") +
    ggtitle("Comparisson")
)
## <Figure Size: (640 x 480)>

8.6.1.3 Importancia de variables

Antes de pasar al siguiente paso, es importante validar que hayamos hecho un uso correcto de las variables predictivas. En este momento es posible detectar variables que no estén aportando valor o variables que no debiéramos estar usando debido a que cometeríamos data leakage. Para enfrentar esto, ayuda estimar y ordenar el valor de importancia de cada variable en el modelo.

importance = np.zeros(ames_x_test[columnas_seleccionadas].shape[1])

# Realiza el procedimiento de permutación
for i in range(ames_x_test[columnas_seleccionadas].shape[1]):
    ames_x_test_permuted = ames_x_test[columnas_seleccionadas].copy()
    ames_x_test_permuted.iloc[:, i] = shuffle(ames_x_test_permuted.iloc[:, i], random_state=42)  
    # Permuta una característica
    y_pred_permuted = final_rf_pipeline.predict(ames_x_test_permuted)
    mse_permuted = mean_squared_error(ames_y_test, y_pred_permuted)
    importance[i] = mse - mse_permuted

# Calcula la importancia relativa
importance = importance / importance.sum()
importance
## array([0.59119845, 0.10125023, 0.18282385, 0.12472747])
importance_df = pd.DataFrame({
  'Variable': columnas_seleccionadas, 
  'Importance': importance
  })

# Crea la gráfica de barras
(
  importance_df >>
  ggplot(aes(x= 'reorder(Variable, Importance)', y='Importance')) + 
  geom_bar(stat='identity', fill='blue', color = "black") + 
  labs(title='Importancia de las Variables', x='Variable', y='Importancia') +
  coord_flip()
)
## <Figure Size: (640 x 480)>

La gráfica anterior muestra la importancia de una variable cuando se lleva a cabo una permutación. Dado que este resultado fue aleatorio, resulta vital contar con un conjunto de permutaciones que permitan conocer la caída promedio y desviación esperada en el desempeño cuando se elimina una variable o esta se vuelve inservible.

n_permutations = 50
performance_losses = []

for i in range(ames_x_test[columnas_seleccionadas].shape[1]):
    loss = []
    for j in range(n_permutations):
        ames_x_test_permuted = ames_x_test[columnas_seleccionadas].copy()
        ames_x_test_permuted.iloc[:, i] = np.random.permutation(ames_x_test_permuted.iloc[:, i])
        y_pred_permuted = final_rf_pipeline.predict(ames_x_test_permuted)
        mse_permuted = mean_squared_error(ames_y_test, y_pred_permuted)
        loss.append(mse_permuted)
    performance_losses.append(loss)

performance_losses = performance_losses/np.sum(performance_losses, axis=0)
mean_losses = np.mean(performance_losses, axis=1)
std_losses = np.std(performance_losses, axis=1)

importance_df = pd.DataFrame({
  'Variable': columnas_seleccionadas, 
  'Mean_Loss': mean_losses, 
  'Std_Loss': std_losses
  })
(
  importance_df >>
  mutate(
    ymin = _.Mean_Loss - _.Std_Loss,
    ymax = _.Mean_Loss + _.Std_Loss) >>
  ggplot(aes(x = 'reorder(Variable, Mean_Loss)', y = "Mean_Loss")) +
  geom_errorbar(aes(ymin='ymin', ymax='ymax'),
    width=0.2, position=position_dodge(0.9)) +
  geom_point(alpha = 0.65) +
  labs(title='Importancia de las Variables', x='Variable', y='Importancia') +
  coord_flip()
)
## <Figure Size: (640 x 480)>

importances = final_rf_pipeline.named_steps['regressor'].feature_importances_
columns = final_rf_pipeline["preprocessor"].get_feature_names_out()

(
  pd.DataFrame({'feature': columns, 'Importance': importances}) >>
  ggplot(aes(x= 'reorder(feature, Importance)', y='Importance')) + 
  geom_bar(stat='identity', fill='blue', color = "black") + 
  labs(title='Importancia de las Variables', x='Variable', y='Importancia') +
  coord_flip()
)
## <Figure Size: (640 x 480)>

Si en la gráfica anterior notamos algo raro en cuanto a la(s) variable(s) más importante(s) y la factibilidad de conseguirla(s) o usarla(s)… ¡¡HAY QUE EMPEZAR DESDE CERO SIN CONSIDERAR ESA VARIABLE!!