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.
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
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")
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
# 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 …
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
# 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)
# 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
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()
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")