TF.P2.Machine Learning en Aplicaciones Espaciales. Docente: Roxana Villafañe

1. Objetivo y contexto:

El objetivo es construir un modelo de clasificación binaria que estime la probabilidad de que una medición horaria de PM2.5 en el área de Atenas exceda un umbral de referencia, utilizando variables meteorológicas y espaciales (latitud/longitud), junto con información temporal (mes, hora, día de la semana).

La motivación es que la contaminación atmosférica presenta variación espacial (zonas) y temporal (hora/estación), por lo que el enfoque de aprendizaje supervisado es adecuado para capturar relaciones no lineales entre condiciones meteorológicas y eventos de alta concentración.

2. Dataset:

Se utiliza el dataset “Air Quality Monitoring in European Cities” (Kaggle)

Se trabaja con el archivo correspondiente a Atenas, con observaciones horarias y variables:

• Contaminantes: PM2.5, PM10, NO2, O3
• Meteorología: temperatura, humedad relativa, viento (componentes u y v), precipitación
• Espaciales: latitud, longitud
• Temporal: fecha/hora (derivando mes, hora y día de la semana)
• Identificación de “estación” o punto de medición

Las variables wind_speed_u y wind_speed_v son las componentes del vector viento (ejes ortogonales). Se conservan ambas porque la dispersión del material particulado depende no solo de la intensidad sino también de la dirección del viento.

3. Librerías:

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(janitor)

## 
## Attaching package: 'janitor'

## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(tidymodels)

## ── Attaching packages ────────────────────────────────────── tidymodels 1.4.1 ──

## ✔ broom        1.0.12     ✔ rsample      1.3.1 
## ✔ dials        1.4.2      ✔ tailor       0.1.0 
## ✔ ggplot2      4.0.0      ✔ tidyr        1.3.1 
## ✔ infer        1.1.0      ✔ tune         2.0.1 
## ✔ modeldata    1.5.1      ✔ workflows    1.3.0 
## ✔ parsnip      1.4.1      ✔ workflowsets 1.1.1 
## ✔ purrr        1.2.1      ✔ yardstick    1.3.2 
## ✔ recipes      1.3.1

## ── Conflicts ───────────────────────────────────────── tidymodels_conflicts() ──
## ✖ purrr::discard() masks scales::discard()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ✖ recipes::step()  masks stats::step()

library(ggplot2)
library(vip)

## 
## Attaching package: 'vip'

## The following object is masked from 'package:utils':
## 
##     vi

4 Carga y estandarización de nombres:

dataset_qa_air_athens <- read.csv("data/athens_data.csv")

quality_air_athens <- dataset_qa_air_athens %>% clean_names()

5 Limpieza de datos:

Duplicados

dup_n <- sum(duplicated(quality_air_athens))
dup_n

## [1] 0

Si existieran duplicados, se eliminan para evitar que el modelo “memorice” observaciones repetidas:

quality_air_athens <- quality_air_athens %>%
  distinct()

Valores faltantes (NA)

Se contabilizan faltantes por variable:

na_by_col <- sapply(quality_air_athens, function(x) sum(is.na(x)))
na_by_col

##                date            latitude           longitude        station_name 
##                   0                   0                   0                   0 
##        wind_speed_u        wind_speed_v       dewpoint_temp           soil_temp 
##                8573                8573                8573                8573 
## total_percipitation     vegitation_high      vegitation_low                temp 
##                8573                8573                8573                7565 
##   relative_humidity                pm10               pm2_5                 no2 
##                8573                   0                   0              143293 
##                  o3                code                  id 
##              148064                   0                   0

No se eliminan filas masivamente para no perder volumen; en el modelado se aplicará imputación mediana para predictores numéricos y manejo de niveles nuevos/desconocidos para variables categóricas.

Conversión de fecha y variables temporales

quality_air_athens <- quality_air_athens %>%
  mutate(date = as_datetime(date),
    month = factor(month(date)),
    hour = factor(hour(date)),
    weekday = factor(wday(date, label = TRUE)))

Tratamiento de valores extremos en PM2.5

Se observan valores máximos muy altos (potencialmente errores o eventos extremos). Para reducir influencia de outliers, se elimina PM2.5 > 500.

summary(quality_air_athens$pm2_5)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    8.00   11.65   15.07   17.76 1748.72

quality_air_athens <- quality_air_athens %>%
  filter(pm2_5 <= 500)

summary(quality_air_athens$pm2_5)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    8.00   11.65   15.04   17.76  494.20

Justificación: el objetivo es un modelo estable; outliers extremos pueden distorsionar la separación entre clases y el entrenamiento.

6 Variable objetivo (clasificación):

Se define como “alta contaminación” si PM2.5 > 15 µg/m³ (valor establecido por la OMS). Este umbral se usa como referencia práctica para discriminar episodios de mayor exposición (y, además, genera un balance de clases razonable para clasificación).

quality_air_athens <- quality_air_athens %>%
  mutate(
    pollution_high = if_else(pm2_5 > 15, "yes", "no"),
    pollution_high = factor(pollution_high))
prop.table(table(quality_air_athens$pollution_high))

## 
##        no       yes 
## 0.6654677 0.3345323

7 Selección de observaciones recientes:

Para simular un escenario predictivo más realista (entrenar con datos recientes), se trabaja con las 100.000 observaciones más actuales.

quality_air_athens <- quality_air_athens %>%
  arrange(desc(date))

athens_model <- quality_air_athens %>%
  select(
    pollution_high,
    latitude,
    longitude,
    temp,
    relative_humidity,
    wind_speed_u,
    wind_speed_v,
    total_percipitation,
    month,
    hour,
    weekday,
    station_name) %>%
  filter(!is.na(pollution_high)) %>%
  slice_head(n = 100000)

dim(athens_model)

## [1] 100000     12

prop.table(table(athens_model$pollution_high))

## 
##      no     yes 
## 0.70953 0.29047

8 Partición Train/Test:

Se separa en entrenamiento y test usando estratificación para mantener proporciones de clase.

set.seed(123)
aq_split <- initial_split(athens_model, strata = pollution_high)
aq_train <- training(aq_split)
aq_test  <- testing(aq_split)

aq_train %>% count(pollution_high)

##   pollution_high     n
## 1             no 53214
## 2            yes 21785

aq_test  %>% count(pollution_high)

##   pollution_high     n
## 1             no 17739
## 2            yes  7262

9 Preprocesamiento:

• Imputación mediana para numéricos.
• Manejo de niveles nuevos/unknown para categóricas.
• Dummies para modelos que lo requieren (RF y LR en este pipeline).

aq_recipe <-
  recipe(pollution_high ~ ., data = aq_train) %>%
  step_impute_median(all_numeric_predictors()) %>%
  step_novel(all_nominal_predictors()) %>%
  step_unknown(all_nominal_predictors()) %>%
  step_dummy(all_nominal_predictors())

10 Modelo final: Random Forest:

Se utiliza Random Forest por su capacidad de: - capturar relaciones no lineales, - manejar interacciones entre variables, - ser robusto ante ruido.

rf_spec <- rand_forest(
  mtry = 6,
  trees = 500,
  min_n = 10) %>%
  set_engine("ranger", importance = "impurity") %>%
  set_mode("classification")

aq_workflow <- workflow() %>%
  add_recipe(aq_recipe) %>%
  add_model(rf_spec)

11 Validación cruzada (5-fold):

set.seed(234)
aq_folds <- vfold_cv(aq_train, v = 5, strata = pollution_high)

rf_rs <- fit_resamples(
  aq_workflow,
  resamples = aq_folds,
  control = control_resamples(save_pred = TRUE),
  metrics = metric_set(accuracy, roc_auc, sens, spec, bal_accuracy))

rf_metrics <- rf_rs %>% collect_metrics()
rf_metrics

## # A tibble: 5 × 6
##   .metric      .estimator  mean     n std_err .config        
##   <chr>        <chr>      <dbl> <int>   <dbl> <chr>          
## 1 accuracy     binary     0.883     5 0.00149 pre0_mod0_post0
## 2 bal_accuracy binary     0.824     5 0.00189 pre0_mod0_post0
## 3 roc_auc      binary     0.951     5 0.00103 pre0_mod0_post0
## 4 sens         binary     0.965     5 0.00142 pre0_mod0_post0
## 5 spec         binary     0.684     5 0.00328 pre0_mod0_post0

12 Curva ROC (Random Forest):

preds_rf <- rf_rs %>% collect_predictions()

preds_rf %>%
  roc_curve(pollution_high, .pred_yes, event_level = "second") %>%
  autoplot() +
  labs(title = "Curva ROC - Random Forest")

13 Evaluación final en conjunto de test

final_fit <- last_fit(
  aq_workflow,
  aq_split,
  metrics = metric_set(accuracy, roc_auc, sens, spec, bal_accuracy))

final_fit %>% collect_metrics()

## # A tibble: 5 × 4
##   .metric      .estimator .estimate .config        
##   <chr>        <chr>          <dbl> <chr>          
## 1 accuracy     binary         0.885 pre0_mod0_post0
## 2 sens         binary         0.965 pre0_mod0_post0
## 3 spec         binary         0.691 pre0_mod0_post0
## 4 bal_accuracy binary         0.828 pre0_mod0_post0
## 5 roc_auc      binary         0.953 pre0_mod0_post0

14 Importancia de variables (Random Forest)

Se reportan las variables más relevantes según importancia por impureza.

final_wf <- extract_workflow(final_fit)
rf_fit   <- extract_fit_parsnip(final_wf)

imp <- vip::vi(rf_fit$fit) %>%
  arrange(desc(Importance)) %>%
  slice_head(n = 15) %>%
  mutate(rank = row_number(),
         color_grp = if_else(rank %% 2 == 0, "Naranja", "Azul"))

ggplot(imp, aes(x = Importance, y = reorder(Variable, Importance), fill = color_grp)) +
  geom_col() +
  scale_fill_manual(values = c("Azul" = "#1f77b4", "Naranja" = "#ff7f0e")) +
  labs(
    title = "Importancia de variables - Random Forest",
    subtitle = "Top 15 variables según reducción de impureza (Gini).\nValores mayores implican mayor aporte predictivo.",
    x = "Importancia (impureza)",
    y = "Variable",
    fill = NULL) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold"),
    plot.margin = margin(10, 20, 10, 10),
    legend.position = "none")

Interpretación esperable: variables meteorológicas (viento, humedad, precipitación, temperatura) y ubicación (lat/long) pueden explicar episodios de mayor PM2.5.

15 Modelo comparativo: Regresión Logística:

Se incluye un modelo lineal (baseline) para comparar desempeño.

log_spec <- logistic_reg() %>%
  set_mode("classification")

log_workflow <- workflow() %>%
  add_recipe(aq_recipe) %>%
  add_model(log_spec)

set.seed(7443)
log_rs <- fit_resamples(
  log_workflow,
  resamples = aq_folds,
  control = control_resamples(save_pred = TRUE),
  metrics = metric_set(accuracy, roc_auc, bal_accuracy))

log_metrics <- log_rs %>% collect_metrics()
log_metrics

## # A tibble: 3 × 6
##   .metric      .estimator  mean     n std_err .config        
##   <chr>        <chr>      <dbl> <int>   <dbl> <chr>          
## 1 accuracy     binary     0.780     5 0.00186 pre0_mod0_post0
## 2 bal_accuracy binary     0.697     5 0.00281 pre0_mod0_post0
## 3 roc_auc      binary     0.831     5 0.00165 pre0_mod0_post0

16 Comparación RF vs LR (métricas y ROC superpuesta):

Métricas ROC AUC

bind_rows(
  rf_metrics  %>% mutate(modelo = "Random Forest"),
  log_metrics %>% mutate(modelo = "Regresión Logística")) %>%
  filter(.metric == "roc_auc") %>%
  select(modelo, mean, std_err)

## # A tibble: 2 × 3
##   modelo               mean std_err
##   <chr>               <dbl>   <dbl>
## 1 Random Forest       0.951 0.00103
## 2 Regresión Logística 0.831 0.00165

ROC superpuesta (validación cruzada):

pred_log <- log_rs %>% collect_predictions() %>% mutate(modelo = "Regresión Logística")
pred_rf  <- rf_rs  %>% collect_predictions() %>% mutate(modelo = "Random Forest")

roc_df <- bind_rows(pred_rf, pred_log) %>%
  group_by(modelo) %>%
  roc_curve(pollution_high, .pred_yes, event_level = "second")

autoplot(roc_df) +
  labs(title = "Curva ROC comparativa",
    subtitle = "Random Forest vs Regresión Logística (validación cruzada 5-fold)",
    x = "1 - especificidad",
    y = "sensibilidad")

17 Conclusiones:

El objetivo del trabajo fue modelar la probabilidad de que una medición horaria de PM2.5 en Atenas supere un umbral operativo de 15 µg/m³, utilizando variables meteorológicas, espaciales y temporales.

Desempeño del modelo Random Forest

En validación cruzada (5-fold), el modelo Random Forest obtuvo:

• Accuracy: 0.883
• ROC AUC: 0.951
• Sensibilidad: 0.965
• Especificidad: 0.684
• Balanced Accuracy: 0.824

Estos resultados indican una muy buena capacidad discriminativa (AUC ≈ 0.95), con alta sensibilidad. Es decir, el modelo identifica correctamente la gran mayoría de los episodios de alta contaminación (pollution_high = yes).

La especificidad es menor que la sensibilidad, lo que implica que el modelo tiende a clasificar algunos casos “no contaminados” como contaminados. Sin embargo, dado el contexto ambiental, este sesgo hacia mayor sensibilidad puede considerarse aceptable si el objetivo es minimizar falsos negativos.

En el conjunto de test independiente, los resultados fueron:

• Accuracy: 0.885
• ROC AUC: 0.953
• Sensibilidad: 0.965
• Especificidad: 0.691
• Balanced Accuracy: 0.828

El desempeño en test es consistente con la validación cruzada, lo que sugiere que el modelo generaliza adecuadamente y no presenta señales claras de sobreajuste.

Comparación con Regresión Logística

Como modelo base (baseline), se implementó una Regresión Logística bajo el mismo esquema de validación cruzada.

Resultados:

• Accuracy: 0.780
• ROC AUC: 0.831
• Balanced Accuracy: 0.697

Comparando ambos modelos:

Métrica	Random Forest	Regresión Logística
Accuracy	0.883	0.780
ROC AUC	0.951	0.831
Balanced Accuracy	0.824	0.697

La diferencia en ROC AUC (~0.12 puntos) es sustancial y consistente en todas las métricas evaluadas. Esto sugiere que:

• Las relaciones entre variables meteorológicas, espaciales y la contaminación no son puramente lineales.
• Random Forest captura interacciones y no linealidades que la Regresión Logística no logra modelar adecuadamente.

Por este motivo, se selecciona Random Forest como modelo final.

Interpretación general

Los resultados muestran que:

• Las variables meteorológicas (viento, humedad, temperatura, precipitación) y la localización geográfica contienen información suficiente para predecir episodios de mayor concentración de PM2.5.
• El modelo presenta alta sensibilidad, lo cual es deseable en un contexto ambiental donde subestimar episodios de contaminación puede ser más problemático que sobreestimarlos.
• La inclusión de variables espaciales mejora la capacidad predictiva al capturar heterogeneidad territorial.