Random Forest Regression sin Tunning de Hiperparámetros

Evaluación de la calidad del aire con dataset AirQuality

El conjunto de datos Air Quality del repositorio de UCI (https://archive.ics.uci.edu/dataset/360/air+quality) contiene mediciones de la calidad del aire tomadas en una estación de monitoreo de calidad del aire en una zona suburbana de Milán, Italia, en el año 2004. Este dataset es ampliamente utilizado para experimentos de análisis de series temporales, modelos de predicción y problemas de regresión.

Variables del dataset:

Date (Fecha): La fecha de la medición en el formato DD/MM/YYYY.
Time (Hora): La hora de la medición en el formato de 24 horas HH.MM.SS.
CO(GT): Concentración de monóxido de carbono (CO) en mg/m³. Este valor fue registrado por un sensor de gas (GT).
PT08.S1(CO): Señal del sensor 1, que es sensible al CO. Este sensor responde principalmente al CO y está representado en unidades arbitrarias de la señal del sensor.
NMHC(GT): Hidrocarburos no metánicos (NMHC) en microgramos por metro cúbico (µg/m³).
C6H6(GT): Concentración de benceno en microgramos por metro cúbico (µg/m³).
PT08.S2(NMHC): Señal del sensor 2, que responde principalmente a hidrocarburos no metánicos.
NOx(GT): Concentración de óxidos de nitrógeno (NOx) en partes por mil millones (ppb).
PT08.S3(NOx): Señal del sensor 3, que responde principalmente a los óxidos de nitrógeno.
NO2(GT): Concentración de dióxido de nitrógeno (NO₂) en partes por mil millones (ppb).
PT08.S4(NO2): Señal del sensor 4, sensible al dióxido de nitrógeno.
PT08.S5(O3): Señal del sensor 5, sensible al ozono (O₃).
T (Temperatura): Temperatura ambiental en grados Celsius (°C).
RH (Humedad Relativa): Humedad relativa en porcentaje (%).
AH (Humedad Absoluta): Humedad absoluta en gramos por metro cúbico (g/m³).

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(tidymodels)

## ── Attaching packages ────────────────────────────────────── tidymodels 1.1.1 ──
## ✔ broom        1.0.5      ✔ rsample      1.2.1 
## ✔ dials        1.2.1      ✔ tune         1.1.2 
## ✔ infer        1.0.6      ✔ workflows    1.1.4 
## ✔ modeldata    1.3.0      ✔ workflowsets 1.0.1 
## ✔ parsnip      1.2.0      ✔ yardstick    1.3.0 
## ✔ recipes      1.0.10     
## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ scales::discard() masks purrr::discard()
## ✖ dplyr::filter()   masks stats::filter()
## ✖ recipes::fixed()  masks stringr::fixed()
## ✖ dplyr::lag()      masks stats::lag()
## ✖ yardstick::spec() masks readr::spec()
## ✖ recipes::step()   masks stats::step()
## • Use suppressPackageStartupMessages() to eliminate package startup messages

# Cargar el dataset
data <- read.csv("AirQualityUCI.csv", sep = ";", dec = ",", na.strings = c("NA", "-200"))

# Verificar nombres de columnas
names(data)

##  [1] "Date"          "Time"          "CO.GT."        "PT08.S1.CO."  
##  [5] "NMHC.GT."      "C6H6.GT."      "PT08.S2.NMHC." "NOx.GT."      
##  [9] "PT08.S3.NOx."  "NO2.GT."       "PT08.S4.NO2."  "PT08.S5.O3."  
## [13] "T"             "RH"            "AH"            "X"            
## [17] "X.1"

# Seleccionar las variables de interés y eliminar valores faltantes
data <- data %>%
  select(PT08.S5.O3., T, RH, NOx.GT., NO2.GT., NMHC.GT.) %>%
  drop_na()

set.seed(123)  # Para reproducibilidad
split <- initial_split(data, prop = 0.8, strata = PT08.S5.O3.)
train_data <- training(split)
test_data <- testing(split)

# Crear la receta de preprocesamiento
rf_recipe <- recipe(PT08.S5.O3. ~ T + RH + NOx.GT. + NO2.GT. + NMHC.GT., data = train_data) %>%
  step_normalize(all_predictors())  # Estandarizar las variables predictoras

# Especificar el modelo de Random Forest
rf_model <- rand_forest(mode = "regression",
  trees = 500) %>%
  set_engine("ranger")

Workflow

El workflow nos permite trabajar de manera ordenada en MACHINE LEARNING, y es una buena práctica utilizarlo siempre que tengamos más de un paso de pre-procesamiento.

Agregamos la receta mediante add_recipe(), y luego,
agregamos el modelo.

Luego imprimo el workflow

# Crear el flujo de trabajo
rf_workflow <- workflow() %>%
  add_model(rf_model) %>%
  add_recipe(rf_recipe)

# Entrenar el modelo
rf_fit <- rf_workflow %>%
  fit(data = train_data)

rf_fit

## ══ Workflow [trained] ══════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_normalize()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Ranger result
## 
## Call:
##  ranger::ranger(x = maybe_data_frame(x), y = y, num.trees = ~500,      num.threads = 1, verbose = FALSE, seed = sample.int(10^5,          1)) 
## 
## Type:                             Regression 
## Number of trees:                  500 
## Sample size:                      679 
## Number of independent variables:  5 
## Mtry:                             2 
## Target node size:                 5 
## Variable importance mode:         none 
## Splitrule:                        variance 
## OOB prediction error (MSE):       17676.56 
## R squared (OOB):                  0.8904308

Resultados en TRAIN

En este paso lo que vamos a hacer es agregar la columna de la predicción a los datos de TRAIN, mediante la función bind_cols. Esto nos va a permitir ver la diferencia entre la predicción y el valor verdadero.

results <- train_data %>%
    bind_cols(predict(rf_fit, train_data) %>%
                  rename(.pred_tree = .pred))

Ahora vamos a ver las primeras filas de este nuevo dataset que hemos armado.

results %>%
  head()

##   PT08.S5.O3.    T   RH NOx.GT. NO2.GT. NMHC.GT. .pred_tree
## 1         733 11.3 56.8      62      77       31   721.7349
## 2         730 10.7 60.0      62      76       31   733.9542
## 3         445 10.1 60.5      21      34       14   419.0954
## 4         422 11.0 56.2      16      28        8   415.4827
## 5         472 10.5 58.1      34      48       16   534.7845
## 6         730 10.2 59.6      98      82       29   768.6969

Vamos a calcular las metricas correspondientes del modelo. Para ello utilizamos la función metrics() disponible tambien en tidymodels. Agregamos los el dataset, y le decimos cual es el valor verdadero y cual es la predicción.

metrics(results, truth = PT08.S5.O3., estimate = .pred_tree)

## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard      59.4  
## 2 rsq     standard       0.980
## 3 mae     standard      46.3

Error de estimacion del modelo en TEST

Esta etapa es la predicción del modelo de ML, y es el resultado que generalmente se reporta. IMPORTANTE: los resultados se deben reportar SI O SI, puede o no ir acompañado de los resultados de TRAIN.

results_test <- test_data %>%
    bind_cols(predict(rf_fit, test_data) %>%
                  rename(.pred_rf = .pred))

results_test %>%
  head()

##   PT08.S5.O3.    T   RH NOx.GT. NO2.GT. NMHC.GT.  .pred_rf
## 1        1268 13.6 48.9     166     113      150 1173.0146
## 2         972 13.3 47.7     103      92      112  818.2700
## 3         620 10.7 59.7      45      60       24  694.4153
## 4        1409 10.3 64.2     281     151      307 1573.8214
## 5        1285  9.1 64.0     240     136      197 1322.1775
## 6         552  8.2 60.8      47      53       26  703.0286

Métricas

Ahora vemos las métricas.

metrics(results_test, truth = PT08.S5.O3., estimate = .pred_rf)

## # A tibble: 3 × 3
##   .metric .estimator .estimate
##   <chr>   <chr>          <dbl>
## 1 rmse    standard     125.   
## 2 rsq     standard       0.901
## 3 mae     standard     101.

summary(data$PT08.S5.O3.)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     263     755     996    1036    1311    2359

Gráfico de Valores Predichos vs. Valores reales

library(ggplot2)

# Crear el gráfico de valores reales vs valores predichos
ggplot(results_test, aes(x = PT08.S5.O3., y = .pred_rf)) +
  geom_point(color = "blue", alpha = 0.6) +  # Puntos azules para las observaciones
  geom_abline(slope = 1, intercept = 0, color = "red", linetype = "dashed") +  # Línea de igualdad
  labs(
    title = "Predicted vs Real Values of PT08.S5(O3)",
    x = "Real Values (PT08.S5(O3))",
    y = "Predicted Values"
  ) +
  theme_minimal()