rm(list=ls())

Resumen

Que los robos están aumentando en el Valle de Aburrá es una verdad que ya desbordó las cifras. Lo que no deja de sorprender es la manera en la que muchos de estos crímenes ocurren, a plena luz del día, en vías principales o en establecimientos públicos vigilados. ¿A qué se debe que los ladrones estén asumiendo tantos riesgos a la hora de actuar?.

El gran objetivo de este trabajo, es analizar el comportamiento de las geolocalizaciones de robos realizados en la ciudad de Medellín - Colombia, con el fin de detectar las zonas más vulnerables a partir de modelos estadísticos de procesos puntuales

  1. Introducción

Es bastante evidente que la ciudad de Medellín es fuertemente vulnerada por grupos delictivos, en consecuencia, se desea tratar el tema de los robos en la ciudad, que se caracterizan por múltiples modalidades.

En general, la seguridad de la ciudad de Medellín es un tema bastante complejo, sin embargo, la estadística permite analizar comportamientos en este tipo de actos que atentan contra la sociedad, siendo más óptimo al aplicar estadística espacial, con el fin de relacionar estos hechos a partir de datos georeferenciados, para así, poder contribuir de forma significativa en posibles medidas que puedan mitigar estas situaciones.

  1. Descubrimiento y entendimiento de datos

Este conjunto de datos presenta las locaciones de puntos donde ocurrieron robos en diferentes zonas de la ciudad de Medellín, dichos registros, fueron tomados por la policía nacional.

los datos fueron extraídos de la página de MEData, en la cuál yacen datos abiertos suministrados por la alcaldía de Medellín: Hurtos_Medellin.

Cabe resaltar, que la base de datos original cuenta con un total de 275.000 registros, por la magnitud de la misma, se decide depurar los datos, filtrándo únicamente los robos del año 2018, en vía pública y por modalidad cosquilleo, obteniendo así, un total de 1439 registros en la ciudad de Medellín.

Las variables que se tienen en cuenta para el estudio son las siguientes:

Variable Descripción Tipo Detalle
X Número del registro interger 1-1439
fecha_hecho Fecha en la cual se realizó el robo datetime -
cantidad Cantidad robos -
latitud Coordenada latitud de la ubicación del robo float -
longitud Coordenada longitud de la ubicación del robo float -
sexo sexo víctima string -
edad Edad víctima integer -
estado_civil Estado civíl víctima string -
grupo_actor grupo al que pertenece String -
actividad_delictiva historial víctima String -
parentesco parentesco víctima string -
ocupación A que se dedica la víctima string -
discapacidad discapacidad víctima string -
grupo_especial Grupo especial de la víctima string -
medio_transporte Medio transporte usado por la víctima string Clases: metro, caminata, autobús, automóvil, taxi, motocicleta
nivel_académico Nivel académico de la víctima string -
testigo testigo de la víctima string -
conducta conducta víctima string -
modalidad modalidad de robo integer Clases: cosquilleo, atraco, decuido, escopolamina, agentes químicos, abuso de confianza, clonación de tarjeta, comisión de delito, forcejeo, informático, llamada millonaria, llave maestra, paquete chileno, retención de tarjeta, simulando necesidad, suplantación
caracterizacion caracterización víctima string -
conducta_especial - - -
arma_medio Arma utilizada en el robo string corto punzante, arma de fuego, palanca, escopolamina, no
articulo_penal - - -
categoria_penal - - -
nombre_barrio Barrio donde ocurrió el hecho string -
codigo_barrio Código del Barrio donde ocurrió el hecho interger -
codigo_comuna Código de comuna donde ocurrió el hecho interger -
lugar Lugar donde ocurrió el hecho string -
sede_receptora Sede policia donde se registró el hecho string -
bien Artículo robado string -
categoría_bien Categoría de artículo robado string -
grupo_bien Grupo de artículo robado string -
modelo - - -
color - - -
permiso - - -
unidad_medida - - -
fecha_ingestion - - - 
#Lectura librerías
library(leaflet)
library(dplyr)
library(sf)
library(rgdal) 
library(raster)
library(spatstat)
library(rgeos)
library(mapview)
library(maptools)
library(ggplot2)
library(kableExtra)

Lectura Base De Datos 

#Lectura Datos
crimen.med <- read.csv("crimen_med.csv", header = TRUE,
                  encoding="UTF-8")
head(crimen.med)

Lectura shapefile

#lectura de .csv y .shp
catastral <- read.csv("Limite_Barrio_Vereda_Catastral.csv", encoding="UTF-8")
med.shp <- st_read("Limite_Barrio_Vereda_Catastral.shp")
## Reading layer `Limite_Barrio_Vereda_Catastral' from data source 
##   `D:\Desktop\Estadística\EstadisticaEspacial\Robos_Medallo\Limite_Barrio_Vereda_Catastral.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 371 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.71931 ymin: 6.162904 xmax: -75.47185 ymax: 6.374872
## Geodetic CRS:  WGS 84
barrio_vereda <- read.csv("Barrio_Vereda_2014.csv", encoding="UTF-8")

Nombres de columnas shapefile

colnames(med.shp)
## [1] "OBJECTID"   "SHAPEAREA"  "SHAPELEN"   "COMUNA"     "BARRIO"    
## [6] "CODIGO"     "NOMBRE_BAR" "NOMBRE_COM" "geometry"

Mapa y locaciones 

mapa_med <- leaflet() %>%
  addProviderTiles(providers$OpenStreetMap)
mapa_med %>% 
  addMarkers(lat=~latitud, lng =~longitud, data = crimen.med,
                   icon = makeIcon(iconUrl = "robber.png",
                            iconWidth = 10, iconHeight = 10),
             popup = paste(
               "<h3>Datos del robo</h3>",
               
               paste("Edad:",
                     (crimen.med$edad)),
               
               paste("Arma de por medio:",
                     (crimen.med$arma_medio)),
               sep = "<br>"
             )
  )

Transformación coordenadas a UTM

get_utm <- function(longitud, latitud, zona, loc){
  points = SpatialPoints(cbind(longitud, latitud), proj4string = CRS("+proj=longlat +datum=WGS84"))
  points_utm = spTransform(points, CRS(paste0("+proj=utm +north=T +zone=",zona[1]," +ellps=WGS84")))
  if (loc == "longitud") {
    return(coordinates(points_utm)[,1])
  } else if (loc == "latitud") {
    return(coordinates(points_utm)[,2])
  }
}

crimen.med %<>% 
  mutate(zona = (floor((longitud + 180)/6) %% 60) + 1, keep = "all"
  ) %>% 
  group_by(zona) %>% 
  mutate(utm_x = get_utm(longitud, latitud, zona, loc = "longitud"),
         utm_y = get_utm(longitud, latitud, zona, loc = "latitud"))

Puntos en el plano

plot(crimen.med$utm_x,crimen.med$utm_y)

Geometría Medellín longlat

plot(med.shp$geometry)

Transformado coordenadas de Shapefile a UTM

medshputm <- st_transform(med.shp, "+proj=utm +zone=18 ellps=WGS84")

Geometría Medellín UTM

plot(medshputm$geometry)

Se toman las especificaciones del cuadro delimitador de la variable medshputm para establecer la extensión geográfica para el proceso de generación de datos.

xmin <- as.vector(st_bbox(medshputm)[1])
ymin <- as.vector(st_bbox(medshputm)[2])
xmax <- as.vector(st_bbox(medshputm)[3])
ymax <- as.vector(st_bbox(medshputm)[4])

Análisis descriptivo 

#Mapa para todos los barrios, usando 'innerjoin' con el .shp de Limite_Barrio_Vereda_Catastral
Unido <- inner_join(catastral, crimen.med, by = c("COMUNA" = "codigo_comuna"))
nueva_base <- Unido %>% filter(YEAR >= 2018 & YEAR <= 2018) %>% 
  group_by(CODIGO) %>%
  dplyr::summarise(robos = n()) %>%
  dplyr::ungroup()
#Se realiza la conversion de CODIGO a formato numerico
med.shp$CODIGO <- as.numeric(as.character(med.shp$CODIGO))
#Se utiliza 'inner join' para unir dos bases y para luego generar mapa
mapa <- inner_join(med.shp, nueva_base, by = c("CODIGO" = "CODIGO"))
mypal <- colorNumeric(palette = c("#000000","#280100","#3D0201","#630201","#890100","#B00100","#DD0100","#F50201",
                                   "#FF5F5E","#FF7A79","#FF9796","#FEB1B0","#FDC9C8", "#FFE5E4"), domain = mapa$accidentes, reverse = T)
# Creaci?n del mapa
leaflet() %>% addPolygons(data = mapa, color = "#0A0A0A", opacity = 0.6, weight = 1, fillColor = ~mypal(mapa$robos),
                          fillOpacity = 0.6, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "black", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", mapa$NOMBRE_BAR, "<br>", "Robos: ", mapa$robos, "<br>")) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", pal = mypal, values = mapa$robos, title = "Robos", opacity = 0.6)
medellin_comuna <- crimen.med %>% filter(YEAR >= 2018 & YEAR <= 2018) %>% 
  group_by(sede_receptora) %>% 
  dplyr::summarize(robos = n())
ggplot(data = medellin_comuna, aes(x = reorder(sede_receptora,+robos), y = robos)) +
  geom_bar(stat = "identity",  fill = "red", color = "black", alpha = 0.6) +
  geom_text(aes(y = robos, label = robos),
                position = position_dodge(width = 0.7), size = 3.5, vjust = 0.5, hjust = -0.1, col = "black") +
  xlab("Comuna") + 
  ylab("Total de Robos") +
  ggtitle("Total de Robos por Comuna") +
  ylim(c(0,1000)) +
  theme_minimal() +
  coord_flip()

Se puede evidenciar la vulnerabilidad de la comuna La Candelaria, la cual corresponde más específicamente al centro de Medellín.

kable(table(crimen.med$sexo),col.names = NULL) %>% 
  kable_styling(full_width = F) %>% 
  kable_minimal()
Hombre 464
Mujer 975 
kable(table(crimen.med$categoria_bien),col.names = NULL) %>% 
  kable_styling(full_width = F) %>% 
  kable_minimal()
Alimento 51
Dinero, joyas, piedras preciosas y título valor 302
Documentos 283
Electrodomésticos 3
Herramientas 1
Librería, Papelería y útiles 3
Otros elementos 5
Prendas de vestir y accesorios 153
Sin dato 1
Tecnología 637

Objeto ppp

Finalmente, construimos objetos ppp a partir de nuestros tres conjuntos de datos diferentes, crimen.med, medellin.regular y medellin.random usando la función ppp().

class(crimen.med)
## [1] "grouped_df" "tbl_df"     "tbl"        "data.frame"

Coordenadas UTM

crimen.med[, 41:42]

Ahora, construyamos un objeto ppp, denominado ppp.locations.

crimen.med$utm_x <- as.numeric(crimen.med$utm_x)
crimen.med$utm_y <- as.numeric(crimen.med$utm_y)
win <- as(as(medshputm, 'Spatial'), 'owin')
coordinate.units <- c("metre", "metres")
ppp.locations <- ppp(x = crimen.med$utm_x,
                     y = crimen.med$utm_y,
                     window = win)

Ahora asignamos el nombre de la unidad de longitud para las coordenadas espaciales al objeto ppp y lo trazamos

unitname(ppp.locations) <- coordinate.units
plot(ppp.locations)

  1. Análisis de patrones de puntos espaciales (Intensidad)

La intensidad es la densidad promedio de puntos o, en otras palabras, el número esperado de puntos por unidad de área. La intensidad puede ser constante (uniforme) o puede variar de un lugar a otro (no uniforme o no homogéneo).

Un enfoque para evaluar la intensidad es dividir el área de estudio en cuadrantes y contar el número de puntos que caen en cada cuadrante. Si los puntos tienen una intensidad uniforme y son completamente aleatorios, entonces los recuentos de cuadrantes deben ser números aleatorios de Poisson con media constante. Podemos probar esa hipótesis usando el estadístico de prueba de bondad de ajuste χ2

En R hacemos uso de las funciones quadratcount() y quadrat.test() del paquete spatstat.

qc.loc <- quadratcount(ppp.locations, nx=3, ny=4)
plot(ppp.locations, pch=3, cex=0.6)
plot(qc.loc, add=T, textargs = list(col='red'))

Prueba de chi-cuadrado de CSR utilizando recuentos de cuadrantes

quadrat.test(qc.loc)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 5407.4, df = 9, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 10 tiles (irregular windows)

El valor p muy bajo indica que la intensidad no es constante. Dicho proceso puntual se conoce como IPP y su intensidad se puede estimar de forma no paramétrica, por ejemplo, con suavizado del núcleo (Bivand et al. 2008).

Repitamos este enfoque para el conjunto de datos ppp.random.

Kernel density plot (2-D)

Un diagrama de densidad del kernel es la representación gráfica del suavizado del kernel. En R y paquetes contribuidos, hay varias funciones disponibles para la visualización de gráficos de densidad del kernel. Aquí usamos la función de density.ppp() del paquete spatstat.

La density.ppp() proporciona una estimación de la función de intensidad del proceso de puntos que generó los datos del patrón de puntos. El argumento sigma se toma como la desviación estándar del kernel gaussiano isotrópico, correspondiente al parámetro de ancho de banda. Una desviación estándar pequeña da como resultado un gráfico de densidad más puntiagudo y una desviación estándar grande da como resultado un gráfico de densidad más uniforme.

Codificamos un bucle para aplicar la función de density.ppp() en nuestros tres conjuntos de datos, ppp.locations, ppp.regular y ppp.random con diferentes valores de árbol para sigma (1000,2500,3000).

par(mfrow=c(1,3))
sigma <- c(1000, 2500, 5000)
main <- c('Localizaciones')
for (j in 1:3){
    ds <- density.ppp(ppp.locations, sigma=sigma[j])
    plot(ds, 
         main = paste0('Localizaciones', ', sigma: ', sigma[j]))
    plot(ppp.locations, add=T, cex=0.01, regular=F)
    }

Aplicando diferentes anchos de banda

Ahora, interesa aplicar diferentes anchos de banda, con el fin de analizar sus respectivos comportamientos, para finalmente concluir cuál o cuales serían los más indicados para éste contexto en Medellín.

cvldensity <- density(ppp.locations,sigma=bw.CvL(ppp.locations))
plot(cvldensity,main="CVL")

bw.CvL(ppp.locations)
##    sigma 
## 7407.238
scottdensity <- density(ppp.locations,sigma=bw.scott(ppp.locations))
plot(scottdensity,main="scott")

bw.scott(ppp.locations)
##  sigma.x  sigma.y 
## 402.3170 608.0769
fracdensity <- density(ppp.locations,sigma=bw.frac(ppp.locations))
plot(fracdensity,main="frac")

bw.frac(ppp.locations)
## [1] 6559.358
ppldensity <- density(ppp.locations,sigma=bw.ppl(ppp.locations))
plot(ppldensity,main="ppl")

bw.ppl(ppp.locations)
##    sigma 
## 327.2871