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()
## • Search for functions across packages at https://www.tidymodels.org/find/

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

Modelo de Random Forest

rf_tune <- rand_forest(
  mtry = tune(),
  trees = 1000,
  min_n = tune()
) %>%
  set_mode("regression") %>%
  set_engine("ranger", importance = "permutation")

Workflow

tune_wf <- workflow() %>%
  add_recipe(rf_recipe) %>%
  add_model(rf_tune)

Entrenamiento del Modelo

# Define the metrics, ensuring "mae" is included
metrics <- metric_set(mae, rmse, mape)
set.seed(8577)
doParallel::registerDoParallel()
ranger_tune <-tune_grid(tune_wf,
    resamples = p_folds,
    metrics=metrics,
    grid = 11
  )
## i Creating pre-processing data to finalize unknown parameter: mtry
show_best(ranger_tune, metric = "mae")
## # A tibble: 5 × 8
##    mtry min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
## 1     5     6 mae     standard    104.     3   1.37  Preprocessor1_Model06
## 2     4    10 mae     standard    106.     3   1.36  Preprocessor1_Model05
## 3     2     4 mae     standard    106.     3   0.945 Preprocessor1_Model10
## 4     3    18 mae     standard    110.     3   1.50  Preprocessor1_Model01
## 5     4    24 mae     standard    111.     3   1.67  Preprocessor1_Model03
ranger_tune %>%
  collect_metrics() %>%
  filter(.metric == "mae") %>%
  dplyr::select(mean, min_n, mtry) %>%
  pivot_longer(min_n:mtry,
    values_to = "value",
    names_to = "parameter"
  ) %>%
  ggplot(aes(value, mean, color = parameter)) +
  geom_point(show.legend = FALSE) +
  facet_wrap(~parameter, scales = "free_x") +
  labs(x = NULL, y = "mae")

show_best(ranger_tune, "mae")
## # A tibble: 5 × 8
##    mtry min_n .metric .estimator  mean     n std_err .config              
##   <int> <int> <chr>   <chr>      <dbl> <int>   <dbl> <chr>                
## 1     5     6 mae     standard    104.     3   1.37  Preprocessor1_Model06
## 2     4    10 mae     standard    106.     3   1.36  Preprocessor1_Model05
## 3     2     4 mae     standard    106.     3   0.945 Preprocessor1_Model10
## 4     3    18 mae     standard    110.     3   1.50  Preprocessor1_Model01
## 5     4    24 mae     standard    111.     3   1.67  Preprocessor1_Model03
best_mae <- select_best(ranger_tune, "mae")
best_mae
## # A tibble: 1 × 3
##    mtry min_n .config              
##   <int> <int> <chr>                
## 1     5     6 Preprocessor1_Model06
final_rf <- finalize_workflow(
  tune_wf,
  best_mae
)
final_rf
## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: rand_forest()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 1 Recipe Step
## 
## • step_normalize()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Random Forest Model Specification (regression)
## 
## Main Arguments:
##   mtry = 5
##   trees = 1000
##   min_n = 6
## 
## Engine-Specific Arguments:
##   importance = permutation
## 
## Computational engine: ranger

Predicción en TEST

# Define the metrics, ensuring "mae" is included
metrics <- metric_set(mae, rmse, mape)
final_res <- last_fit(final_rf, p_split, metrics=metrics)
## Warning: package 'ranger' was built under R version 4.3.3
collect_metrics(final_res)
## # A tibble: 3 × 4
##   .metric .estimator .estimate .config             
##   <chr>   <chr>          <dbl> <chr>               
## 1 mae     standard       102.  Preprocessor1_Model1
## 2 rmse    standard       127.  Preprocessor1_Model1
## 3 mape    standard        11.6 Preprocessor1_Model1
# Obtener los valores verdaderos y predichos
predictions <- final_res %>%
  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 variables

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