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)
library(spatstat.geom)

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 
##   `C:\Users\jdbul\OneDrive\Escritorio\UNIVERSIDAD\ASIG.VISTAS\Estadistica Espacial\Trabajo2-EstEspacial\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")
med.shp2 <- st_read("comunasmedellin.shp")
## Reading layer `comunasmedellin' from data source 
##   `C:\Users\jdbul\OneDrive\Escritorio\UNIVERSIDAD\ASIG.VISTAS\Estadistica Espacial\Trabajo2-EstEspacial\comunasmedellin.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 17 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -75.63243 ymin: 6.175719 xmax: -75.52478 ymax: 6.312849
## Geodetic CRS:  WGS 84

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 = 8, iconHeight = 8),
             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.

Género más vulnerable a robos en Medellín por cosquilleo

kable(table(crimen.med$sexo), col.names = NULL) %>%
  kable_paper("hover", full_width = F) %>%
  column_spec(2, bold = T,color = "#D7261E")
Hombre 464
Mujer 975

Categoría de artículos más robados por modalidad cosquilleo en Medellín

table(crimen.med$categoria_bien) %>% as.data.frame()%>%
  arrange(desc(Freq)) %>% kbl(col.names = NULL) %>%
  kable_paper("hover", full_width = F) %>%
  column_spec(2, bold = T,color = "#D7261E")
Tecnología 637
Dinero, joyas, piedras preciosas y título valor 302
Documentos 283
Prendas de vestir y accesorios 153
Alimento 51
Otros elementos 5
Electrodomésticos 3
Librería, Papelería y útiles 3
Herramientas 1
Sin dato 1

*Polígono de estudio estimado con función **ripras*

med.ripras <- ripras(crimen.med$utm_x, crimen.med$utm_y)
plot(med.ripras)

Shape área metropolitana de Medellín

medshputm2 <- st_transform(med.shp2, "+proj=utm +zone=18 ellps=WGS84")
plot(medshputm2$geometry)

Puntos aleatorios

Primero generamos números uniformes aleatorios usando la función runif().

set.seed(1111) # set seed for reproducibility
## generate random uniform numbers
x.random <- runif(n = nrow(crimen.med)*0.4, 
                  min = xmin, 
                  max = xmax)
y.random <- runif(n = nrow(crimen.med)*0.4, 
                  min=ymin, 
                  max=ymax)

Luego combinamos los datos del vector de números aleatorios en un objeto data.frame y lo trazamos junto con la geometría de Medellín.

## add to dataframe
medellin.random  <- data.frame('x' = x.random, 'y' = y.random)

plot(x=medellin.random$x, 
     y=medellin.random$y, 
     main="Random Points", 
     xlab="x", ylab="y", cex=.5, 
     asp=T)
plot(medshputm2$geometry, add=T, border='blue', col=NA)

Puntos regulares

Primero se genera una secuencia de números usando la función seq().

x.regular <- seq(from=xmin, to=xmax,
                 length.out = sqrt(nrow(crimen.med))*0.35)
y.regular <- seq(from=ymin,to=ymax,
                 length=sqrt(nrow(crimen.med))*0.35)

Luego se crea un marco de datos a partir de todas las combinaciones de variables usando la función expand.grid().

regular <- expand.grid(x.regular, y.regular)
plot(regular, asp=T)

Finalmente, agregamos algo de ruido a los puntos usando la función jitter(), combinamos los datos en un objeto data.frame, medellin.regular, y lo trazamos junto con la geometría de medellin.

x.regular <- jitter(regular[,1], factor = 1)
y.regular <- jitter(regular[,2], factor = 1)
medellin.regular <- data.frame('x' = x.regular, 'y' = y.regular)

plot(x=medellin.regular$x,
     y=medellin.regular$y, 
     main="Regular Points", 
     xlab="x",ylab="y", cex=.5, 
     asp=T)
plot(medshputm2$geometry, add=T, border='blue', col=NA)

Creación de objetos spatstat

Comenzamos con la asignación de una variable de marcador de posición que contiene la representación de cadena de las unidades de coordenadas.

coordinate.units <- c("metre", "metres")

Ventana

Un objeto de clase owin (ventana de observación) representa una región o ventana en un espacio bidimensional. La ventana puede ser un rectángulo, un polígono o polígonos, con agujeros poligonales o una forma irregular representada por una máscara de imagen de píxeles binarios.

class(medshputm2)
## [1] "sf"         "data.frame"

Para establecer la ventana en el polígono medshputm2, primero debemos convertir medshputm2 en un objeto espacial, que luego se convierte en un objeto owin. Esto puede sonar complicado, pero al usar la función as() en secuencia, el objetivo se vuelve bastante simple.

class(as(medshputm2, 'Spatial'))
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
class(as(as(medshputm2, 'Spatial'), 'owin'))
## [1] "owin"

Ahora asignamos el objeto owin a una variable denotada como win.

win <- as(as(medshputm2, 'Spatial'), 'owin')
win
## window: polygonal boundary
## enclosing rectangle: [430040.7, 441947.5] x [682661.9, 697821.7] units

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"
class(medellin.regular)
## [1] "data.frame"
class(medellin.random)
## [1] "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)
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)

Evidentemente, se aprecian observaciones fuera del área metropolitana de Medellín. Si queremos deshacernos de los puntos rechazados, debemos aplicar un pequeño truco: aplicar la función as.ppp() en el objeto ppp elimina los puntos rechazados.

ppp.locations <- as.ppp(ppp.locations)
plot(ppp.locations)

Continuamos con la construcción de un objeto ppp para el conjunto de datos medellin.random y lo denotamos como ppp.random. Como medellin.random es un objeto data.frame, simplemente indexamos sus columnas con la notación $.

ppp.random <- ppp(x = medellin.random$x,
                  y = medellin.random$y,
                  window = win)
unitname(ppp.random) <- coordinate.units
ppp.random <- as.ppp(ppp.random)
plot(ppp.random)

Finalmente, construimos un objeto ppp a partir del conjunto de datos medellin.regular y lo denotamos como ppp.regular.

ppp.regular <- ppp(x = medellin.regular$x,
                   y = medellin.regular$y,
                   window = win)
unitname(ppp.regular) <- coordinate.units
ppp.regular  <- as.ppp(ppp.regular)
plot(ppp.regular)

save(file = "ppp_data_medellin.RData", 
     ppp.locations, ppp.random, ppp.regular, crimen.med)

  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 = 3036.5, df = 11, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 12 tiles (irregular windows)

El valor p muy bajo indica que se rechaza la hipótesis de aleatoriedad completa. 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.

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

quadrat.test(qc.rand)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 18.243, df = 11, p-value = 0.1522
## alternative hypothesis: two.sided
## 
## Quadrats: 12 tiles (irregular windows)

El valor p alto (p>0.05) indica un HPP, donde no se rechaza la hipótesis de aleatoriedad completa, lo cual era de esperarse por la naturaleza de los puntos simulados.

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(3,3), mar=c(0,0,1,2))
sigma <- c(1000, 2500, 5000)
data <- list(ppp.locations, ppp.regular, ppp.random)
main <- c('Locations', 'Regular', 'Random')
for (i in 1:3){
  for (j in 1:3){
    ds <- density.ppp(data[[i]], sigma=sigma[j])
    plot(ds, 
         main = paste0(main[i], ', sigma: ', sigma[j]))
    plot(data[[i]], 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 
## 3656.692
scottdensity <- density(ppp.locations,sigma=bw.scott(ppp.locations))
plot(scottdensity,main="scott")

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

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

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

  1. Métodos de probabilidad para modelar la intensidad

Un modelo de Poisson homogéneo, p. λθ((x,y))=exp(θ0) puede ser ajustado por:

base.loc <- ppm(ppp.locations, ~1)
base.loc
## Stationary Poisson process
## Intensity: 1.417378e-05
##              Estimate       S.E.   CI95.lo   CI95.hi Ztest      Zval
## log(lambda) -11.16412 0.02644429 -11.21595 -11.11229   *** -422.1749

El modelo log-lineal

Un modelo de Poisson no homogéneo con una intensidad logarítmica lineal en las coordenadas cartesianas, p. λθ((x,y))=exp(θ0+θ1x+θ2y) se ajusta por:

loglin.loc <- ppm(ppp.locations, ~x+y)
coef(loglin.loc)
##   (Intercept)             x             y 
## -1.038934e+01  5.741105e-05 -3.740594e-05

Este modelo puede escribirse específicamente como:

λθ((x,y))=exp(−1.038934e+01+5.741105e-05x−3.740594e-05y)

El modelo log-cuadrático

Un modelo de Poisson no homogéneo con una intensidad que es log-cuadrática, tal que es cuadrática en x e y; p.ej.

λθ((x,y))=exp(θ0+θ1x+θ2y+θ3x2+θ4y2) se ajusta mediante

logquad.loc <- ppm(ppp.locations, ~polynom(x, y, 2))
loglin.rand <- ppm(ppp.random, ~x+y)
logquad.rand <- ppm(ppp.random, ~polynom(x, y, 2))

Predicción del modelo

Ahora que ajustamos un montón de modelos en dos conjuntos de datos diferentes, creamos estos modelos para la predicción.

El valor devuelto por la función de ajuste de modelos ppm es un objeto de clase ppm que representa el modelo ajustado. Esto es análogo a otros objetos modelo en R y, por lo tanto, se pueden aplicar operaciones estándar, como trazar, predecir o coef, entre otras.

Combinamos el comando plot() con el comando predict() para visualizar los modelos log-lineal y log-cuadrático para el conjunto de datos ppp.locations. Agregamos las coordenadas del punto para mejorar la visualización.

par(mfrow=c(1,2),mar=c(0,0,1,2))
plot(predict(loglin.loc), main='log-linear')
points(ppp.locations, pch=16, cex=0.2)
plot(predict(logquad.loc), main='log-quadratic')
points(ppp.locations, pch=16, cex=0.2)

Comprobación de un modelo de Poisson ajustado

Después de ajustar un modelo de proceso de punto a un conjunto de datos de patrón de punto, debemos verificar que el modelo se ajuste bien y que cada suposición de componente del modelo sea apropiada.

Un enfoque informal, no probabilístico, consiste en examinar los residuos del modelo.

par(mfrow=c(1,2),mar=c(1,1,2,2))
diagnose.ppm(loglin.loc, which = "smooth")
## Model diagnostics (raw residuals)
## Diagnostics available:
##  smoothed residual field
## range of smoothed field =  [-1.652e-05, 4.727e-05]
diagnose.ppm(logquad.loc, which = "smooth")

## Model diagnostics (raw residuals)
## Diagnostics available:
##  smoothed residual field
## range of smoothed field =  [-9.812e-06, 1.787e-05]

Creación objeto ppp.locations2 con ventana entregada por función ripras()

ppp.locations2 <- ppp(x = crimen.med$utm_x,
                     y = crimen.med$utm_y,
                     window = med.ripras)

Kest y pcf 

k <- Kest(ppp.locations2)
pcf <- pcf(ppp.locations2)
pcfinh <- pcfinhom(ppp.locations2)
kinh <- Kinhom(ppp.locations2)

par(mfrow=c(2,2))
plot(k, main = "K homogéneo")
plot(pcf, main = "PCF homogéneo")
plot(pcfinh, main = "PCF inhomogéneo")
plot(kinh, main = "K inhomogéneo")

Función envelope()

Al realizar para el caso inhomogéneo

Al realizar para el caso homogéneo

Con la ventana entregada por la función ripras

env <- envelope(ppp.locations2, fun = Kinhom, sigma = bw.scott)
## Generating 99 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40,
## 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80,
## 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98,  99.
## 
## Done.
plot(env, main="Envelope con k inhomogénea \n para las localizaciones (ripras)")

env1 <- envelope(ppp.locations2,fun = pcfinhom,nsim=39)
## Generating 39 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,  39.
## 
## Done.
env2 <- envelope(ppp.locations2,fun = Kinhom,nsim=39)
## Generating 39 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,  39.
## 
## Done.
plot(env1, main="Pcfinhom (owin ripras)")

plot(env2, main="Kinhom (owin ripras)")

Conclusiones

  • Evidentemente, la ventana de observación elaborada con el shape del área metropolitana, implicaba un costo computacional gigante, por lo que se opta por aplicar los procesos para calcular los envelopes con el objeto owin entregado por la función ripras, el cual es mucho más sobrio y permite el computo tranquilamente.

  • Los resultados obtenidos fueron bastante coherentes, pues se espera que para distancias cortas de localizaciones donde se presenten robos por cosquilleo, muy probablemente se presenten muchos más.