Random Forest en AirQuality Dataset

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 tidymodels_prefer() to resolve common conflicts.

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³).

data <- read.csv("https://raw.githubusercontent.com/data-datum/datasets_ml/refs/heads/main/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
p_split <- initial_split(data, prop = 0.8, strata = PT08.S5.O3.)
train_data <- training(p_split)
test_data <- testing(p_split)

Los datos de TRAIN voy a dividirlos en 3-folds

p_folds <- vfold_cv(train_data, v=3)
# 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

Especificamos el modelo de ML

En esta etapa se especifica el modelo que vamos a implementar, en este caso, Random Forest Baseline.

rf_spec <- rand_forest() %>% 
  set_engine("ranger", importance = "permutation") %>% 
  set_mode("regression")

Workflow

Inicializo el workflow para trabajar de manera ordenada, luego visualizo los pasos.

rf_wf <- workflow() %>%
    add_recipe(rf_recipe) %>%
    add_model(rf_spec)
rf_wf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_normalize()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (regression)
## 
## Engine-Specific Arguments:
##   importance = permutation
## 
## Computational engine: ranger

Entrenamiento del Modelo

# Define the metrics, ensuring "mae" is included
metrics <- metric_set(mae, rmse, mape)
set.seed(123)
rf_res <- rf_wf %>%
    fit_resamples(
        p_folds,
        metrics = metrics, 
        control = control_resamples(save_pred = TRUE)
    )
## Warning: package 'ranger' was built under R version 4.3.3
glimpse(rf_res)
## Rows: 3
## Columns: 5
## $ splits       <list> [<vfold_split[452 x 227 x 679 x 6]>], [<vfold_split[453 …
## $ id           <chr> "Fold1", "Fold2", "Fold3"
## $ .metrics     <list> [<tbl_df[3 x 4]>], [<tbl_df[3 x 4]>], [<tbl_df[3 x 4]>]
## $ .notes       <list> [<tbl_df[0 x 3]>], [<tbl_df[0 x 3]>], [<tbl_df[0 x 3]>]
## $ .predictions <list> [<tbl_df[227 x 4]>], [<tbl_df[226 x 4]>], [<tbl_df[226 x …

Estos son los resultados de TRAIN

rf_res %>%
  collect_metrics()
## # A tibble: 3 × 6
##   .metric .estimator  mean     n std_err .config             
##   <chr>   <chr>      <dbl> <int>   <dbl> <chr>               
## 1 mae     standard   107.      3   0.639 Preprocessor1_Model1
## 2 mape    standard    12.2     3   0.240 Preprocessor1_Model1
## 3 rmse    standard   136.      3   1.26  Preprocessor1_Model1

Finalizamos el workflow

# Select the best parameters based on MAE and finalize the workflow
best_params <- select_best(rf_res, "mae")
final_wf <- finalize_workflow(rf_wf, best_params)

Predicción en TEST set

# Define the metrics, ensuring "mae" is included
metrics <- metric_set(mae, rmse, mape)
final_rs <- last_fit(final_wf, p_split, metrics=metrics)
final_rs
## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits            id               .metrics .notes   .predictions .workflow 
##   <list>            <chr>            <list>   <list>   <list>       <list>    
## 1 <split [679/172]> train/test split <tibble> <tibble> <tibble>     <workflow>
final_rs %>%
  collect_metrics()
## # A tibble: 3 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 mae     standard       101.  Preprocessor1_Model1
## 2 rmse    standard       125.  Preprocessor1_Model1
## 3 mape    standard        11.7 Preprocessor1_Model1
final_rs %>%
  collect_predictions()
## # A tibble: 172 × 5
##    id               .pred  .row PT08.S5.O3. .config             
##    <chr>            <dbl> <int>       <int> <chr>               
##  1 train/test split 1202.     1        1268 Preprocessor1_Model1
##  2 train/test split  837.     2         972 Preprocessor1_Model1
##  3 train/test split  692.     9         620 Preprocessor1_Model1
##  4 train/test split 1591.    24        1409 Preprocessor1_Model1
##  5 train/test split 1335.    27        1285 Preprocessor1_Model1
##  6 train/test split  697.    29         552 Preprocessor1_Model1
##  7 train/test split  425.    33         384 Preprocessor1_Model1
##  8 train/test split  902.    36         748 Preprocessor1_Model1
##  9 train/test split 1079.    42        1061 Preprocessor1_Model1
## 10 train/test split 1285.    43        1139 Preprocessor1_Model1
## # ℹ 162 more rows

Gráfico de Valores Predichos vs. Valores reales

library(ggplot2)

# Obtener los valores verdaderos y predichos
predictions <- final_rs %>%
  collect_predictions()

# Crear el gráfico de valores reales vs valores predichos
ggplot(predictions, aes(x = PT08.S5.O3., y = .pred)) +
  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()

Importancia de las variables

Vamos a plotear las variables más importantes del MODELO, mediante la librería vip. Más información sobre esta librería en https://koalaverse.github.io/vip/

library(vip)
# Fit the model to the entire training data to get the final model
final_fit <- rf_wf %>%
  fit(data = train_data)
# Plot variable importance
vip(pull_workflow_fit(final_fit)$fit, geom = "point")