Rafael Zambrano: Predicción de delitos en CABA

Rafael Zambrano

⚠️ Antes de empezar …

Aunque tengo la fortuna de trabajar en el área de ciencia de datos, no soy un experto en el tema, estoy aprendiendo constantemente 📖, ¿El código mostrado puede escribirse de manera más limpia y eficiente? ¿Se pudo aplicar un mejor enfoque? ¿La función de pérdida, las métricas de evaluación u optimización seleccionadas podrían ser otras?.

A todo lo anterior, Sí, muy probablemente Sí. Siéntase libre de contribuir comentando o sugiriendo cualquier idea que considere ✌️. Exponer nuestro trabajo, intercambiar y compartir conocimiento es la forma más linda de aprender 📚, este es el espíritu de todas las publicaciones de este Blog Posts 💪🏻.

🔴 Aunque esta publicación tiene un origen académico, especificamente como TIF del postgrado de Métodos Cuantitativos para la Gestión y Análisis de Datos en Organizaciones de la Facultad de Económicas de la Universidad de Buenos Aires. Se presenta de una manera más distendida con el objetivo de hacer más agradable su lectura, para ello los fragmentos de código extensos tienen la opción de plegarse, y 👉🏻 se acompañan las descripciones con emojis y memes 🤪😚👻.

Resumen 📖

🔹 La recopilación de datos referentes a la ocurrencia de delitos según una fecha y localización determinada, brinda un amplio abanico de oportunidades para la explotación de información en diversos tipos de organizaciones. El presente trabajo tiene como objetivo 👉🏻👉🏻👉🏻 implementar modelos predictivos de aprendizaje automático, que permita relacionar una ubicación 🌍 y un momento determinado 📆 con la ocurrencia de delitos 🚨, con base en los casos registrados entre los años 2017 y 2019 en la Ciudad Autónoma de Buenos Aires (CABA). Se plantea el estudio dentro de una organización hipotética que pueda usar este proyecto para apoyar los procesos de toma de decisiones, en este caso se plantea como institución La Policía de la Ciudad de Buenos Aires.¹

🔹 Específicamente se busca implementar modelos que permitan predecir la ocurrencia de delitos en un conjunto de esquinas 📍 de CABA, agregando elementos meteorológicos 🌦️ y del entorno cercano a la ocurrencia del crimen. En base a los delitos registrados, se busca predecir la ocurrencia para el mes de diciembre del año 2019, se aplica un enfoque de regresión, clasificación multiclase y clasificación binaria múltiple. Se utilizará el modelo de series de tiempo Croston y los modelos de machine learning XGBoost y LightGBM, a través del framework de modelado tidymodels (Kuhn and Wickham 2020b) .

Palabras clave: Ocurrencia de delitos, predicción, modelos predictivos, entorno-tiempo, toma de decisiones.

Principales bibliotecas 📚

Show code

library(skimr)
library(tidyverse)
library(rgdal)
library(leaflet)
library(geosphere)
library(sp)
library(readxl)
library(ggsci)
library(ggforce)
library(tidymodels)
library(rcrimeanalysis)
library(psych)
library(plyr)
library(DataExplorer)
library(DescTools)
library(ggpubr)
library(timetk)
library(parallel)
library(janitor)
library(spatialEco)
library(osmdata)
library(lubridate)
library(patchwork)
library(tsintermittent)
library(treesnip)
library(rangeBuilder)
library(themis)
library(zoo)
library(hrbrthemes)
library(cowplot)
library(echarts4r)
library(erer)
library(gganimate)

Paquete 📦 especifico para este proyecto

🔎Dentro del flujo de trabajo implementado, surgió la necesidad de crear un conjunto de funciones referentes al cálculo de distancias entre puntos geolocalizados, y a la evaluación de predicciones. Además de la implementación de técnicas de “ventanas deslizantes” (sliding window) , con adaptaciones específicas de este estudio. Para esto fue creado el paquete 📦 sknifedatar² , para consultar más detalles visite el website , el mismo puede descargarse e instalarse de la siguiente manera:

#devtools::install_github("rafzamb/sknifedatar")
library(sknifedatar)

Importación de datos 🔌

Todos los datos se obtuvieron del repositorio público del Gobierno de la Ciudad Autonoma de Buenos Aires ³ . A modo ilustrativo, se mencionan algunos de los elementos de entorno relevantes usados en el estudio, por ejemplo, estaciones de metro 🚇, comisarías 👮🏻, clínicas 🏥, hoteles 🏨, universidades 🏫, locales nocturnos 🍻, entre otros. Los nombres son bastante descriptivos, por lo que es fácil asociarlos con el elemento que representan. Se analizan 34 elementos de entorno, el listado total puede verificarse en los dos siguientes fragmentos de código.

❗A continuación, se procede a realizar la extracción de los conjuntos de datos desde dicho repositorio, en caso de querer reproducir este estudio se debe tomar la precaución de verificar posibles cambios en las rutas de las fuentes de datos.

Datos tabulares 📓

Corresponde a datos que al ser descargados del repositorio tienen un formato “.csv” y una estructura tabular.

Show code

delitos_2019 <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/mapa-del-delito/delitos_2019.csv")

delitos_2018 <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/mapa-del-delito/delitos_2018.csv")

delitos_2017 <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/mapa-del-delito/delitos_2017.csv")

universidades <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/universidades/universidades.csv")

local_bailables <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/locale
                            s-bailables/locales-bailables.csv") 

centros_de_salud_privados <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/dataset
                                      s/centros-de-salud-privados/centros-de-salud-privados.csv")

hospitales <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/salud/hospitales/hospitales.csv")

comisarias <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/comisarias-polici
                       a-de-la-ciudad/comisarias-policia-de-la-ciudad.csv")

embajadas <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/embajadas/embajadas.csv")

bocas_de_subte <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/bocas-de-subt
                           e/bocas-de-subte.csv")

atm <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/cajeros-automatico
                s/cajeros-automaticos.csv")

bancos <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/bancos/accesibilidad-en-bancos.csv")

gasolina <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/estaciones-d
                     e-servicio/estaciones_servicio_caba.csv")

gastronomica<- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/ofert
                        a-gastronomica/oferta_gastronomica.csv") 

hotel<- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/alojamientos-turistico
                 s/alojamientos-turisticos.csv") %>% 
  drop_na(long,lat)

metrobus <-  read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/transporte/metrobu
                     s/estaciones-de-metrobus.csv") 

taxi <-  read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/paradas-de-taxis/paradas-de-taxis.csv")

carteles <-  read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/carteles-luminoso
                     s/carteles-luminosos.csv")

colectivo <-  read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/colectivos/paradas-de-colectivo.csv") 

farmacias <-  read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/farmacias/farmacias.csv")

centros_salud_comunitarios <-  read_csv("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/salud/centros-salud
                                       -accion-comunitaria-cesac/centros_de_salud_nivel_1_BADATA_WGS84.csv") 

centros_medicos_barriales <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/centros-medicos-barriales
                                      /centros-medicos-barriales.csv")

estadios <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/estadios/estadios.csv")

skate <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/pistas-de-skate/pistas-de-skate.csv")

gimnasios <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/gimnasios/gimnasios.csv")

wifi <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/puntos-de-wi-fi-publicos/sitios-de-wifi.csv")

garajes <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/garajes-comerciales/garajes-comerciales.csv") 

educativos <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/establecimientos-educativos/
                      establecimientos-educativos.csv")

bibliotecas <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/bibliotecas/bibliotecas.csv")

basilicas <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/basilicas/basilicas.csv")

espacios_culturales <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/espacios-culturales/
                                espacios-culturales.csv")

organizaciones_sociales <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/organizaciones-sociales/
                                    organizaciones-sociales.csv")

estaciones_ferrocarril <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/estaciones-de-ferrocarril/
                                   estaciones-de-ferrocarril.csv")

temperatura_mes_histoico <- read_csv("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/direccion-general-de
                                     -estadisticas-y-censos/registro-temperatura-ciudad/historico_temperaturas.csv") 

lluvia_mes_histoico <- read_csv("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/direccion-general-de
                                -estadisticas-y-censos/registro-precipitaciones-ciudad/historico_precipitaciones.csv") 

viento_mes_histoico <- read_csv("https://data.buenosaires.gob.ar/dataset/13df72f9-0080-4480-ad89-68baa02e38b6/resource
                                /6e20f472-3e39-4a5c-93b2-9d5d58cdfc04/download/historico_velocidad_maxima_viento.csv", 
                                skip = 1) 

hogares_paradores <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/hogares-y-paradores/hogares-
                              y-paradores.csv") 

bus_turistico <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/paradas-y-recorridos-bus-turistico/
                          paradas_bus_turistico.csv") 

consulados <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/consulados/consulados.csv") 

bomberos <- read_csv("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/cuarteles-y-destacamentos-de-
                     bomberos/cuarteles-y-destacamentos-de-bomberos-de-policia-federal-argentina.csv")

Datos Geoespaciales 🌎

Estos datos tiene un proceso de lectura no tan directo como los anteriores, se debe crear una carpeta previamente para alojar el archivo “.zip,” luego descargarlo y descomprimirlo. Posteriormente, leer el objeto con la función readOGR y transformarlo mediante la función spTransfor ,estableciendo una proyección espacial.

Show code

proyeccion <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

# Barrios Caba
dir.create(file.path(getwd(), "shape_barrios"))  

download.file("http://cdn.buenosaires.gob.ar/datosabiertos/datasets/barrios/barrios-zip.zip",
              destfile="shape_barrios/barrios-zip.zip")

unzip("shape_barrios/barrios-zip.zip",exdir= "shape_barrios")

shape_barrios <- readOGR("./shape_barrios", layer = "barrios_badata") #Limites de Barrios CABA

shape_barrios_caba <- spTransform(shape_barrios, CRS(proyeccion))

# barrios-populares
dir.create(file.path(getwd(), "shape_barrios_populares"))  

download.file("https://cdn.buenosaires.gob.ar/datosabiertos/datasets/desarrollo-humano-y-habitat/barrios-populares
              /barrios-populares-badata.zip", destfile<-"shape_barrios_populares/barrios-populares-badata.zip")

unzip("shape_barrios_populares/barrios-populares-badata.zip",exdir= "shape_barrios_populares")

shape_barrios_populares <- readOGR("./shape_barrios_populares", layer = "barrios_vulnerables") 

# distritos economicos
shape_distritos_economicos <- readOGR("./distritos-economicos-zip",layer = "distritos_economicos")

# inclusion_urbana
shape_inclusion_urbana <- readOGR("./unidades-territoriales-de-inclusion-urbana-zip",
                           layer = "unidades-territoriales-de-inclusion-urbana")

Visualización de objetos espaciales 🗺️

A modo de ejemplo se visualizan algunos objetos considerados como elementos del entorno cercano a la ocurrencia del delito, mediante el paquete 📦 leaflet⁴ .

barrios_populares <- spTransform(shape_barrios_populares, CRS(proyeccion))

leaflet(barrios_populares) %>%
  addTiles() %>% 
  addPolygons()

Figure 1: Barrios populares de CABA. Fuente: (OpenStreetMap contributors 2017)

inclusion_urbana <- spTransform(shape_inclusion_urbana, CRS(proyeccion))

leaflet(inclusion_urbana) %>%
  addTiles() %>% 
  addPolygons()

Figure 2: Zonas de inclusión urbana de CABA. Fuente: (OpenStreetMap contributors 2017)

Limpieza de datos 🧹

🔹Existen problemas de consistencia relacionados con el nombre de los campos, nombrando los atributos latitud y longitud de diversas formas en distintos datasets, inclusive en algunas situaciones con formatos diferentes. Para estos dos atributos, se estandariza el nombre y formato en todos los conjuntos de datos. Serán excluidos los elementos de entorno cercano al delito que presenten valores faltantes en la geolocalización. A continuación, se muestra este proceso inicial de limpieza.

Show code

delitos_2019 <- delitos_2019 %>% drop_na(long,lat)

delitos_2018 <- delitos_2018 %>% drop_na(long,lat)

delitos_2017 <- delitos_2017 %>% drop_na(long,lat)

universidades <- universidades %>% distinct(universida,long,lat)

local_bailables <- local_bailables %>%
  rename(long = longitud,lat = latitud) %>% 
  drop_na(long,lat)

wifi <- wifi %>% 
  filter(!tipo %in% c("Biblioteca","Deporte","Cultura","Hospitales y Centros de Salud",
                      "Metrobus","Museos","Salud"))

hotel_baja <- hotel %>% filter(categoria %in% c("1*","2*","3*"))

hotel_alta <- hotel %>% filter(categoria %in% c("4*","5*"))

cine <- espacios_culturales %>% 
  filter(funcion_principal == "SALA DE CINE") %>% 
  rename(long = longitud,lat = latitud)

teatros <- espacios_culturales %>% 
  filter(subcategoria == "SALA DE TEATRO") %>% 
  rename(long = longitud,lat = latitud)

hogares_paradores <- hogares_paradores %>% 
  rename(long = X, lat = Y) %>% 
  drop_na(long,lat)

gimnasios <- gimnasios %>% drop_na(long,lat)

hospitales <- hospitales %>% 
  mutate(coordenadas = str_trim(gsub("\\(|\\)|[[:alpha:]]","",WKT))) %>% 
  select(coordenadas) %>% 
  separate(coordenadas,c("long","lat"), sep = " ") %>% 
  mutate_all(as.double)

metrobus <- metrobus %>% 
  mutate(coordenadas = str_trim(gsub("\\(|\\)|[[:alpha:]]","",WKT))) %>% 
  select(coordenadas) %>% 
  separate(coordenadas,c("long","lat"), sep = " ") %>% 
  mutate_all(as.double)

centros_salud_comunitarios <- centros_salud_comunitarios %>% 
  mutate(coordenadas = str_trim(gsub("\\(|\\)|[[:alpha:]]","",WKT))) %>% 
  select(coordenadas) %>% 
  separate(coordenadas,c("long","lat"), sep = " ") %>% 
  mutate_all(as.double)

colectivo <- colectivo %>% rename(long = stop_lon, lat = stop_lat)

barrios_populares_centroides <- coordinates(barrios_populares) %>% 
  as.data.frame() %>% 
  `colnames<-`(c("long","lat"))

Construcción del Dataset 🛠️

Para este estudio se considerarán los delitos ocurridos en los años 2017, 2018 y 2019, los tres años suman en total 354.220 delitos, a continuación se unifican los data frames de delitos y se analiza su estructura.

delitos <- bind_rows(delitos_2019,
                    delitos_2018,
                    delitos_2017 %>% mutate(franja_horaria = as.double(franja_horaria)))

skim(delitos)

Table 1: Data summary
Name	delitos
Number of rows	354220
Number of columns	10
_______________________
Column type frequency:
character	3
Date	1
numeric	6
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	n_unique
tipo_delito	0	1.00	8	21	4
subtipo_delito	306758	0.13	6	15	4
barrio	0	1.00	4	17	48

Variable type: Date

skim_variable	n_missing	complete_rate	min	max	median	n_unique
fecha	0	1	2017-01-01	2019-12-31	2018-07-05	1095

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
id	0	1	307893.60	105197.95	125088.00	216197.75	309027.50	399796.25	488541.00	▇▇▇▇▇
franja_horaria	40	1	13.78	6.33	0.00	10.00	15.00	19.00	23.00	▃▅▅▇▇
cantidad_registrada	0	1	1.00	0.01	1.00	1.00	1.00	1.00	4.00	▇▁▁▁▁
comuna	0	1	7.30	4.65	1.00	3.00	7.00	12.00	15.00	▇▆▅▅▆
lat	0	1	-34.61	0.03	-34.70	-34.63	-34.61	-34.59	-34.53	▁▃▇▅▁
long	0	1	-58.44	0.04	-58.53	-58.47	-58.43	-58.40	-58.34	▃▆▇▇▃

Tipos de delitos 🚨

👉🏻 No se observan valores faltantes en los atributos de geolocalización y fecha, sí en las franjas horarias pero para este estudio no serán consideradas. La columna "cantidad_registrada" refleja el número de crímenes ocurridos en una localización (lat,long) y franja horaria determinada, para la mayoría de los casos el valor es 1, las pocas observaciones que reflejan un valor superior serán tratados como un delito único.

Show code

tabla <- as.data.frame.matrix(table(delitos$tipo_delito,delitos$subtipo_delito, exclude = F))
names(tabla)[5] = "NA"

knitr::kable(tabla %>% rownames_to_column("Tipos / Subtipos de Delitos"),
             caption = "Tipos y subtipos de delitos")

Table 2: Tipos y subtipos de delitos
Tipos / Subtipos de Delitos	Doloso	Hurto Automotor	Robo Automotor	Siniestro Vial	NA
Homicidio	351	0	0	363	0
Hurto (sin violencia)	0	16320	0	0	115834
Lesiones	0	0	0	22823	0
Robo (con violencia)	0	0	7605	0	190924

🔎 Para este estudio serán excluidos los crímenes relacionados con siniestros viales y lesiones 🚗, esto bajo el spuesto de asumir que estos crímenes responden a patrones muy distintos del resto de los delitos .

Show code

delitos %>%
  dplyr::filter(subtipo_delito != "Siniestro Vial" | is.na(subtipo_delito)) %>%
  dplyr::count(tipo_delito) %>% 
  mutate(prop = round(n/sum(n)*100,2)) %>% 
  ggplot(aes(x = tipo_delito, y = prop, color = tipo_delito)) +
  scale_color_jco() +
  geom_segment(aes(xend = tipo_delito, yend = 0), show.legend = F) +
  geom_point(aes(size = prop), show.legend = F) +
  geom_label(aes(label = paste0(prop,"% / n = ",n)),
             fill = "white", 
             hjust = "inward",
             show.legend = F) +
  labs(y = "%",
       x = "Tipo de delitos") +
  coord_flip() +
  theme_minimal()

Figure 3: Tipos y subtipos de delitos

❗En la Figura 3 , se observa la distribución de los delitos definitivos sobre los cuales se modelará, el enfoque de este estudio no reside en discriminar entre los tipos de delitos, todos los delitos se contabilizarán como un crimen ocurrido. En el siguiente fragmento de código, se realizan las exclusiones comentada en esta sección.

delitos <- delitos %>% 
  filter(subtipo_delito != "Siniestro Vial" | is.na(subtipo_delito)) %>% 
  uncount(cantidad_registrada)

Análisis exploratorio de datos 🔭

❗Se presenta un análisis descriptivo de la ocurrencia histórica de delitos, su distribución geográfica y la evolución de los factores meteorológicos. Para efectos de legibilidad de esta publicación, se excluye el análisis de los elementos de entorno del delito, sin embargo, también es necesario visualizar y analizar estos factores.

Ocurrencia histórica de delitos 📆

Show code

ggarrange(delitos %>%
            group_by(fecha) %>% 
            tally() %>% 
            ungroup() %>% 
            plot_time_series(fecha,n, .interactive = F, .title = "Frecuencia diaria"),
          delitos %>% 
            group_by(fecha) %>% 
            tally() %>% 
            summarise_by_time(
              fecha, .by = "month",
              casos = sum(n)
              ) %>%
            plot_time_series(
              fecha,
              casos, 
              .interactive = F,
              .color_lab = "Month", 
              .title = "Frecuencia mensual"), 
          ncol = 1) %>% 
  annotate_figure(top = text_grob("Delitos Ocurridos Dic 2017 - Dic 2019", color = "black", face = "bold", size = 14))

Histórico de delitos ocurridos en CABA. Elaborado a través del paquete timetk [@timetk]

Figure 4: Histórico de delitos ocurridos en CABA. Elaborado a través del paquete timetk (Dancho and Vaughan 2020)

❗En la Figura 4, se observa un patrón estacional en la ocurrencia de delitos, presentando saltos y caídas en los mismos meses durante los años estudiados, este estudio se limitará temporalmente al periodo diciembre 2017 - diciembre 2019.

delitos <- delitos %>% 
  filter(fecha > "2017-11-30")

Distribución geográfica 🌐

delitos %>% 
  dplyr::rename(latitude = lat, longitude = long) %>% 
  kde_map(pts = F)

Figure 5: Distribución geográfica de los delitos ocurridos. Elaboradoa través del paquete rcrimeanalysis (Spaulding and Morris 2020)

👀En la Figura 5, se muestra geográficamente una mayor concentración de delitos en la zona este y central de la ciudad, las zonas que no están marcadas en color rojo también registran delitos, pero con una menor frecuencia. Esta visualización fue realizada a través del paquete 📦 rcrimeanalysis⁵.

Show code

delitos %>% 
  dplyr::count(fecha) %>% 
  filter(year(fecha) != 2017) %>% 
  dplyr::mutate(year = format(fecha, "%Y")) %>% 
  group_by(year) %>% 
  e_charts(fecha) %>% 
  e_calendar(range = "2018") %>% 
  e_calendar(range = "2019",top="350") %>% 
  e_heatmap(n, coord_system = "calendar") %>% 
  e_visual_map(max = 470, min = 130, bottom = 140) %>% 
  e_title("Calendario de delitos: Mes /Día de la semana", "Número de delitos ocurridos") %>% 
  e_tooltip("item")

Figure 6: Deltios ocurridos por días de la semana y meses del año. Elaborado mediante el paquete echarts4r (Coene 2020)

🚔 La figura 6 , muestra el historial de delitos ocurridos por días de la semana y meses del año. El calendario es interactivo ✨, la barra de ocurrencia de la izquierda permite que se ajuste interactivamente para verificar patrones en los datos. Los fines de semana se registra la incidencia más baja de delictividad, los meses de junio, julio y agosto presentan un bajo número de días que registren más de 300 crímenes, siendo este un patrón común para el 2018 y 2019.

delitos_def <- delitos %>% 
  mutate(fecha_split = as.yearmon(fecha, "%Y-%m")) %>% 
  select(id,  lat, long, fecha_split)

Variables climáticas 🌦️

El número de días de lluvia 🌧️, la cantidad de milímetros de agua ☔, la temperatura 🌡️ y velocidad del viento 🌬️ promedio mensual, corresponden a los elementos meteorológicos que se estudiarán.

temperatura_mes_histoico <- temperatura_mes_histoico %>% 
  clean_names() %>% 
  bind_cols(meses = match(.$mes,meses)) %>% 
  mutate(fecha = paste0(ano,"-",meses)) %>%
  mutate(fecha = as.yearmon(fecha, "%Y-%m")) %>% 
  select(fecha, maxima,minima,media)

lluvia_mes_histoico <- lluvia_mes_histoico %>% 
  clean_names() %>% 
  bind_cols(meses = match(.$mes,meses)) %>% 
  mutate(fecha = paste0(ano,"-",meses)) %>%
  mutate(fecha = as.yearmon(fecha, "%Y-%m")) %>% 
  select(-c(ano,meses,mes))

viento_mes_histoico <- viento_mes_histoico %>% 
  clean_names() %>% 
  select(-3) %>% 
  setNames((c("ano","mes","veloc_viento"))) %>% 
  bind_cols(meses = match(.$mes,meses)) %>% 
  mutate(fecha = paste0(ano,"-",meses)) %>%
  mutate(fecha = as.yearmon(fecha, "%Y-%m")) %>% 
  select(-c(ano,meses,mes))

Show code

p1 <- lluvia_mes_histoico %>%
  filter(year(fecha) > 2016) %>% 
  pivot_longer(!fecha, names_to = "Date", values_to = "valores") %>% 
  group_by(Date) %>% 
  plot_time_series(fecha, valores, .interactive = F, 
                   .title = "LLuvia: mm agua/días de lluvia")

p2 <- ggarrange(
  viento_mes_histoico %>% 
    filter(year(fecha) > 2016) %>% 
    plot_time_series(fecha, veloc_viento, .interactive = F, .title = "Viento"),
  temperatura_mes_histoico %>% 
    filter(year(fecha) > 2016) %>% 
    plot_time_series(fecha, media, .interactive = F, .title = "Temperatura"), 
  ncol = 1)

ggarrange(p1,p2,ncol = 2) %>% 
  annotate_figure(top = text_grob("Dic 2017 - Dic 2019", color = "black",
                                  face = "bold", size = 14))

Figure 7: Histórico de factores meteorológicos

❗Los factores meteorológicos también requieren una transformación en los datos para poder incorporarlos en el estudio, dicho proceso se muestra en el primer fragmento de código. Se observan patrones estacionales propios de la natruraleza de estos los factores, aunque también puntos extremos en algunas de las series.

Selección de esquinas de CABA 📍

👉🏻 Las esquinas de la ciudad cumplen la función de permitir la construcción del factor “entorno-tiempo” 📐/⏰, que busca relacionarse con la incidencia de la criminalidad, mediante la contabilización del número de delitos que ocurren en las cercanías de las esquinas en distintos espacios de tiempo. Para esto, a través del paquete 📦 osmdata⁶ y del repositorio OpenStreetMap (OpenStreetMap contributors 2017), se extraen las calles y avenidas de la ciudad, a continuación se muestra el procedimiento.

CABA_OS <- getbb("Ciudad Autonoma de Buenos Aires")

CABA_OS_poly <- getbb("Ciudad Autonoma de Buenos Aires", format_out = "sf_polygon")

CABA <- opq(CABA_OS) %>% 
  add_osm_feature(key = "highway") %>% 
  osmdata_sf()

ggplot() +
  geom_sf(data = CABA$osm_lines) +
  theme_minimal() +
  ggtitle("Principales calles y avenidas de CABA")

Figure 8: Principales calles y avenidas de CABA

Intercepción de calles y avenidas 🛣️

📌 El único dato de interés de las avenidas y calles consiste en sus intersecciones 🚦 , los puntos de intersección resultantes hacen referencias a las esquinas, de estas esquinas se extrae la longitud y latitud, siendo estos los campos claves para poder incorporar esta información al estudio 💡. A continuación, se muestra el procedimiento.

calles <- CABA$osm_lines

calles <- st_intersection(calles, CABA_OS_poly)

intercepcion_calles <- st_cast(calles, "POINT") %>% 
  st_coordinates() %>% 
  as.data.frame() %>% 
  distinct() %>% 
  rename(lat = Y, long = X) %>% 
  rownames_to_column("id")

intercepcion_calles_1 <- intercepcion_calles

coordinates(intercepcion_calles_1) <- c('long', 'lat')
proj4string(intercepcion_calles_1) <- CRS(proyeccion)
intercepcion_calles_1 <- spTransform(intercepcion_calles_1, CRS(proyeccion))

Selección óptima de esquinas en función a la proximidad ⚛️

📌 Se observa que muchas esquinas están muy próximas entre sí, en este caso, al realizar el conteo de delitos dentro de un radio de metros determinado, se podría llegar a contabilizar un delito en más de una esquina 🚫. Por lo tanto, se propone seleccionar las esquinas que estén a una distancia de más de 200 metros entre sí, para esto se implemneta un algoritmo de filtro de proximidad desarrollado en el paquete 📦 rangeBuilder⁷ . A continuación, se implementa el algoritmo de optimización.

intercepcion_calles_2 <- filterByProximity(intercepcion_calles_1,dist=0.20, mapUnits=F)

Show code

ggarrange(
data.frame(x = coordinates(intercepcion_calles_1)[,1],
           y = coordinates(intercepcion_calles_1)[,2]) %>% 
  ggplot(aes(x = x, y = y)) +
  geom_point(alpha = 0.5, size = 0.5) +
  theme_minimal() +
  ggtitle("Esquinas Originales") +
  xlab("") +
  ylab(""),
data.frame(x = coordinates(intercepcion_calles_2)[,1],
           y = coordinates(intercepcion_calles_2)[,2]) %>% 
  ggplot(aes(x = x, y = y)) +
  geom_point(alpha = 0.5, size = 0.5) +
  theme_minimal() +
  ggtitle("Esquinas Óptimas") +
  xlab("") +
  ylab(""), 
ncol = 2) %>% 
  annotate_figure(top = text_grob("Intercepciones de calles de CABA", color = "black", face = "bold", size = 14))

Figure 9: Esquinas de CABA antes y después de la aplicación del algoritmo

❗En la Figura 9, se muestran las esquinas sobre el mapa de CABA, comparando las intercepciones originales y las resultantes luego de la aplicación del algoritmo de selección. Originalmente se habían encontrado 22.048 esquinas, luego de la aplicación del algoritmo se seleccionaron 2.417 esquinas. Las intercepciones definitivas son almacenadas como un data frame, extrayendo únicamente las coordenadas.

intercepcion_calles_def <- intercepcion_calles_2 %>% 
  coordinates %>% 
  as.data.frame() %>% 
  rownames_to_column("id")

Cálculo mensual de los delitos ocurridos en cada esquina de CABA 🚔

Función de distancia entre delitos y esquinas

🔎 El siguiente paso consiste en contabilizar cuántos delitos ocurrieron en las cercanías de cada esquina, para ello, se fija una distancia de 150 metros y se implementa un algoritmo almacenado en la función pertenencia_punto. Sin embargo, se aclara que debido al gran número de delitos y esquinas a analizar, puede demorar varias horas en ejecutarse ⏳, se propone realizar una implementación por lotes, como se muestra en el siguiente ejemplo:

pertenencia_punto <- function(data, referencia, metros){

  crime <- data %>%
    dplyr::select(.data$long,.data$lat) %>%
    split(1:nrow(data))

  esquinas <- referencia[,c("long","lat")]

  vector_esquinas<- c()

  Esquina <- parallel::mclapply(crime, function(x){

    for (i in 1:nrow(esquinas)) {

      vector_localizacion <- geosphere::distm(x,  esquinas[i,], fun = geosphere::distHaversine)

      vector_esquinas[i] <- vector_localizacion
    }

    vector <- which(vector_esquinas <= metros)

    n_total <- ifelse(length(vector) == 0, 0, vector)
  })
  return(Esquina)
}

esquinas_lote_1 <- pertenencia_punto(data = delitos_def[1:25000,],
                                     referencia = intercepcion_calles_def,
                                     metros = 150)

👉🏻 El argumento data es un data frame con los delitos ocurridos, referencia es otro data frame con las esquinas, y metros es la distancia en metros del radio de cada esquina sobre los cuales se contarán los delitos. A continuación se muestra la implementación completa de la función para todos los datos.

esquinas_delitos <- pertenencia_punto(delitos_def, intercepcion_calles_def, 150)

❗La salida es una lista donde cada elemento corresponde a un delito, el valor almacenado representa elID de la esquina asociada a cada crimen 🤔. Por ejemplo, el segundo elemento indica que el segundo delito almacenado en “delitos_def” , ocurrió dentro del radio de 150 metros de la esquina “418.” Si un delito tiene el valor “0,” significa que no tiene una esquina en las cercanías.

esquinas_delitos %>% head(3)

$`1`
[1] 515

$`2`
[1] 418

$`3`
[1] 0

👉🏻 Finalmente, la salida del algoritmo es agregada como columna al objeto “delitos_def,” permitiendo asociar a cada delito una esquina de cercanía 👌🏻. Se eliminan los delitos que no tengan ninguna esquina asociada, posteriormente se agrupan las esquinas con los distintos meses de estudio, con el objetivo de contabilizar cuántos delitos ocurrieron en cada esquina durante cada mes de estudio. Inicialmente, luego de limitar el periodo de estudio y filtrar los tipos de delitos, habían 231.514 delitos, de los cuales 182.346 tienen asociada una esquina.

delitos_def_1 <- delitos_def %>% 
  bind_cols(esquina = esquinas_delitos %>% unlist()) %>% 
  filter(esquina != 0) %>% 
  mutate(esquina = paste0("esquina_",esquina)) %>% 
  group_by(esquina,fecha_split) %>% 
  tally(name = "delitos")

knitr::kable(delitos_def_1 %>% head(),
             caption = "Conteo de delitos ocurridos en cada combinación de esquinas y fechas")

Table 3: Conteo de delitos ocurridos en cada combinación de esquinas y fechas
esquina	fecha_split	delitos
esquina_1	Dec 2017	1
esquina_1	Jan 2018	8
esquina_1	Feb 2018	3
esquina_1	Mar 2018	10
esquina_1	Apr 2018	4
esquina_1	May 2018	6

Indexación de variables meteorológicas 🌤️≫🚔

☔Para cada fecha se agregan los 4 factores meteorológicos , los valores de estos atributos serán iguales para todas las esquinas, el único elemento que genera cambio en los valores es la fecha. A continuación, se muestra el histórico mensual de delitos ocurridos y de registro meteorológicos para cada esquina, pero almacenando toda la información en una fila, es decir, cada observación representa una esquina uncia.

Show code

delitos_def_2 <- delitos_def_1 %>% 
  pivot_wider(names_from = fecha_split,
              values_from = delitos, 
              values_fill = 0) %>% 
  pivot_longer(!esquina, names_to = "fecha_split", values_to = "delitos") %>% 
  mutate(fecha_split = as.yearmon(fecha_split, "%b %Y")) %>% 
  left_join(
    temperatura_mes_histoico %>% 
      select(fecha, media), by = c("fecha_split"="fecha")
    ) %>% 
  left_join(
    lluvia_mes_histoico, by = c("fecha_split"="fecha")
    ) %>% 
  left_join(
    viento_mes_histoico, by = c("fecha_split"="fecha")
  ) %>% 
  ungroup() %>% 
  dplyr::rename(temperatura = media, mm_agua = mm, dias_lluvia = dias)

delitos_def_3 <- delitos_def_2 %>% 
  pivot_wider(names_from = fecha_split,
              values_from = c(delitos, temperatura, mm_agua, dias_lluvia, veloc_viento), 
              values_fill = 0) %>% 
  clean_names()

rmarkdown::paged_table(delitos_def_3 %>% head())

Depuración de esquinas con baja incidencia delictiva 👁️

Las esquinas presentan 24 meses de observación, pero no en todos los meses registran delitos ocurridos, se debe prestar atención a aquellas esquinas que registran una baja ocurrencia de delitos en los meses observados.

Show code

delitos_def_4 <- delitos_def_3 %>% 
 dplyr::mutate(meses_sin_delitos = rowSums(dplyr::select(., delitos_dec_2017:delitos_dec_2019) == 0)) %>% 
  relocate(meses_sin_delitos)

delitos_def_4$meses_sin_delitos %>% 
  table() %>% 
  as.data.frame() %>% 
  ggplot(aes(x=., y=Freq)) +
  geom_segment( aes(x=., xend=., y=0, yend=Freq)) +
  geom_point(size=4, color="red", alpha=0.7, shape=21, stroke=2) +
  theme_minimal() +
  labs(x = "Meses sin delitos", y = "Número de esquinas",
       title = "Distribución de esquinas con meses sin delitos registrados")

Figure 10: Esquinas con meses sin delitos registardos

❗Mediante la distribución mostrada en la Figura 10, se observa que 571 esquinas registraron delitos en todos los meses. El 88% de las esquinas no mostraron delitos como máximo en 15 de los 24 meses analizados, aquellas esquinas que no presentan delitos en más de 15 meses se excluyen del estudio. Esto debido a que pueden acarrear un desbalance excesivo sobre la variable objetivo ⚖️, realizando la depuración quedan 2023 esquinas definitivas.

delitos_def_5 <- delitos_def_4 %>% 
  filter(meses_sin_delitos < 16) %>% 
  select(-meses_sin_delitos)

Indexación de variables de entorno 🌏≫🚔

👉🏻 Para indexar todas las variables de entorno, se recupera la longuitud y latitud de las esquinas, ya que es el insumo primcipal para vincular los datos. El factor de zonas de inclusion urbana, corresponde a poligonos espaciales que delimitan zonas con determiandas carencias dentro de la ciudad. A continuación, se procede a determinar para cada esquina, si esta geolocalizada sobre alguno de estos poligonos.

delitos_def_6 <- delitos_def_5 %>% 
  left_join(
    intercepcion_calles_def %>% 
      mutate(esquina = paste0("esquina_",id)) %>% 
      select(-id)
    ) %>% 
  relocate(esquina,long,lat)

#Trnasformacion de coordenadas de esquinas en objetos espaciales
puntos_esquinas <- SpatialPointsDataFrame(coords = delitos_def_6[c("long", "lat")],
                                          data = delitos_def_6[c("long", "lat")],
                                          proj4string = CRS(proyeccion)) 

# Puntos dentro de los poligonos
punto_poligono <- point.in.poly(puntos_esquinas, inclusion_urbana)

punto_poligono <- punto_poligono@data

delitos_def_7 <- delitos_def_6 %>% 
  bind_cols(
    punto_poligono %>% 
      select(Id)
  ) %>% 
  mutate(inclu_urbana = ifelse(is.na(Id),0,1)) %>% 
  select(-Id)

Distancias entre esquinas y factores de entorno físico 📐

👉🏻 El siguiente grupo de variables a indexar corresponde al resto de los elementos de entorno, para ello se debe calcular el número de observaciones de cada variable de entorno, que están dentro de un radio inferior a 250 metros de cada esquina . Esto requiere la construcción de un algoritmo ⚙️, el mismo debe automatizar el cálculo de distancias y el conteo de eventos. A continuación se muestra la función a implementar.

F.distancia.punto <- function(referencia, data,metros){
  
  conteo <- lapply(data, function(x){
    
    vector_localizacion <- distm(x, referencia, fun = distHaversine)
    n_total <- sum(vector_localizacion <= metros)
    n_total
  })
  conteo  %>% unlist()
}

❗La función recibe 3 argumentos, las esquinas en formato de lista con la geolocalización, los conjuntos de datos de entorno almacenados en una lista y la distancia en metros.

split_esquina <- delitos_def_7 %>% 
  select(long,lat) %>% 
  split(1:nrow(delitos_def_7))

datas <- erer::listN(
  universidades,
  local_bailables,
  wifi,
  hotel_baja,
  hotel_alta,
  cine,
  teatros,
  hogares_paradores,
  gimnasios,
  hospitales,
  metrobus,
  centros_salud_comunitarios,
  colectivo,
  barrios_populares_centroides,
  centros_de_salud_privados,
  comisarias,
  embajadas,
  bocas_de_subte,
  atm,
  bancos,
  gasolina,
  gastronomica,
  taxi,
  carteles,
  farmacias,
  centros_medicos_barriales,
  estadios,
  skate,
  garajes,
  educativos,
  bibliotecas,
  basilicas,
  organizaciones_sociales,
  estaciones_ferrocarril,
  bus_turistico,
  consulados,
  bomberos
  )

datas <- lapply(datas, function(x){
  
  x %>% select(c("long","lat"))
  
  })

variables_entorno <- lapply(datas, F.distancia.punto, metros = 250, data = split_esquina) 

variables_entorno2 <- bind_cols(variables_entorno %>% as.data.frame())

Data_set_final_1 <- bind_cols(delitos_def_7,variables_entorno2)

rmarkdown::paged_table(Data_set_final_1 %>% head())

🔎 Por cada esquina se puede observar el conteo del número de elementos de entorno.

👉🏻🔝🔝🔝En el caso de querer reproducir el estudio y saltarse toda la etapa de preprocesamiento mostrada, se puede acceder al dataset preprocesado guardado en el paquete 📦 sknifedatar, mediante la siguiente sintaxis.

# install.packages("devtools")
devtools::install_github("rafzamb/sknifedatar")
library(sknifedatar)

data(data_crime_clime)

Aplicación del “método” de ventanas deslizantes 🤔🧐

💡 Inspiración

🔹El proceso de transformación de información para llegar al conjunto de datos definitivo, fue inspirado en gran medida sobre el estudio que presentaron Lin, Yen, and Yu (2018) . Estos investigadores incorporan el factor temporal a los elementos demográficas y de entorno, pero evaluando los tres factores de manera simultánea para construir un modelo de aprendizaje automático, con el objetivo de capturar desde una perspectiva más completa el entorno criminal.

🔹El análisis fue realizado en una ciudad de Taiwán, dividiendo el mapa de la ciudad en cuadrículas de 200 x 200 metros, posteriormente calculaban el número de delitos ocurridos en el polígono. Sin embargo, el presente estudio toma este principio, pero en lugar de cuadrículas son seleccionadas esquinas, para luego calcular los delitos que ocurrieron en un radio de metros de cercanía. El enfoque que se utiliza en este proyecto se considera mucho más realista 👀, ya que a priori las esquinas están localizadas en zonas habitadas o mínimamente transitadas, el enfoque de cuadrículas contabiliza polígonos únicamente por estar dentro del espacio geográfico de la ciudad.

📖 Funcionamiento

🔹Aunque en este estudio el método fue denominado “ventana deslizante” , la idea fue tomada de un planteamiento similar mostardo en Lin, Yen, and Yu (2018). Este consiste en seleccionar el mes intermedio “t” de todo el periodo de estudio, posteriormente se calcula el número de delitos ocurridos para cada esquina en el mes anterior, en los últimos 3 meses, 6 meses, 12 meses y el mismo mes del año anterior. La variable objetivo sería el número de delitos ocurridos en el mes de estudio, además este mismo procedimiento es aplicado con las variables climáticas.

Figure 11: Funcionamiento del enfoque de ventana deslizante. Fuente: Elaboración propia

🔹El procedimiento descrito anteriormente es replicado de manera móvil, es decir, rodando la ventana temporal desde \(t+1\) hasta \(n\), siendo \(n\) el último mes de estudio. Estos \(n-t\) conjuntos de datos se unifican en un dataset único, representando el conjunto de entrenamiento, el conjunto número \(n\) hace referencia a los datos sobre los cuales se hartan las predicciones finales. Para este estudio, el último dataset está representado por el mes de diciembre del año 2019.

⏯ Implementación

🔎 Este enfoque está automatizado mediante la creación de la función sliding_window almacenada dentro del paquete 📦 sknifedatar. La función toma 4 argumentos, que se citan a continuación:

Data: data frame con el conteo historico por variables (Estructura como la mostrada en la Figura 11 .
inicio: representa el mes inicial donde \(t = 1\).
pliegues: vector que inicia en 1 y termina en el número de periodos que se desee recorrer.
variables: una palabra o vector que permita seleccionar las variables en conjunto e implementar la función para cada grupo.

A continuación, se aplica la función sliding_window sobre los delitos y los 4 factores meteorológicos.

library(sknifedatar)
pliegues <- 1:13
names(pliegues) <- pliegues

variables <- c("delitos", "temperatura", "mm_agua", "lluvia", "viento")
names(variables) <- variables

Data_set_final_2 <- sliding_window(data = Data_set_final_1 %>% select(-c(long,lat)),
                                  inicio = 13,
                                  pliegues = pliegues,  
                                  variables = variables)

rmarkdown::paged_table(Data_set_final_2 %>% head())

Unificación y seleccion final de variables

👉🏻 Al conjunto de esquinas con los recuentos originados por la función sliding_window, se agregan las variables de entorno previamente calculadas. Adicionalmente, para los factores meteorológicos fueron seleccionados únicamente los conteos del mes anterior y de los últimos 3 meses, esto debido a que la inclusión del resto de los meses resultó en una disminución en las métricas de evaluación de los modelos.

Data_set_modelos <- Data_set_final_2 %>% 
  left_join(
    Data_set_final_1 %>% 
      select(1,names(datas),c(long,lat)), by =c("id" = "esquina")
    ) %>% 
  select(
    c(id,pliegue,contains("delitos"), mm_last_3, mm_last_1, temperatura_last_3, temperatura_last_1, dias_last_3,
      dias_last_1, veloc_last_3, veloc_last_1, 33:(ncol(.)))
  )

☑️ Finalmente está terminada la etapa de manipulación y transformación de datos 👌🏻👌🏻👌🏻. Aunque en el último fragmento de código se realiza una selección de variables, dicha selección surge de una gran cantidad de iteraciones, evaluando el rendimiento y “sentido de negocio” de las predicciones de los modelos.

Modelado 📚💻🎯

🔸El estudio será abordado desde tres enfoques en el siguiente orden: series de tiempo, regresión y clasificación. El conjunto de validación es el mes de diciembre del año 2019, en la práctica este mes sería la observación con \(t = 13\) para cada una de las esquinas, y el conjunto de entrenamiento viene dado por las observaciones de \(t = 1\) hasta \(t= 12\) , siguiendo el ejemplo de la Figura 11 .En el siguiente diagrama se muestra el flujo de trabajo para modelar el problema.

Figure 12: Diagrama refrente al flujo de trabajo. Fuente: Elaboración propia

Baseline: Series de tiempo 📉

🔶Este enfoque consiste en el establecimiento de un modelo base, con esto se implementa un algoritmo poco complejo que establezca el límite inferior de rendimiento, se aplicará el modelo de series de tiempo Croston. Este es un modelo especial para regresiones de conteo, está diseñado para abordar problemas de predicciones de series intermitentes, tomando en cuenta únicamente el histórico de la ocurrencia de delitos, Hyndman (2018) (sección 12.2 Timeseries of counts, párr.4) expone el funcionamiento del modelo. A continuación, se muestra la preparación de los datos para introducirlos en el modelo.

base_line_ts <-  Data_set_final_1 %>% 
  select(esquina, contains("delitos")) %>% 
  pivot_longer(!esquina, names_to = "Tiempo", values_to = "Delitos") %>% 
  mutate(Tiempo =  gsub("delitos_","",Tiempo)) %>% 
  mutate(Tiempo = as.Date(parse_date_time2(Tiempo, orders = "mY"))) %>% 
  filter(Tiempo < "2019-12-01")

Table 4: Preparación de datos para el modelado de series de tiempo
esquina	Tiempo	Delitos
esquina_1	2017-12-01	1
esquina_1	2018-01-01	8
esquina_1	2018-02-01	3
esquina_1	2018-03-01	10
esquina_1	2018-04-01	4
esquina_1	2018-05-01	6

Transformación de datos en objetos de series de tiempo

list_ts <- base_line_ts %>% 
  split(base_line_ts$esquina)

ts_object <- lapply(list_ts, function(x) ts(x$Delitos, start=c(2017,12),frequency=12))

tail(ts_object,2)

$esquina_998
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2017                                               5
2018   4   4   5   3   0   2   3   2   4   3   6   3
2019   3   3   1   1   2   2   4   2   4   3   1    

$esquina_999
     Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
2017                                               0
2018   2   2   2   0   0   0   0   0   0   1   0   1
2019   1   1   0   0   1   5   0   0   1   1   2

❗El objeto “ts_object”es una lista que almacena en cada elemento el registro histórico de los delitos ocurridos por esquina.

Figure 13: Registro histórico mensual de delitos ocurridos en la esquina -999-

📎 En la Figura 13 ,se observan muchos ceros presentes en los datos, es decir, las esquinas presentan meses sin registrar delitos. Este tipo de series no suelen presentar una tendencia definida, además de ser bastante volátiles, lo que dificulta realizar pronósticos con métodos convencionales. Según Hyndman (2018), aunque el modelo Croston fue diseñado originalmente para el pronóstico de demanda en un período determinado, se usa habitualmente para tratar este tipo de series. La implementación se hará mediante el paquete 📦 tsintermittent⁸ .

🔹 Modelos de series de tiempo

Para todas las esquinas se utilizará la media de las series y el modelo Croston para predecir el mes de diciembre del 2019, el modelo se aplicará con los parámetros que tiene por defecto.

Media de las series

ts_mean <- lapply(ts_object , mean)

Modelo Croston

ts_crost <- lapply(ts_object , function(x) crost(x,type='croston',h=1)$components$c.out[,1])

Unificación de predicciones

Baseline_TS <- bind_cols(
  delitos_dec_2019 = Data_set_final_1$delitos_dec_2019,
  data.frame(mean_predict = unlist(ts_mean)),
  bind_rows(ts_crost)
  ) %>% 
  rename(crost_predict = Demand)

rmarkdown::paged_table(Baseline_TS)

Evaluación de Errores 🎯

A través del paquete 📦 yardstick⁹ perteneciente a tidymodels, se evaluaran las predicciones mostradas, sin embargo cuando se comparan predicciones generadas fuera del flujo de tidymodels, no es posible aplicar de manera automatizada las funciones de yardstick (al menos no encontre alguna solución en la documentación🤔).

❗Para solventar este problema fue creada la función multieval está permite utilizar las métricas de yardstick sobre las predicciones de los modelos. Una parte importante de la codificación de multieval , esta almacenada en el paquete 📦 sknifedatar.Para consultar en detalle el proceso de construcción de la función y ver algunos tópicos de programación funcional en R, consultar la publicación previa “Evaluación de predicciones mediante programación funcional en R” Zambrano (2021) .

library(sknifedatar)
Metricas_ts <- multieval(.dataset = Baseline_TS,
                         .observed = "delitos_dec_2019",
                         .predictions = c("mean_predict","crost_predict"),
                         .metrics = listn(rmse, rsq, mae))

Metricas_ts$summary_table %>% 
  knitr::kable(caption = "Métricas modelo Croston y medias de las series") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 5: Métricas modelo Croston y medias de las series
modelo	rmse	rsq	mae
mean_predict	2.458516	0.7520300	1.574353
crost_predict	2.429252	0.7668506	1.697229

👉🏻👉🏻👉🏻 Estos resultados se contrastan con las predicciones del modelo de regresión de la siguiente sección, evaluando cual de las predicciones presenta un mejor resultado. Existen métodos de series de tiempo más avanzados para estudiar este tipo de fenómenos, pero están fuera del alcance de este proyecto, una introducción acerca del tema puede encontrarse en Christou and Fokianos (2013) .

Modelos de aprendizaje automático 🚀

🔹Para los modelos de aprendizaje automático, se implementan dos modelos basados en árboles, siendo el Extreme Gradient Boosting (XGB) (Chen and Guestrin 2016) y el Light Gradient Boosting Machine (LightGBM) (Ke et al. 2017) .

Partición de datos ➗

📌 En este estudio no hace falta realizar una partición aleatoria en datos de entrenamiento y validación, debido a que ambos conjuntos de datos se obtienen a través de las fechas.

train <- Data_set_modelos %>% filter(pliegue != 13)

test <- Data_set_modelos %>% filter(pliegue == 13)

ind <- list(analysis = seq(nrow(train)),assessment = nrow(train) + seq(nrow(test)))

splits <- make_splits(ind, Data_set_modelos)
splits$id <- tibble(id = "Resample1")

train <- training(splits)
test <- testing(splits)

splits

<Analysis/Assess/Total>
<24276/2023/26299>

Validación cruzada 🧮

🔎Se aplicará una validación cruzada de diez particiones utilizando la técnica “k-Fold-Cross-Validation (CV),” adicionalmente los hiperparámetros de los modelos se ajustarán a través de una optimización bayesiana, con un máximo de 15 iteraciones y una restricción de finalización de 10 iteraciones sin mejora.

set.seed(123)
cv_folds <- vfold_cv(train, v = 10, strata = delitos)
ctrl <- control_bayes(no_improve = 10, verbose = T, save_pred = T)

Enfoque de regresión 📈

🔎 Para cada esquina se busca pronosticar el número de delitos que ocurrirán, este enfoque se puede contrastar directamente con el modelo base anterior. Se requiere previamente analizar la distribución de la variable objetivo, con la finalidad de evaluar la aplicación alguna transformación en la misma, y verificar si dicha transformación mejora el rendimiento del modelo.

Show code

g1 <- ggplot(Data_set_modelos, aes(x=delitos)) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  geom_density(alpha=.2, fill="#FF6666")+ 
  theme_minimal() +
  labs(x = "", y = "")

# qqplot
g3 <- ggplot(Data_set_modelos, aes(sample = delitos))+ 
  stat_qq()+
  stat_qq_line()+  
  theme_minimal() +
  labs(x = "", y = "")

plot_original <- g3 | g1 

g1 <- ggplot(Data_set_modelos, aes(x= log1p(delitos))) + 
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  geom_density(alpha=.2, fill="#FF6666")+ 
  theme_minimal() +
  labs(x = "", y = "")

# qqplot
g3 <- ggplot(Data_set_modelos, aes(sample = log1p(delitos)))+ 
  stat_qq()+
  stat_qq_line()+  
  theme_minimal() +
  labs(x = "", y = "")

plot_trnasform <- g3 | g1 

ggpubr::ggarrange(plot_original + plot_annotation(title = 'Deltios originales'),
                  plot_trnasform + plot_annotation(title = 'Deltios con la transformación log1p'),
                  ncol = 1)

Figure 14: Trandsformación de la variable de ocurrencia de delitos

⚠️ Analizando la Figura 14 ,se observa una fuerte asimetría en la distribución de la variable, en la aplicación de modelos de aprendizaje automático suele ser una práctica común utilizar transformaciones sobre la variable respuesta, con el objetivo de intentar mejorar el poder predictivo de los modelos minimizando la asimetría de la variable.

❗Se aplica la transformación logarítmica con la adición de 1. Aunque la transformación minimizó la asimetría de la variable, aún se observa un importante sesgo en la misma. Además dicha transformación no reflejó una mejora en las métricas de evaluación de los modelos, por lo tanto, esta transformación no será aplicada sobre los datos.

Receta de preprocesamiento 🥘

Creación de las variables bancos_atm y teatros_cine unificando bacnos con atm y teatros con cine , esto debido a que dichas variables por sí solas presentan una baja varianza.
Se remueven temporalmente las 4 variables utilizadas en el paso anterior, pero también se elimina el id, pliegue, lat y long, la exclusión de la geolocalización resultó en una leve mejora en el ajuste de los modelos de “regresión.”
Finalmente se eliminan los predictores con varianza cercana a 0, es decir variables muy escasas y desequilibradas. A excepción de los atributos estaciones_ferrocarril, bocas_de_subte, hotel_baja, bancos_atm, teatros_cine, ya que estaban muy cerca de los umbrales de tolerancia para la eliminación y resultaron contribuir positivamente a mejorar las predicciones. Para ver en detalle el funcionamiento de este paso de preprocesamiento, puede consultar la función step_nzv en el paquete 📦 recipes¹⁰ de tidymodels.

Receta_regresion <- 
  
  recipe(formula = delitos ~ ., data =  train
         ) %>%
  
  step_mutate(bancos_atm = bancos + atm, teatros_cine = teatros + cine
              ) %>% 
  
  step_rm(bancos, atm, teatros, cine, id, pliegue, lat, long
          ) %>% 
  
  step_nzv(all_predictors(), -c(estaciones_ferrocarril,bocas_de_subte,hotel_baja,bancos_atm,teatros_cine))

❗Los 3 pasos de preprocesamiento descritos anteriormente, se reflejan en las 3 step_* functions del fragmento de código.

Definición del modelo ✍️

xgb_spec_reg <- boost_tree(
  trees = tune(), 
  tree_depth = tune(), 
  min_n = tune(), 
  loss_reduction = tune(),                     
  sample_size = tune(),
  mtry = tune(),         
  learn_rate = tune()                         
) %>% 
  set_engine("xgboost") %>% 
  set_mode("regression")

xgb_param_reg <- 
  xgb_spec_reg %>%
  parameters() %>% 
  update(mtry = mtry(c(28L, 35L)),
         trees = trees(c(900L, 1500L)),
         sample_size = sample_prop(c(.6, .99))
  ) 

xgb_wf_reg <- workflow() %>%
  add_recipe(Receta_regresion) %>%
  add_model(xgb_spec_reg)

Optimización de hiperparámetros 🕹️

🔹El mae es la métrica que se optimizará en la búsqueda de los mejores hiperparámetros, esto debido a que se sabe de la existencia de esquinas que presentan una alta volatilidad en el registro de delitos. Aunque son pocas las observaciones que presentan esta característica, afectan de manera severa la evaluación global de los modelos, por lo tanto, no interesa penalizar excesivamente los errores en estos datos atípicos.

set.seed(1234)
doParallel::registerDoParallel(cores = parallel::detectCores() -1)

xgb_reg <- tune_bayes(xgb_wf_reg, 
                      resamples = cv_folds,
                      iter = 15,
                      param_info = xgb_param_reg,
                      metrics = metric_set(mae),
                      control = ctrl,
                      initial = 8)

doParallel::stopImplicitCluster()

Show code

autoplot(xgb_reg, type = "performance") +
  labs(title = "Evolución del error durante el ajuste de hiperpárametros") +
  theme_minimal()

Figure 15: XGB / Evolución del error durante el ajuste de hiperpárametros

Ajuste sobre los mejores hiperparámetros

Show code

knitr::kable(select_best(xgb_reg, "mae"),
             caption = "XGB / Mejores hiperpárametros") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 6: XGB / Mejores hiperpárametros
mtry	trees	min_n	tree_depth	learn_rate	loss_reduction	sample_size	.config
35	1017	6	4	0.0028825	14.34346	0.6036219	Iter8

modelo_xgb_reg <- finalize_workflow(xgb_wf_reg,
                                    parameters = select_best(xgb_reg, "mae"))

Errores de valdiacion cruzada

validacion_xgb_reg  <- fit_resamples(
  object       = modelo_xgb_reg,
  resamples    = cv_folds,
  metrics      = metric_set(rmse, mae,rsq),
  control      = control_resamples(save_pred = FALSE)
)

Errores sobre los datos de prueba

final_xgb_reg <- last_fit(modelo_xgb_reg, splits)

Comparación de errores 🔎

Show code

bind_rows(
  collect_metrics(validacion_xgb_reg) %>% 
    dplyr::rename(valor = mean) %>% 
    mutate(Error = "validación cruzada"),
  collect_predictions(final_xgb_reg) %>%
    metrics(truth = delitos, estimate = .pred) %>% 
    dplyr::rename(valor = .estimate) %>% 
    mutate(Error = "Diciembre 2019")) %>% 
  mutate(modelo = "XGB regresión") %>% 
  dplyr::select(modelo,.metric,valor,std_err,Error) %>% 
  arrange(.metric) %>% 
  knitr::kable(caption = "XGB regresión / Comparación de errores")

Table 7: XGB regresión / Comparación de errores
modelo	.metric	valor	std_err	Error
XGB regresión	mae	1.5624788	0.0117539	validación cruzada
XGB regresión	mae	1.5034106	NA	Diciembre 2019
XGB regresión	rmse	2.3306645	0.0352261	validación cruzada
XGB regresión	rmse	2.1774435	NA	Diciembre 2019
XGB regresión	rsq	0.7796652	0.0069891	validación cruzada
XGB regresión	rsq	0.8026857	NA	Diciembre 2019

Show code

predicciones_reg <- final_xgb_reg %>%
  collect_predictions() %>% 
  dplyr::select(-.row)

plot_reg <- predicciones_reg %>% 
  ggplot(aes(x = .pred, y = delitos))+
  geom_point()+ 
  geom_abline(intercept = 0, col = "red") +
  theme_minimal()

densidad <- predicciones_reg %>% 
  pivot_longer(!c(id,.config),names_to = "Variable", values_to = "n") %>% 
  ggplot(aes(x=n, volor = Variable, fill = Variable)) + 
  geom_density(alpha=.2)+ 
  theme_minimal() +
  labs(x = "", y = "")

plot_reg / densidad

Figure 16: XGB / Visualización de predicciones

💡En la Tabla 7 , se observa una mejoría de las tres métricas al comparar los errores de validación cruzada con respecto a los del mes de diciembre del 2019, los cual es un signo de que el modelo a priori no está sobreajustado ✅. La diferencia existente en las métricas rmse y mae, se debe principalmente a 👉🏻👉🏻esquinas que presentaron saltos abruptos en el registro de delitos para el mes de diciembre.

Comparación con el Baseline 🏁

👉🏻A continuación, se muestran los gráficos comparando los registros observados con los predichos a través de los modelos Croston , XGB y las medias de las series, además, una tabla resumen con las principales métricas.

Show code

p1 <- final_xgb_reg %>%
  collect_predictions() %>%
  ggplot(aes(delitos, .pred)) +
  geom_abline(intercept = 0, col = "red") +
  geom_point(alpha = 0.6, color = "midnightblue") +
  theme_minimal() +
  ggtitle("XGB") +
  scale_y_comma() +
  scale_x_comma()

p2 <- Baseline_TS %>%
  ggplot(aes(delitos_dec_2019, mean_predict)) +
  geom_abline(intercept = 0, col = "red") +
  geom_point(alpha = 0.6, color = "midnightblue") +
  theme_minimal() +
  ggtitle("Mean") +
  scale_y_comma() +
  scale_x_comma() +
  xlab("delitos")

p3 <- Baseline_TS %>%
  ggplot(aes(delitos_dec_2019, crost_predict)) +
  geom_abline(intercept = 0, col = "red") +
  geom_point(alpha = 0.6, color = "midnightblue") +
  theme_minimal() +
  ggtitle("Croston") +
  scale_y_comma() +
  scale_x_comma() +
  xlab("delitos")

ggdraw() +
  draw_plot(p2, x = 0, y = .5, width = .5, height = .5) +
  draw_plot(p3, x = .5, y = .5, width = .5, height = .5) +
  draw_plot(p1, x = .25, y = 0, width = .5, height = .5)

Figure 17: Evaluacion de las predicciones del enfoque de regresión

Show code

bind_rows(
  final_xgb_reg %>%
    collect_predictions() %>%
    metrics(truth = delitos, estimate = .pred) %>% 
    mutate(model = "XGB") %>% 
    dplyr::select(-2) %>% 
    pivot_wider(names_from = .metric, values_from = .estimate),
  Metricas_ts$summary_table) %>% 
  knitr::kable(caption = "Comparación de métricas del enfoque de regresión") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 8: Comparación de métricas del enfoque de regresión
model	rmse	rsq	mae	modelo
XGB	2.177443	0.8026857	1.503411	NA
NA	2.458516	0.7520300	1.574353	mean_predict
NA	2.429252	0.7668506	1.697229	crost_predict

📌 En la Figura 17 se aprecia que el modelo XGB fue superior en las tres métricas de evaluación, visualmente también se observa un mejor ajuste. Existen esquinas que presentaron cambios abruptos en el registro de delitos para el mes de evaluación, afectando el rendimiento de los modelos. Las regresiones de conteo son una alternativa para modelar este tipo de problemas 💡, se podrían contrastar los resultados con este tipo de enfoque alternativo, Hilbe (2017) presenta una explicación del tema.

Ajuste final del modelo ☑️

👉🏻Para fines didácticos, se ajusta el modelo sobre todo los datos, simulando el ejercicio de tener un modelo listo para predecir sobre nuevos datos, por ejemplo, para predecir la ocurrencia de delitos en el mes de enero del 2020. En este caso será simplemente sobre una muestra de 6 esquinas, elimnando la columna de delitos (no es necesario), se recuerda que el modelo es en realidad un objeto de tipo “workflow,” por lo tanto los datos a predecir no necesitan ser pre procesados, ya que el “workflow” tiene entrenado el preprocesamiento realizado.

modelo_fit_reg <- modelo_xgb_reg %>% fit(Data_set_modelos)

predict(modelo_fit_reg, Data_set_modelos %>% head() %>% dplyr::select(-delitos))

# A tibble: 6 x 1
  .pred
  <dbl>
1 5.85 
2 7.35 
3 1.85 
4 0.988
5 2.23 
6 1.06

Enfoque de clasificación 📊

❗El enfoque de clasificación será abordado y contrastado desde dos perspectivas, contraponiendo las mismas para evaluar cual resulta una mejor aproximación a la problemática que se estudia. Se implementará el modelo LightGBM para ambas pespectivas, mediante el paquete 📦 treesnip¹¹ que permite utilizar el paquete 📦 lightgbm¹² dentro del marco de tidymodels .

Criterios de discretización

Table 9: Criterios de discretización de la ocurrencia de delitos
Condición	Categoría	Frequencia
Si los delitos son <= 1	Bajo	9190
Si los delitos son <= 4 y >= 2	Medio	10511
Si los delitos son >= 5	Alto	6598

📌En la Tabla 9 se agrupan las esquinas según el nivel de riesgo, distribuido en tres niveles, el objetivo reside en implementar modelos que puedan discriminar las tres categorías.

Evolución del nivel de riesgo de las esquinas de CABA. Elaborado con el paquete gganimate [@gganimate]

Figure 18: Evolución del nivel de riesgo de las esquinas de CABA. Elaborado con el paquete gganimate (Pedersen and Robinson 2020)

👉🏻 📅🚨En la Figura 18 se observa la evolución de la variación del nivel de riesgo de las esquinas durante el periodo de estudio .

❗Es importante destacar 2 limitaciones planteadas sobre el estudio, primero, las zonas periféricas de la ciudad probablemente registren tasas similares o superiores a la media de la ciudad, sin embargo, solo se cuenta con los delitos denunciados y registrados 🧐, pudiendo ser una aproximación alejada de la realidad para estas zonas.

👉🏻Segundo, cuando se capturan los elementos de entorno (estaciones de bombero, restaurantes, farmacias cercanas, …), es posible (muy raro) que algunos cambien o dejen de estar presentes físicamente en algún intervalo del periodo de estudio, se asume que es invariante la ubicación física de los elementos.

Clasificación multiclase 🔴⚫🔵

🔹 Utilizando el criterio de discretización presentado en la Tabla 9, se ajustará un modelo de clasificación para generar predicciones multiclase referentes al nivel de riesgo de las esquinas.

Receta 🥘

🔹 Se utilizará la misma receta de preprocesamiento aplicada en el enfoque de regresión, aunque con una ligera modificación. Se debe agregar un paso extra al final de la receta, esto con el objetivo de discretizar los delitos ocurridos. Se probó realizar un balanceo de clases a través de un método de sobremuestreo, con la finalidad de mejorar las predicciones.

👉🏻 Se implemento el algoritmo de Synthetic Minority Oversampling Technique (SMOTE), para mayores detalles consultar a Chawla et al. (2002) . Sin embargo, las observaciones sintéticas generadas no mejoraron el rendimiento del modelo, por lo tanto, se ajustó el modelo definitivo sin balancear las clases.

Receta_multiclass <- Receta_regresion %>% 
  step_mutate(delitos = case_when(delitos <= 1 ~ "Bajo",
                                 delitos <= 4  ~ "Medio",
                                 delitos >= 5 ~ "Alto"))

📌PD: Observadno el fragmento de código anterior ,la creación de la receta para la clasificación multiclase se genera reutilizando la receta de regresión, simplemente agregando un "%>%" seguido de la discretización de la variable dependiente

Especificación del modelo ✍️

lgb_spec_mclass <- boost_tree(
  trees = tune(), 
  tree_depth = tune(),
  min_n = tune(), 
  loss_reduction = tune(),                     
  sample_size = tune(),
  mtry = tune(),       
  learn_rate = tune() 
) %>% 
  set_engine("lightgbm", verbose=-1) %>% 
  set_mode("classification")

lgb_param_mclass <- 
  lgb_spec_mclass %>%
  parameters() %>% 
  update(mtry = mtry(c(20L, 30L)),
         trees = trees(c(900L, 1500L)),
         sample_size = sample_prop(c(.6, .99))) 

lgb_wf_mclass <- workflow() %>%
  add_recipe(Receta_multiclass) %>%
  add_model(lgb_spec_mclass)

Optimización de hiperparámetros 🕹️

🔹 La métrica optimizada para la búsqueda de hiperparámetros corresponde a la roc_auc en la implementación Hand-Till, debido a que se cuenta con una distribución de clases desequilibrada (pero no excesiva). Dentro del marco tidymodels se expone que 👉🏻👉🏻este método no hace suposiciones específicas sobre la distribución de clases o los costos de clasificación erróneas, es decir, es insensible a distribuciones de clases. El objetivo es realizar una evaluación posterior de los modelos a través de matrices de confusión, agregando costos a los errores. Para más detalles, Hand and Till (2001) exponen la implementación del método.

set.seed(1234)
doParallel::registerDoParallel(cores = parallel::detectCores() -1)

lgb_mclass <- tune_bayes(lgb_wf_mclass, 
                         resamples = cv_folds,
                         iter = 15,
                         param_info = lgb_param_mclass,
                         metrics = metric_set(roc_auc),
                         control = ctrl,
                         initial = 8)

doParallel::stopImplicitCluster()

Show code

autoplot(lgb_mclass, type = "performance") +
  labs(title = "Evolución del error durante el ajuste de hiperpárametros") +
  theme_minimal()

Figure 19: Lgbm multiclase / Evolución del error durante el ajuste de hiperpárametros

Ajuste sobre los mejores hiperparámetros

Show code

knitr::kable(select_best(lgb_mclass, "roc_auc"),
             caption = "Lgbm multiclase / Mejores hiperpárametros")

Table 10: Lgbm multiclase / Mejores hiperpárametros
mtry	trees	min_n	tree_depth	learn_rate	loss_reduction	sample_size	.config
20	906	22	2	0.0521606	1.129974	0.8498181	Iter5

modelo_lgb_mclass <- finalize_workflow(lgb_wf_mclass,
                                       parameters = select_best(lgb_mclass, "roc_auc"))

Errores de valdiación cruzada

validacion_lgb_mclass  <- fit_resamples(
  object       = modelo_lgb_mclass,
  resamples    = cv_folds,
  metrics      = metric_set(roc_auc),
  control      = control_resamples(save_pred = FALSE)
)

Errores sobre los datos de prueba

final_lgb_mclass <- last_fit(modelo_lgb_mclass, splits)

Comparación de errores 🔎

Show code

bind_rows(
  collect_metrics(validacion_lgb_mclass) %>%
    dplyr::rename(valor = mean) %>% 
    mutate(Error = "validación cruzada"),
  collect_metrics(final_lgb_mclass) %>%
    dplyr::rename(valor = .estimate) %>% 
    mutate(Error = "Diciembre 2019")) %>% 
    filter(.metric == "roc_auc") %>% 
    mutate(modelo = "Lgbm multiclase") %>% 
    dplyr::select(-.config) %>% 
    arrange(.metric) %>% 
  knitr::kable(aption = "Lgbm multiclase / Comparación de errores")

.metric	.estimator	valor	n	std_err	Error	modelo
roc_auc	hand_till	0.8113618	10	0.0024115	validación cruzada	Lgbm multiclase
roc_auc	hand_till	0.8147452	NA	NA	Diciembre 2019	Lgbm multiclase

Show code

Matriz_multiclase <- final_lgb_mclass %>%
  collect_predictions() %>% 
  mutate(across(c(.pred_class, delitos),
                factor, levels = c("Alto","Medio","Bajo"))) %>% 
  conf_mat(delitos, .pred_class)

autoplot(Matriz_multiclase, type = "heatmap")

Figure 20: Matriz de confusión multiclase

Show code

  collect_predictions(final_lgb_mclass) %>% 
  roc_curve(delitos, .pred_Alto:.pred_Medio) %>%
  autoplot()

Figure 21: Curvas Roc para los tres niveels de riesgo

❗Analizando los resultados , se observa un ligero aumento de la métrica de evaluación roc_auc entre la validación cruzada y los datos de prueba 🔔, significando que el modelo tuvo una buena generalización de los datos. A través de la Figura 20 , se aprecia en la matriz de confusión que los errores entre las categorías extremas son bastante bajos, siendo de 14 y 9 errores en ambos pares 👌🏻.

🔎 En la última sección del enfoque de clasificación, se agregan costos a las matrices resultantes, sin embargo, se percibe que el modelo logra discriminar satisfactoriamente una esquina peligrosa de una de bajo riesgo. Los errores se acumulan principalmente en la categoría intermedia, siendo más difícil discriminar en muchos casos donde el número de delitos está ubicado en los límites de las categorías 👈🏻, aunque se destaca que a pesar de lo mencionado se logró clasificar correctamente al 63% de esta categoría.

Predicción sobre nuevos datos 🔮

🔹Siguiendo la metodología del enfoque de regresión, se ajusta el modelo con todos los datos y se predice para un hipotético nuevo set de datos sin la variable respuesta (Esto se hace de manera intencional para analizar cómo funcionan en estas situaciones losworkflows). En este caso, la receta del workflow realiza la discretización de la variable dependiente, sin embargo, al eliminar dicha variable el proceso de predicción generará un error, ya que no encuentra este atributo.

tidy(Receta_multiclass) %>% 
  knitr::kable(caption = "Receta multiclase") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 11: Receta multiclase
number	operation	type	trained	skip	id
1	step	mutate	FALSE	FALSE	mutate_SMSbs
2	step	rm	FALSE	FALSE	rm_MR9Ts
3	step	nzv	FALSE	FALSE	nzv_IusHH
4	step	mutate	FALSE	FALSE	mutate_pjl57

🔎Para solventar esta situación, se actualiza el argumento skip a TRUE para el paso de preprocesamiento número 4 correspondiente a la discretización, esto permite que esta step_* function no se aplique al momento de generar predicciones.

Receta_multiclass$steps[[4]] <- update(Receta_multiclass$steps[[4]], skip = T)

tidy(Receta_multiclass) %>% 
  knitr::kable(caption = "Receta multiclase con skip = TRUE en la step 4") %>%
  kableExtra::kable_styling(full_width = FALSE)

Table 12: Receta multiclase con skip = TRUE en la step 4
number	operation	type	trained	skip	id
1	step	mutate	FALSE	FALSE	mutate_SMSbs
2	step	rm	FALSE	FALSE	rm_MR9Ts
3	step	nzv	FALSE	FALSE	nzv_IusHH
4	step	mutate	FALSE	TRUE	mutate_pjl57

❗Se observa el step 4 actualizado con el parámetro skip = TRUE.

model_fit_multiclass <- update_recipe(modelo_lgb_mclass, Receta_multiclass) %>% 
  fit(Data_set_modelos)

predict(model_fit_multiclass, Data_set_modelos %>% head() %>% dplyr::select(-delitos))

# A tibble: 6 x 1
  .pred_class
  <fct>      
1 Alto       
2 Alto       
3 Medio      
4 Bajo       
5 Medio      
6 Bajo

Clasificación binaria múltiple

🔴⚫🔵 ➜ 🔴|⚫🔵 🔵|🔴⚫

🔹El objetivo de este enfoque es implementar dos clasificadores binarios, para luego combinar las predicciones. Se aplica el mismo proceso de discretización utilizado en el enfoque multiclase (ver Tabla 9), pero generando las categorías “Alto” y “Medio_Bajo” para un primer clasificador, “Bajo” y “Medio_Alto” para un segundo clasificador. Con esto se intenta capturar el orden jerárquico de las clases (Bajo < Medio < Alto), luego se unifican las predicciones de ambos modelos de forma tabular.

🔹Modelo A
- P(Bajo) ➡️ True
- P(Medio_Alto) ➡️ False
🔸Modelo B
- P(Alto) ➡️ True
- P(Medio_Bajo) ➡️ False

Table 13: Lgbm multiclase / Errores sobre datos de valdiación
Modelo.A	Modelo.B	Predicción.combinada
True (Baja)	False (Media_Baja)	Baja
True (Baja)	True (Alta)	No sabe
False (Media_Alta)	True (Alta)	Media
False (Media_Alta)	False (Media_Baja)	Alta

❗Luego de ajustar cada clasificador de manera individual, se combinan las predicciones siguiendo los criterios de la Tabla 13 . Para ambos clasificadores se utiliza una receta de preprocesamiento similar a la aplicada en el modelo multiclase, pero en este caso se utiliza la técnica de SMOTE, igualando la clase minoritaria a igual proporción que la mayoritaria, ya que se está comparando una categoría con dos categorías combinadas, produciendo un desbalance en los datos.

📎 La métrica por optimizar en el ajuste de hiperparametros también es la roc_auc, en este caso no se desea minimizar un tipo de error de clasificación, ya que para los clasificadores individuales es igual de importante clasificar correctamente los dos tipos de esquinas. Además el entrenamiento se hará con clases balanceadas, gracias a la aplicación del método SMOTE, el mismo será implementado a través del paquete 📦 themis¹³ .

📌 El ajuste de estos dos clasificadores sigue el mismo flujo de trabajo, a excepción de la discretización, por lo que pudiera escribirse una función para hacer ambos ajustes. Sin embargo, para mantener la legibilidad y transparencia de la publicación, se mostrará el ajuste individual de ambos modelos.

🔷 Modelo A 🔵|🔴

🔵P(Bajo) ➡️ True 🔴P(Medio_Alto) ➡️ False

Receta 🥘

Receta_class_A <- recipe(
  formula = delitos ~ .,
  data =  train
  ) %>%
  step_mutate(
    bancos_atm = bancos + atm,
    teatros_cine = teatros + cine,
    embajadas_cons = embajadas + consulados
    ) %>% 
  step_rm(bancos, atm, teatros, cine, id, pliegue,embajadas, consulados) %>% 
  step_nzv(all_predictors(), -c(estaciones_ferrocarril,bocas_de_subte,hotel_baja)) %>%
  step_mutate(
    delitos = factor(case_when(delitos <= 1 ~ "Bajo",
                               delitos >= 2 ~ "Medio_Alto"),
                     levels = c("Bajo","Medio_Alto"))
    ) %>% 
  step_smote(delitos, over_ratio = 1, seed = 1234)

❗Se utilizan pasos de preprocesamiento similares a las recetas anteriores, con la salvedad de incluir las variables de geolocalización, la nueva discretización de la variable respuesta y la inclusión del método de SMOTE.

Definción del modelo ✍️

  lightgbm_spec_A <- parsnip::boost_tree(
    mode = "classification",
    trees = tune(),
    mtry = tune(),
    min_n = tune(),
    tree_depth = tune(), 
    sample_size = tune(),
    loss_reduction = tune(), 
    learn_rate = tune()
  ) %>%
    set_engine("lightgbm", verbose=-1)
  
  # PARAMETROS
  lightgbm_params_A <- 
    lightgbm_spec_A %>% 
    dials::parameters() %>% 
    update(
      trees = trees(c(1000L, 1500L)),
      sample_size = sample_prop(c(.7, .99)),
      mtry = mtry(c(30L, 35))
    )
  
  lightgbm_wf_A <- workflow() %>%
    add_recipe(Receta_class_A) %>%
    add_model(lightgbm_spec_A)

Optimización de hiperparámetros 🕹️

  set.seed(123)
  doParallel::registerDoParallel(cores = parallel::detectCores() -1)
  
  lightgbm_res_A <- tune_bayes(lightgbm_wf_A, 
                             resamples = cv_folds,
                             iter = 10,
                             param_info = lightgbm_params_A,
                             metrics = metric_set(roc_auc),
                             control = ctrl, 
                             initial = 10)
  
  doParallel::stopImplicitCluster()

Show code

autoplot(lightgbm_res_A, type = "performance") +
  labs(title = "Evolución del error durante el ajuste de hiperpárametros") +
  theme_minimal()

Figure 22: Lgbm binario A / Evolución del error durante el ajuste de hiperpárametros

select_best(lightgbm_res_A) %>% knitr::kable() %>%  kableExtra::kable_styling(full_width = FALSE)

modelo_lightgbm_A <- finalize_workflow(x = lightgbm_wf_A,  
                                      parameters = select_best(lightgbm_res_A))

validacion_lgbm_A <- fit_resamples(
  object       = modelo_lightgbm_A,
  resamples    = cv_folds,
  metrics      = metric_set(roc_auc),
  control      = control_resamples(save_pred = FALSE)
)

final_lgbm_A <- last_fit(modelo_lightgbm_A, splits)

Show code

tabla <- 
  bind_rows(
  collect_metrics(validacion_lgbm_A) %>%
    dplyr::rename(valor = mean) %>% 
    mutate(Error = "validación cruzada"),
  collect_metrics(final_lgbm_A) %>%
    dplyr::rename(valor = .estimate) %>% 
    mutate(Error = "Diciembre 2019")) %>% 
    mutate(modelo = "Lgbm Binario A") %>% 
    dplyr::select(-.config) %>% 
    filter(.metric == "roc_auc")

knitr::kable(tabla,
  caption = "Lgbm binario A / Comparación de errores") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 14: Lgbm binario A / Comparación de errores
.metric	.estimator	valor	n	std_err	Error	modelo
roc_auc	binary	0.8201688	10	0.0023982	validación cruzada	Lgbm Binario A
roc_auc	binary	0.8128123	NA	NA	Diciembre 2019	Lgbm Binario A

Show code

grafico <-
collect_predictions(final_lgbm_A) %>% 
  roc_curve(delitos, .pred_Bajo) %>%
  autoplot() | 
  final_lgbm_A %>%
  collect_predictions() %>%
  mutate(across(c(.pred_class, delitos), factor, levels =  c("Bajo", "Medio_Alto"))) %>% 
  conf_mat(delitos, .pred_class) %>% 
  autoplot(type = "heatmap")

grafico

Figure 23: Lgbm binario A / Comparación de errores

El “Modelo A” busca discriminar la categoría “Bajo” del resto, sin embargo, en la Tabla 14 ,la métrica de error sufre una baja de rendimiento ⚠️ al comparar el error de validación cruzada con los datos de prueba, aunque es leve la diferencia, este es un síntoma de que el modelo no pudo generalizar de la mejor manera.

🔮 PREDICCIÓN SOBRE NUEVOS DATOS

Receta_class_A$steps[[4]] <- update(Receta_class_A$steps[[4]], skip = T)

model_fit_A <- update_recipe(modelo_lightgbm_A, Receta_class_A) %>% 
  fit(Data_set_modelos)

predict(model_fit_A, Data_set_modelos %>% head() %>% dplyr::select(-delitos))

# A tibble: 6 x 1
  .pred_class
  <fct>      
1 Medio_Alto 
2 Medio_Alto 
3 Medio_Alto 
4 Bajo       
5 Bajo       
6 Bajo

🔷 Modelo B 🔴|🔵

🔴 P(Alto) ➡️ True 🔵 P(Medio_Bajo) ➡️ False

Receta 🥘

Receta_class_B <- recipe(
  formula = delitos ~ .,
  data =  train
  ) %>%
  step_mutate(
    bancos_atm = bancos + atm,
    teatros_cine = teatros + cine,
    embajadas_cons = embajadas + consulados
    ) %>% 
  step_rm(bancos, atm, teatros, cine, id, pliegue,embajadas, consulados) %>% 
  step_nzv(all_predictors(), -c(estaciones_ferrocarril,bocas_de_subte,hotel_baja)) %>%
  step_mutate(
    delitos = factor(case_when(delitos <= 4 ~ "Media_Bajo",
                               delitos >= 5 ~ "Alto"),
                     levels = c("Alto","Media_Bajo"))
    ) %>% 
  step_smote(delitos, over_ratio = 1, seed = 1234)

Definción del modelo ✍️

lightgbm_spec_B <- parsnip::boost_tree(
    mode = "classification",
    trees = tune(),
    mtry = tune(),
    min_n = tune(),
    tree_depth = tune(), 
    sample_size = tune(),
    loss_reduction = tune(), 
    learn_rate = tune()
  ) %>%
    set_engine("lightgbm", verbose=-1)
  
  lightgbm_params_B <- 
    lightgbm_spec_B %>% 
    dials::parameters() %>% 
    update(
      trees = trees(c(1000L, 1500L)),
      sample_size = sample_prop(c(.7, .99)),
      mtry = mtry(c(30L, 35L))
    )
  
  lightgbm_wf_B <- workflow() %>%
    add_recipe(Receta_class_B) %>%
    add_model(lightgbm_spec_B)

Optimización de hiperparámetros 🕹️

set.seed(123)
doParallel::registerDoParallel(cores = parallel::detectCores() -1)

lightgbm_res_B <- tune_bayes(Receta_class_B, 
                           resamples = cv_folds,
                           iter = 10,
                           param_info = lightgbm_params_B,
                           metrics = metric_set(roc_auc),
                           control = ctrl,
                           initial = 10)

doParallel::stopImplicitCluster()

autoplot(lightgbm_res_B, type = "performance") +
  labs(title = "Evolución del error durante el ajuste de hiperpárametros") +
  theme_minimal()

Figure 24: Lgbm binario B / Evolución del error durante el ajuste de hiperpárametros

select_best(lightgbm_res_B) %>% knitr::kable() %>%  kableExtra::kable_styling(full_width = FALSE)

mtry	trees	min_n	tree_depth	learn_rate	loss_reduction	sample_size	.config
34	1461	8	14	0.0699555	26.46622	0.8966313	Model1

modelo_lightgbm_B = finalize_workflow(x = lightgbm_wf_B,  
                                      parameters = select_best(lightgbm_res_B))

validacion_lgbm_B  <- fit_resamples(
  object       = modelo_lightgbm_B,
  resamples    = cv_folds,
  metrics      = metric_set(roc_auc),
  control      = control_resamples(save_pred = FALSE)
)

final_lgbm_B <- last_fit(modelo_lightgbm_B, splits)

Show code

tabla <- 
  bind_rows(
  collect_metrics(validacion_lgbm_B) %>%
    dplyr::rename(valor = mean) %>% 
    mutate(Error = "validación cruzada"),
  collect_metrics(final_lgbm_B) %>%
    dplyr::rename(valor = .estimate) %>% 
    mutate(Error = "Diciembre 2019")) %>% 
    mutate(modelo = "Lgbm Binario B") %>% 
    dplyr::select(-.config) %>% 
    filter(.metric == "roc_auc")

knitr::kable(tabla,
  caption = "Lgbm binario B / Comparación de errores") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 15: Lgbm binario B / Comparación de errores
.metric	.estimator	valor	n	std_err	Error	modelo
roc_auc	binary	0.8981792	10	0.002967	validación cruzada	Lgbm Binario B
roc_auc	binary	0.9108769	NA	NA	Diciembre 2019	Lgbm Binario B

Show code

grafico <-
collect_predictions(final_lgbm_B) %>% 
  roc_curve(delitos, .pred_Alto) %>%
  autoplot() | 
  final_lgbm_B %>%
  collect_predictions() %>%
  mutate(across(c(.pred_class, delitos), factor, levels = c("Alto", "Media_Bajo"))) %>% 
  conf_mat(delitos, .pred_class) %>% 
  autoplot(type = "heatmap")

Figure 25: Lgbm binario B / Comparación de errores

En la Tabla 15, el “Modelo B” muestra el caso opuesto del “Modelo A,” mejorando su rendimiento, se recuerda que este modelo busca discriminar la categoría “Alto” del resto.

🔮 PREDICCIÓN SOBRE NUEVOS DATOS

Receta_class_B$steps[[4]] <- update(Receta_class_B$steps[[4]], skip = T)

model_fit_B <- update_recipe(modelo_lightgbm_B, Receta_class_B) %>% 
  fit(Data_set_modelos)

predict(model_fit_B, Data_set_modelos %>% head() %>% dplyr::select(-delitos))

# A tibble: 6 x 1
  .pred_class
  <fct>      
1 Alto       
2 Alto       
3 Media_Bajo 
4 Media_Bajo 
5 Media_Bajo 
6 Media_Bajo

Unificación de modelos 🧲

📝 Aplicando la metodología planteada en la Tabla 13 , se obtiene una combinación de las predicciones de ambos clasificadores, mostrando dicha combinación en una de las matrices de confusión que se muestran a continuación . Se observa una mejoría en la predicción de las categorías extremas por parte de la combinación de los clasificadores binarios 👍🏻, la incorporación del orden jerárquico de las categorías a través de este enfoque, arrojó resultados satisfactorios ,✅ sin embargo, se observó una menor precisión en la categoría media.

Show code

data_original <- Data_set_modelos %>% 
  mutate(
    delitos = case_when(delitos <= 1 ~ "Bajo",
                        delitos <= 4  ~ "Medio",
                        delitos >= 5 ~ "Alto"))

tabla_comparativa <- bind_cols(
  collect_predictions(final_lgbm_A) %>% 
    dplyr::rename(A = .pred_class) %>% 
    dplyr::select(A),
  collect_predictions(final_lgbm_B) %>% 
    dplyr::rename(B = .pred_class) %>% 
    dplyr::select(B),
  data_original %>% 
    filter(
      pliegue == 13) %>% 
    dplyr::select(delitos)
  ) %>% 
  mutate_all(as.character) %>% 
  mutate(B = gsub("Media","Medio", B)) %>% 
  mutate(across(c(A,B), ~ gsub("_"," ",.))) %>% 
  rowwise() %>% 
  mutate(Categoria_final = paste0(intersect(unlist(str_split(A, " ")),unlist(str_split(B, " "))),collapse = "")) %>% 
  mutate(across(c(Categoria_final, delitos),factor, levels =  c("Alto","Medio","Bajo")))

autoplot(Matriz_multiclase, type = "heatmap") | 
  autoplot(conf_mat(table(tabla_comparativa$Categoria_final,tabla_comparativa$delitos)), type = "heatmap")

Figure 26: Matrices de confusión de los modelos de clasificación

📌Analizando la problemática desde la perspectiva de una hipotética organización policial, se puede tolerar de mejor manera un error de predicción de un esquina considerada de riesgo medio y que se clasifique como de bajo o alto riesgo ✅.

👉🏻 Sin embargo, una esquina considerada de bajo riesgo, que reciba una clasificación de alto riesgo o viceversa, puede resultar un error sumamente costoso 🚫 desde el punto de vista logístico, social y económico. Aunque ambos enfoques presentan buenos rendimientos discriminando las categorías extremas, se propone asignar un mayor costo a este tipo de errores a través de una matriz de costos.

Matriz de costos ⚖

🔹 Se comparan ambos enfoques de clasificación asignando costos a los errores de las matrices.

Matriz_costos <- data.frame(Alto = c(-1,1,10),
                           Medio = c(1,0,1),
                           Bajo = c(10,1,-1))

rownames(Matriz_costos) = c("Predict -> Alto","Predict -> Medio","Predict -> Bajo")

knitr::kable(Matriz_costos, caption = "Matriz de costos") %>% 
  kableExtra::kable_styling(full_width = FALSE)

Table 16: Matriz de costos
	Alto	Medio	Bajo
Predict -> Alto	-1	1	10
Predict -> Medio	1	0	1
Predict -> Bajo	10	1	-1

🔹 A modo ilustrativo, en la 16 se muestra la matriz de costos, asignado los costos de manera arbitraria 🎲, con el objetivo de servir de guía para la mejora de las decisiones por parte de la organización policial hipotética. Para los errores extremos fue asignado un costo de 10 unidades, para el resto de los errores el costo fue de 1, los aciertos de los casos extremos generan un beneficio de 1 unidad y la predicción del riesgo medio no presenta costo.

Figure 27: Evaluación de costos

📌En la Figura 26 ,se observa que los clasificadores binarios presentan un costo más bajo para la organización, a excepción del par “Bajo-Medio”, el resto de los errores fueron sustancialmente menores que el clasificador multiclase. En definitiva, el enfoque de clasificadores binario superó al modelo de clasificación multiclase, sin embargo, 👉🏻👉🏻 esto puede variar en función de las prioridades que tenga la organización. Si la predicción correcta del riesgo medio se convierte en una prioridad, asignándole a la misma un valor positivo en la matriz de costos, el modelo multiclase pasaría a ser el de menor costo.

Deployment de modelos 🎬 📽

🔎 Simplemente como un ejemplo didáctico, se propone guardar 💾 los modelos entrenados a través de métodos y formatos que permitan utilizarlos en un futuro, de manera independiente al flujo de trabajo 🔄 utilizado, además de poder guardar distintas versiones de los modelos ajustados . Si bien existen muchos enfoques para guardar modelos, a continuación, se muestra el método RDS, y el método binario que permite la portabilidad entre los lenguajes Python y R ↔︎.

Muestra de datos para predecir 🔮

🔹 Para los cuatro modelos ajustados, se toma una muestra aleatoria de seis observaciones de los datos originales, eliminando además la variable delitos, con el objetivo de “simular” 🔢 la generación de predicciones aisladas del flujo del trabajo.

set.seed(123)
# Datos a predecir
Datos_sample <- Data_set_modelos %>% sample_n(6) %>% 
  dplyr::select(-delitos)

rmarkdown::paged_table(Datos_sample)

Método RDS 💾

🔹 Los cuatro modelos ajustados se guardan en formato “.rds” , posteriormente pueden cargarse en una nueva sesión independiente.

# Save y Load de modelos ajustados
# Regresion
saveRDS(modelo_fit_reg, "modelo_fit_reg.rds")
modelo_fit_reg = readRDS("modelo_fit_reg.rds")

# Multiclass
saveRDS(model_fit_multiclass, "model_fit_multiclass.rds")
model_fit_multiclass <- readRDS("model_fit_multiclass.rds")

# Binario A
saveRDS(model_fit_A, "model_fit_A.rds")
model_fit_A <- readRDS("model_fit_A.rds")

# Binario B
saveRDS(model_fit_B, "model_fit_B.rds")
model_fit_B <- readRDS("model_fit_B.rds")

🔹El objeto guardado es un “workflow” que contiene dos elementos, la receta de preprocesamiento y el modelo ajustado, estos se ejecutarán de forma secuencial al momento de predecir sobre nuevos datos. Una de las ventajas que ofrece tidymodels reside en lo siguiente, los datos a predecir que están almacenados dentro del objeto “Datos_sample” están en su formato bruto 🤔, es decir, la variable dependiente no está discretizada y no tiene ningún tipo de preprocesamiento o selección de variables.

📌Pero a través de las recetas de preprocesamiento guardadas dentro de los “workflows,” se pueden aplicar modelos que requieren transformaciones y distintas aplicaciones como regresión, clasificación multiclase y binaria, sin necesidad de modificar la data original o el flujo de trabajo 🙌🏻.

bind_cols(regreesion = predict(modelo_fit_reg, Datos_sample)$.pred,
          multiclase = predict(model_fit_multiclass, Datos_sample)$.pred_class,
          clasificdor_A = predict(model_fit_A, Datos_sample)$.pred_class,
          clasificdor_B = predict(model_fit_B, Datos_sample)$.pred_class) %>% 
  rmarkdown::paged_table()

🔹Se observan las predicciones de los 4 modelos para las mismas 6 observaciones almacenadas en su formato original.

Portabilidad de modelos R❤️ Python 🐍

👉🏻👉🏻 Este enfoque consiste en guardar los modelos como objetos binarios, muchos modelos tienen un método asociado que permite guardarlos en este formato ✅. El objeto generado tiene la ventaja de ser universal entre las diversas interfaces del modelo 🌐, se menciona como desventaja el hecho de tener que descomponer el “workflow” de tidymodels. A continuación, a modo de ejemplo se muestra la implementación para el modelo XGB.

Modelo <- extract_model(modelo_fit_reg)

Receta <- extract_recipe(modelo_fit_reg)

🔹 El primer paso consiste en utilizar las funciones extract_model y extract_recipe, con la finalidad de extraer el modelo ajustado y la receta de preprocesamiento del “workflow.”

xgboost::xgb.save(Modelo, 'modelo_xgb_deploy.model')

[1] TRUE

regresion_data <- bake(Receta, new_data = Datos_sample) %>% as.matrix()

Modelo_load <- xgboost::xgb.load('modelo_xgb_deploy.model')

predict(Modelo_load ,regresion_data)

[1] 1.7409922 2.3842442 2.8863437 5.2717023 0.9696367 4.2494283

🔹 Posteriormente se guarda el modelo como un objeto binario a través del método “save” y mediante la receta se aplica el preprocesamiento de los datos a predecir con la función bake 🍳. Finalmente se carga el modelo binario y se realiza la predicción 🔮 sobre los datos previamente preprocesados.

🔹 Se puede observar que se obtienen los mismos resultados mostrados en el método “.rds,” con la desventaja de no poder usar directamente el workflow 🔄, sin embargo, este formato binario permite la portabilidad del modelo con Python 🐍. A continuación, se exporta la data preprocesada por la receta en formato “.csv” para luego importarla en Python y realizar la predicción.

write.csv(regresion_data, 'regresion_data.csv', row.names=FALSE)

🔹El siguiente paso consiste en abrir una sesión en Python, cargar el objeto binario ‘modelo_xgb_deploy.model’ exportado desde R, además del archivo “.csv,” se transforman los datos en una matriz “sparse,” para luego realizar la predicción. A continuación, el procedimiento correspondiente en Python.

Código Python

!pip install xgboost==1.2.1
import xgboost
import pandas as pd

model = xgboost.Booster()
model.load_model('modelo_xgb_deploy.model')
data = pd.read_csv('regresion_data.csv')

data=xgboost.DMatrix(data.values)

model.predict(data)

array([1.7409934, 2.3842454, 2.8863425, 5.271702, 0.96963775, 4.249423],

dtype=float32)

🔹 Se observan los mismos resultados para las tres predicciones sobre el mismo conjunto de datos. Aunque el código Python se ejecutó en una instancia fuera de R, el mismo puede implementarse dentro de R a través del paquete 📦 reticulate,¹⁴ leyendo el dataframe de R como un pandas dataframe sin necesidad de exportar los datos en formato “.csv.”

🔹 El proceso inverso también es posible, entrenar un modelo en Python, exportarlo como binario y finalmente cargarlo en R para hacer predicciones.

💡 Comentarios finales

Para el enfoque de regresión se pueden aplicar modelos de conteo de eventos 📝, de esta manera se estaría tratando el conteo de delitos ocurridos de una manera más realista.
El algoritmo diseñado para realizar el cálculo que determina cuántos delitos ocurrieron en las cercanías de cada esquina, es excesivamente lento si se analiza un número muy grande de delitos. Si se lograra optimizar su implementación 🤔, se pudiera realizar primero el cálculo de delitos para todas las esquinas y posteriormente eliminar aquellas esquinas con bajos registros de delitos, dejando como última etapa la selección óptima de esquinas en base al filtro de proximidad. Por limitaciones computacionales 💻, el presente estudio tuvo que implementarse de manera inversa.
Probar con distintas ventanas de entrenamiento puede ser una alternativa interesante para validar los resultados 🗓, evaluando cuál sería la ventana de meses más óptima, Lin, Yen, and Yu (2018) en su implementación aplicaron esta estrategia.
Considerar la ocurrencia de delitos en esquinas cercanas 📍en los distintos meses de evaluación, puede ayudar a mejorar las predicciones incorporando estos insights a los modelos , Lin, Yen, and Yu (2018) también incorporó este factor en el estudio.
Probar distintos enfoques de discretización para ocurrencia de delitos , implementar múltiples algoritmos 🤖, ajustar un mayor número de hiperparámetros con una mayor cantidad de iteraciones, utilizar otras métricas 🎯 y funciones de pérdida , además de abordar el estudio con modelos sin supervisión, puede darle una perspectiva más profunda al análisis 💡.

🌐 Siéntase libre de comentar , sugerir o compartir cualquier idea.

Muchas gracias 👏🏻👏🏻👏🏻

Contacto ✉

Rafael Zambrano, Linkedin, Twitter, Github, Blogpost.

Chawla, Nitesh V., Kevin W. Bowyer, Lawrence O. Hall, and W. Philip Kegelmeyer. 2002. “SMOTE: Synthetic Minority over-Sampling Technique.” Journal of Artificial Intelligence Research 16: 321–57.

Chen, Tianqi, and Carlos Guestrin. 2016. “KDD ’16: The 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining.” In. ACM. https://doi.org/10.1145/2939672.2939785.

Cheng, Joe, Bhaskar Karambelkar, and Yihui Xie. 2019. Leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library. https://CRAN.R-project.org/package=leaflet.

Christou, Vasiliki, and Konstantinos Fokianos. 2013. “On Count Time Series Prediction.” Journal of Statistical Computation and Simulation 85 (2): 357–73. https://doi.org/10.1080/00949655.2013.823612.

Coene, John. 2020. Echarts4r: Create Interactive Graphs with ’Echarts JavaScript’ Version 4. https://CRAN.R-project.org/package=echarts4r.

Dancho, Matt, and Davis Vaughan. 2020. Timetk: A Tool Kit for Working with Time Series in r. https://CRAN.R-project.org/package=timetk.

Falbel, Daniel, Athos Damiani, and Roel M. Hogervorst. 2020. Treesnip: Tree Based Model Engines for ’Parsnip’. https://github.com/curso-r/treesnip.

Hand, David J., and Robert J. Till. 2001. Machine Learning 45 (2): 171–86. https://doi.org/10.1023/a:1010920819831.

Hilbe, Joseph M. 2017. “The Statistical Analysis of Count Data / El análisis Estadístico de Los Datos de Recuento.” Cultura y Educación 29 (3): 409–60. https://doi.org/10.1080/11356405.2017.1368162.

Hvitfeldt, Emil. 2020. Themis: Extra Recipes Steps for Dealing with Unbalanced Data. https://CRAN.R-project.org/package=themis.

Hyndman, & Athanasopoulos, R. J. 2018. “Forecasting: Principles and Practice.” Melbourne, Australia: [Ebrary version]. [Online]; OTexts. https://otexts.com/fpp2/.

Ke, Guolin, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. 2017. “LightGBM: A Highly Efficient Gradient Boosting Decision Tree.” In Proceedings of the 31st International Conference on Neural Information Processing Systems, 3149–57. NIPS’17. Red Hook, NY, USA: Curran Associates Inc.

Ke, Guolin, Damien Soukhavong, James Lamb, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. 2020. Lightgbm: Light Gradient Boosting Machine. https://CRAN.R-project.org/package=lightgbm.

Kourentzes, Nikolaos, and Fotios Petropoulos. 2016. Tsintermittent: Intermittent Time Series Forecasting. https://CRAN.R-project.org/package=tsintermittent.

Kuhn, Max, and Davis Vaughan. 2020. Yardstick: Tidy Characterizations of Model Performance. https://CRAN.R-project.org/package=yardstick.

Kuhn, Max, and Hadley Wickham. 2020a. Recipes: Preprocessing Tools to Create Design Matrices. https://CRAN.R-project.org/package=recipes.

———. 2020b. Tidymodels: A Collection of Packages for Modeling and Machine Learning Using Tidyverse Principles. https://www.tidymodels.org.

Lin, Ying-Lung, Meng-Feng Yen, and Liang-Chih Yu. 2018. “Grid-Based Crime Prediction Using Geographical Features.” ISPRS International Journal of Geo-Information 7 (8): 298. https://doi.org/10.3390/ijgi7080298.

OpenStreetMap contributors. 2017. “Planet dump retrieved from https://planet.osm.org .” https://www.openstreetmap.org.

Padgham, Mark, Bob Rudis, Robin Lovelace, and Maëlle Salmon. 2020. Osmdata: Import OpenStreetMap Data as Simple Features or Spatial Objects. https://CRAN.R-project.org/package=osmdata.

Pedersen, Thomas Lin, and David Robinson. 2020. Gganimate: A Grammar of Animated Graphics. https://CRAN.R-project.org/package=gganimate.

Rabosky, AR Davis, CL Cox, DL Rabosky, PO Title, IA Holmes, A Feldman, and JA McGuire. 2016. “Coral Snakes Predict the Evolution of Mimicry Across New World Snakes.” Nature Communications 7: 11484.

Spaulding, Jamie, and Keith Morris. 2020. rcrimeanalysis: An Implementation of Crime Analysis Methods. https://github.com/JSSpaulding/rcrimeanalysis.

Ushey, Kevin, JJ Allaire, and Yuan Tang. 2020. Reticulate: Interface to ’Python’. https://github.com/rstudio/reticulate.

Zambrano, Rafael. 2021. “Rafael Zambrano: Evaluación de Predicciones Mediante Programación Funcional En r.” https://rafael-zambrano-blog-ds.netlify.app/posts/2021-01-03-multieval/.

Zambrano, Rafael, and Karina Bartolome. 2021. Sknifedatar: Swiss Knife of Data for r. https://github.com/rafzamb/sknifedatar.

Página de La Policía de la Ciudad de Buenos Aires https://www.policiadelaciudad.gob.ar/↩︎
Zambrano and Bartolome (2021)↩︎
Repositorio público del GCABA https://data.buenosaires.gob.ar/↩︎
Cheng, Karambelkar, and Xie (2019)↩︎
Spaulding and Morris (2020)↩︎
Padgham et al. (2020)↩︎
Rabosky et al. (2016)↩︎
Kourentzes and Petropoulos (2016)↩︎
Kuhn and Vaughan (2020)↩︎
Kuhn and Wickham (2020a)↩︎
Falbel, Damiani, and Hogervorst (2020)↩︎
Ke et al. (2020)↩︎
Hvitfeldt (2020)↩︎
Ushey, Allaire, and Tang (2020)↩︎

Comment on this article Share:

Predicción de delitos en CABA

⚠️ Antes de empezar …

Resumen 📖

Principales bibliotecas 📚

Importación de datos 🔌

Limpieza de datos 🧹

Construcción del Dataset 🛠️

Tipos de delitos 🚨

Análisis exploratorio de datos 🔭

Selección de esquinas de CABA 📍

Intercepción de calles y avenidas 🛣️

Selección óptima de esquinas en función a la proximidad ⚛️

Cálculo mensual de los delitos ocurridos en cada esquina de CABA 🚔

Indexación de variables meteorológicas 🌤️≫🚔

Depuración de esquinas con baja incidencia delictiva 👁️

Indexación de variables de entorno 🌏≫🚔

Distancias entre esquinas y factores de entorno físico 📐

Aplicación del “método” de ventanas deslizantes 🤔🧐

Modelado 📚💻🎯

Baseline: Series de tiempo 📉

🔹 Modelos de series de tiempo

Evaluación de Errores 🎯

Modelos de aprendizaje automático 🚀

Partición de datos ➗

Validación cruzada 🧮

Enfoque de regresión 📈

Receta de preprocesamiento 🥘

Definición del modelo ✍️

Optimización de hiperparámetros 🕹️

Comparación de errores 🔎

Comparación con el Baseline 🏁

Ajuste final del modelo ☑️

Enfoque de clasificación 📊

Clasificación multiclase 🔴⚫🔵

Receta 🥘

Especificación del modelo ✍️

Optimización de hiperparámetros 🕹️

Comparación de errores 🔎

Predicción sobre nuevos datos 🔮

Clasificación binaria múltiple

🔴⚫🔵 ➜ 🔴|⚫🔵 🔵|🔴⚫

🔷 Modelo A 🔵|🔴

🔷 Modelo B 🔴|🔵

Unificación de modelos 🧲

Matriz de costos ⚖

Deployment de modelos 🎬 📽

Muestra de datos para predecir 🔮

Método RDS 💾

Portabilidad de modelos R❤️ Python 🐍

💡 Comentarios finales

Contacto ✉

References

Reuse

Citation