Trabajo Final Machine Learning

Trabajo Final Machine Learning

Introducción:

Este análisis forma parte de la evaluación final del módulo Machine Learning. La idea es llevar adelante un análisis de los siniestros en la Ciudad de Buenos Aires utilizando la técnica de Regresión Logística para su clasificación. Este trabajo se encuentra disponible en mi blog personal

Objetivo:

Visualizar espacialmente los siniestros considerados de gravedad ocurridos en la Ciudad Autónoma de Buenos Aires utilizando como fuente los datos abiertos publicados por el municipio.

Analizar y modelar bajo la lógica de clasificación mediante regresión logística la naturaleza de los accidentes de transito y su mortalidad durante el período 2019-2023.

Metodología:

Para llevar adelante los objetivos planteados y, considerando que lo que se busca es encontrar un modelo de clasificación se utilizará el conjunto de técnicas específicas de regresión logística. Esta técnica se utiliza para clasificar una variable respuesta que, en este caso será dicotómica. Para ello, resultará de gran utilidad los resultados y el análisis de las estadísticas del modelo problematizando el grado de mortalidad en los siniestros acontecidos en CABA.

Fuentes de datos:

A los fines de los ideas planteadas, trabajaremos con los datos de siniestros ocurridos durante el período 2019-2023. Según las notas metodológicas, las categorías de la variable “gravedad” varian en:

GRAVE: Toda persona que cuya lesión exige la hospitalización de al menos 24 hs o una atención especializada, como fracturas, conmoción, shock grave y laceraciones importantes. MORTAL: Víctima que fallece dentro de los 30 días de producido el siniestro vial por causas directa o indirectamente atribuibles al hecho. SD: Sin datos sobre la gravedad de las lesiones provocadas. A efectos analíticos, los casos sin datos se corresponden con una alta probabilidad a casos leves.

Procederemos a considerar sólamente los casos declarados como “graves” o “mortales” a efectos de entrenar un modelo que permita clasificarlos.

Instalamos librerias

rm(list = ls()) 
options(scipen = 999)

###################################################################################
# Instalamos librerías de trabajo
if (!require("pacman")) install.packages("pacman")
pacman::p_load("lubridate",
               "tidyverse",
               "leaflet",
               "sf",
               "tidymodels",
               "spatialsample",
               "themis",
               "ggmap",
               "showtext",
               "ggtext")

###################################################################################                                                                                                               

Limpieza y transformación

A los fines de abordar las distintas formas de analizar los datos generaremos dos datasets, uno espacial para mapear los sucesos y otro dataset sin generar geometrías espaciales para construir el modelo. Se advierte que, debido a los filtros aplicados, el dataset espacial cuenta con 1507 casos, mientras que el dataset sin considerar geometrías espaciales contiene 1844 sinisetros.

Levantamos los datos y realizamos una limpieza y recodificación de las variables de interés haciendo hincapié en la variable objetivo que es “gravedad” la cual en una primera vista, tiene datos en mayúsculas y minúsculas en un claro error de ingreso de datos. Tambien procedemos a recodificar la variable hora (hh) según los momentos del día.

siniestros_geo <- read.csv("data/siniestros_viales_hechos.csv", sep = ";") %>% 
  filter(!latitud %in% c("SD", "","sd"),
         !gravedad %in% c("SD","sd")) %>% 
  mutate(longitud = as.numeric(longitud),
         latitud = as.numeric(latitud),
         hh = as.numeric(hh),
         "HoraRec" = case_when(hh >= 0 & hh < 6 ~ "Madrugada", 
                               hh >= 6 & hh < 12 ~ "Mañana",
                               hh >= 12 & hh < 18 ~ "Tarde",
                               TRUE ~ "Noche" ),
         gravedad = case_when(gravedad %in%  c("GRAVE","grave") ~ "Grave",
                                    gravedad == "MORTAL" ~ "Mortal"))%>%
  filter(!is.na(latitud)) %>% 
  st_as_sf(coords =c("longitud", "latitud"), crs = 4326) %>% 
  select(!c(direccion_normalizada,fecha,calle, altura, cruce, otra_direccion, hora, comuna, tipo_de_calle, geocodificacion_caba,tipo_de_dato, participantes))


siniestros <- read.csv("data/siniestros_viales_hechos.csv", sep = ";") %>%
  filter(!gravedad %in% c("SD","sd")) %>% 
  mutate(longitud = as.numeric(longitud),
         latitud = as.numeric(latitud),
         hh = as.numeric(hh),
         "HoraRec" = case_when(hh >= 0 & hh < 6 ~ "Madrugada", 
                               hh >= 6 & hh < 12 ~ "Mañana",
                               hh >= 12 & hh < 18 ~ "Tarde",
                               TRUE ~ "Noche" ),
         gravedad = case_when(gravedad %in%  c("GRAVE","grave") ~ "Grave",
                                    gravedad == "MORTAL" ~ "Mortal"))%>%
  select(!c(direccion_normalizada,fecha,calle, altura, cruce, otra_direccion, hora, comuna, tipo_de_calle, geocodificacion_caba,tipo_de_dato, participantes))


comunas <- st_read("data/comunas.geoson", 
                  crs = 4326)

## Reading layer `comunas' from data source 
##   `C:\Users\Guille\Documents\Proyectos\FLACSO\Machine Learning\tp final\data\comunas.geoson' 
##   using driver `GeoJSON'
## Simple feature collection with 15 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53152 ymin: -34.70529 xmax: -58.33516 ymax: -34.52649
## Geodetic CRS:  WGS 84

tabla_1 <- siniestros %>% 
  group_by(aaaa, gravedad) %>% 
  summarise(tot = n()) %>% 
  mutate(porcentaje = round(tot/sum(tot),4)) %>% 
  mutate(Total_Acum = cumsum(tot)) %>% 
  mutate(Porcentaje_Acum = cumsum(porcentaje))

ggplot(tabla_1, aes(x = aaaa, y = tot, fill = gravedad))+
  geom_col(position = "dodge")+
  geom_text(aes(label = percent(porcentaje)), vjust = -0.5, check_overlap = TRUE, position = position_dodge(width = .9), size = 3)+
  scale_fill_manual(values = c("#f55951","#000000"))+
  theme_minimal()+
  labs(
    title = "<span style = 'color:#272643;'> Evolución anual de siniestros: <span style = 'color:#f55951;'> Graves <span style = 'color:#543c52;'> - <span style = 'color:#000000;'> Mortales",
    subtitle = "",
    caption = "<b>Fuente</b>: ReNaBaP",
    x = "Año",
    y = "Barrios Populares"
  )+
  theme(
    legend.position = "none",
    plot.title = element_markdown(face = "bold", size = 15),
    plot.caption = element_markdown(size = 12),
  )

ggplot()+
  geom_sf(data=comunas)+
  geom_sf(data=siniestros_geo, aes(color=gravedad), alpha=0.75, show.legend = TRUE)+
  scale_color_manual(values = c("#f55951","#000000"))+
  theme_minimal()+
  theme_void()+
  labs(
    title = "<span style = 'color:#272643;'> Siniestros viales - CABA 2019-2023",
    caption = "<b>Fuente</b>: Datos Abiertos CABA",
    y = "",
    x = ""
  )+
  theme(
    plot.title = element_markdown(size = 20, hjust = 0.5),
    plot.caption= element_markdown(size = 12),
    legend.position = "top",
    legend.title = element_blank(),
    legend.box.margin = margin(11, 16, 6, 6)
  )+
  facet_wrap(~aaaa, nrow = 1)+
  guides(fill=guide_legend(nrow=1,byrow=TRUE))

Al aplicar filtros vinculados a eliminar registros sin georeferenciación, el año 2019 sólo contiene casos mortales. Los años siguientes contienen una mayor cantidad de casos de gravedad manteniendo ciertos patrones para los siniestros fatales.

A continuación, se muestra un mapa interactivo para visualizar la totalidad de casos registrados en el dataset con geometrías espaciales.

factpal <- colorFactor(palette = c("#f55951","#000000"), 
               levels = siniestros_geo$gravedad)

oceanIcons <- iconList(
  ship = makeIcon(
    "https://upload.wikimedia.org/wikipedia/commons/thumb/9/9a/Medical_hospital_emergency_beds.svg/240px-Medical_hospital_emergency_beds.svg.png",
    "https://upload.wikimedia.org/wikipedia/commons/thumb/9/9a/Medical_hospital_emergency_beds.svg/240px-Medical_hospital_emergency_beds.svg.png",
    18,
    18
  ),
  pirate = makeIcon(
    "https://upload.wikimedia.org/wikipedia/commons/thumb/c/c3/Maki2-danger-24.svg/240px-Maki2-danger-24.svg.png",
    "https://upload.wikimedia.org/wikipedia/commons/thumb/c/c3/Maki2-danger-24.svg/24px-Maki2-danger-24.svg.png",
    24,
    24
  )
)


 mapa_interactivo <- leaflet() %>%
     addTiles() %>%
     addMarkers(data = siniestros_geo, icon = ~oceanIcons, clusterOptions = markerClusterOptions(),popup = ~ glue::glue("<b>Estado</b>:{gravedad}<br/>
                                                                    <b>Año:</b>{aaaa}<br/>
                                                                    <b>Cantidad de Víctimas:</b>{n_victimas}<br/>
                                                                    <b>Hora:</b>{hh}<br/>"))

mapa_interactivo

Modelado

En primer lugar procedemos a la confección de set de entrenamientos y testeos considerando 0.8 como proporción y tomando el dataseet que no contiene coordenadas geográficas.

#seteo la semilla
set.seed(123)
#divido los datos
sini_split <- initial_split(siniestros, prop = 0.80, strata = gravedad) #estratificamos segun las clases 
sini_train <- training(sini_split)
sini_test <- testing(sini_split)

Folds con validación cruzada

set.seed(234)
sini_folds <- vfold_cv(sini_train, strata = gravedad)
sini_folds

sini_train

sini_test 

sini_train %>%
  count(gravedad)

sini_test %>%
  count(gravedad)

Se verifica que las proporciones correspondientes a la distribución de las categorías de la variable “gravedad” se conservan.

Receta

El siguiente paso es determinar la receta del modelo tomano “gravedad” como la variable a calsificar.

logistic_recipe <-
  recipe(formula = gravedad ~ ., data = sini_train) %>%
  update_role(id_hecho, new_role = "id") %>%
  step_unknown(all_nominal_predictors()) %>%
  step_other(all_nominal_predictors(), threshold = 0.03) %>%
  step_downsample(gravedad)

Modelización

Especificamos el modelo, en este caso, una regresión logística.

logistic_spec <- logistic_reg()

Workflow

Agregamos al Workflow la receta y la especificación del modelo.

log_workflow <- workflow() %>%
  add_recipe(logistic_recipe) %>%
  add_model(logistic_spec)

Entrenamiento del algoritmo

Con la función fit_resamples() procedemos a hacer el ajuste correspondiente.

doParallel::registerDoParallel() 
set.seed(7443)
logistic_rs <-
  fit_resamples(log_workflow,
    resamples = sini_folds,
    control = control_resamples(save_pred = TRUE)
  )

Performance

logistic_rs %>%
  collect_metrics()

La performance muestra resultados aceptables en el set de entrenamiento. La precisión casi alcanza el 80%, mientras que el área bajo la curva ROC llega al valor de 0.88. El brier_class arroja valores bajos, lo cual es bueno para nuestro modelo

collect_predictions(logistic_rs) %>%
  group_by(id) %>%
  roc_curve(gravedad, .pred_Grave) %>%
  autoplot()

Predicción

En esta etapa con la función last_fit vamos a seleccionar el mejor modelo para predecir en el set de datos de TEST.

final_fitted <- last_fit(log_workflow, sini_split)
final_fitted 

Métricas

final_fitted %>%
  collect_metrics()

Las métricas del modelo final arrojan cifras aceptables para la precisión (78 %) y valores buenos para el área bajo la curva ROC. Tambien muestra valores bajos para brier_class

Predicciones del modelo

final_fitted %>%
  collect_predictions()

Matriz de confusión

Para finalizar y completar la consigna, vamos a crear una matriz de confusión con las predicciones para contrastarla con los casos reales.

final_fitted %>%
  collect_predictions() %>%
  conf_mat(gravedad, .pred_class) 

##           Truth
## Prediction Grave Mortal
##     Grave    155     16
##     Mortal    47     66

Como podemos ver la diagonal principal contiene la mayoría de los casos, dejando en valores bajos los falsos positivos y los falsos negativos.

Conclusión

A partir del empleo de técnicas de modelado, pudimos arribar una problemática severa de la movilidad urbana como lo son los siniestros. Con datos abiertos y técnicas de visualización y regresión logística construimos un modelo que, a pesar de la escasa cantidad de datos, obtuvo cifras de precisión considerables.