Análisis Espacial de los Homicidios en Cali

Fuente de datos

Para este análisis se utilizó una base datos de homicidios de la ciudad de Cali, la cual contiene la georeferenciación de los casos ocurridos durante el año 2017 en el area metropolitana de la Ciudad. Cali es una ciudad con más de 2 millones de habitantes de habitantes, cuya tasa para este año correspondió con 70 casos por cada 100mil habitantes. Los resultados que se muestran a continuación corresponden a los casos ocurridos en el mes de diciembre.

#cargando librerias

require(rgdal)
require(maptools)
require(sp)
require(raster)
require(stpp)
require(KernSmooth)
library(sf) # Development version in the ppp branch
library(spatstat)

#leyendo shape
comunas=st_read(dsn="C:/Users/delia/Desktop/Doctorado/Estadistica_Espacial/Ejercicio_NoParamétrico/Comunas/Comunas/bcs_lim_comunas_WGS84.shp")
## Reading layer `bcs_lim_comunas_WGS84' from data source `C:\Users\delia\Desktop\Doctorado\Estadistica_Espacial\Ejercicio_NoParamétrico\Comunas\Comunas\bcs_lim_comunas_WGS84.shp' using driver `ESRI Shapefile'
## Simple feature collection with 22 features and 5 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -76.59076 ymin: 3.331819 xmax: -76.46125 ymax: 3.505871
## geographic CRS: WGS 84
homicidios=st_read(dsn="C:/Users/delia/Desktop/Doctorado/Estadistica_Espacial/Ejercicio_NoParamétrico/Homicidios2017/Homicidios2017/2017_0_WGS84.shp")
## Reading layer `2017_0_WGS84' from data source `C:\Users\delia\Desktop\Doctorado\Estadistica_Espacial\Ejercicio_NoParamétrico\Homicidios2017\Homicidios2017\2017_0_WGS84.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1242 features and 47 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -1.797693e+308 ymin: -1.797693e+308 xmax: -76.46154 ymax: 3.50128
## geographic CRS: WGS 84
homicidios=st_transform(homicidios,  crs = 6345)
comunas=st_transform(comunas,  crs = 6345)

#creando la ventana del patron puntual espacial
comunas_win=as.owin(comunas)
#Gráfico espacial
par(mfrow=c(1,1))
plot(comunas_win)
plot(homicidios[,1],add=T)

#Creando un objeto de clase ppp (patron puntual)
fin=homicidios$mes=="DEC"
homicidios2=st_coordinates(homicidios[fin,])

homicidiospp <- ppp(x=homicidios2[,1],y=homicidios2[,2],window=comunas_win)

homicidiospp2=ripras(homicidiospp)
plot(homicidiospp2)

homicidiospp <- ppp(x=homicidios2[,1],y=homicidios2[,2],window=homicidiospp2)
#
plot(homicidiospp)

#eliminando duplicados para tener un proceso puntual simple
homicidiospp <- unique(homicidiospp)
plot(homicidiospp)

En el mapa de puntos se observa presencia de casos de homicidios en casi toda la ciudad, sin embargo, una mayor concentración en zonas como el oriente.

Se aplica la prueba Chi-cuadrado para evalaur si la intensidad tiene un comportamiento homogeneo o no, bajo la hipótesis nula de aleatoriedad. Esta prueba se calcula bajo el siguiente estadístico:

\(\chi^{2}=\sum\frac{(observado-Esperado)^{2}}{Esperado}\)

#prueba chi-cuadrado
qc.loc <- quadratcount(homicidiospp, nx=4, ny=5)
plot(homicidiospp, pch=3, cex=0.6)
plot(qc.loc, add=T, textargs = list(col='red'))

quadrat.test(qc.loc)
## 
##  Chi-squared test of CSR using quadrat counts
## 
## data:  
## X2 = 75.889, df = 18, p-value = 8.888e-09
## alternative hypothesis: two.sided
## 
## Quadrats: 19 tiles (irregular windows)

Debido a que el valor p es menor a 0.05, indica que la intensidad no es constante.

Función de Intensidad

Para la estimación de la intensidad se utilizó el método no parámetrico kernel, para los diferentes modelos: “gaussiano” “epanechnikov” y “quartic”. En todos se incluyó la corrección de borde de Diggle. También se variaron los anchos de banda para el modelo “gaussiano” teniendo en cuenta la opción por defecto que usa el software R a través de la función density.ppp, también se incluyeron los métodos de diggle, Scott y ppl para determinar el ancho de banda.

El kernel corregido de Diggle se estima a través de la siguiente fórmula (4):

\(\widetilde{\lambda}^{(D)}(u)=\sum_i^n\frac{1}{e(x_{i})}\kappa(u-x_{i})\)

##intensidad proceso aleatorio##
#intensidad espacial
intensidad<-density.ppp(homicidiospp,diggle=TRUE)
intensidad_homicidios <- density.ppp(homicidiospp,diggle=TRUE,sigma=bw.diggle(homicidiospp))
intensidad_homicidios1 <- density.ppp(homicidiospp,diggle=TRUE,sigma=bw.ppl(homicidiospp))
intensidad_homicidios2 <- density.ppp(homicidiospp,diggle=TRUE,sigma=bw.scott(homicidiospp))
intensidad_homicidios3 <- density.ppp(homicidiospp,diggle=TRUE,sigma=bw.ppl(homicidiospp), kernel="epanechnikov")
intensidad_homicidios4 <- density.ppp(homicidiospp,diggle=TRUE,sigma=bw.ppl(homicidiospp), kernel="quartic")

Al realizar los mapas con las diferentes funciones de intensidad estimadas se observa que la variación es principalmente explicada por el ancho de banda seleccionado que por el modelo, esto coincide con lo explicado por Cronie (5), dado que entre el gaussiano, epanechnikov y quartic no se obervan grandes diferencias, mientras que cuando el sigma es mayor se observa una mayor suavización.

#Gráfico espacial de la función de intensidad usando la opción por defecto para establecer el ancho de banda
par(mfrow=c(1,2))
plot(intensidad,xlab="x",ylab="y",main="Gaussiano sigma=0.016")
contour(density(homicidiospp),add=TRUE)

#Gráfico espacial de la función de intensidad usando usando el método diggle para el ancho de banda (método Diggle)
plot(intensidad_homicidios,xlab="x",ylab="y",main="Gaussiano sigma=0.003")
contour(density(homicidiospp),add=TRUE)

par(mfrow=c(1,2))
#Gráfico espacial de la función de intensidad usando el método ppl
plot(intensidad_homicidios1,xlab="x",ylab="y",main="Gaussiano sigma=0.007")
contour(density(homicidiospp),add=TRUE)

#Gráfico espacial de la función de intensidad usando el método scott
plot(intensidad_homicidios2,xlab="x",ylab="y",main="Gaussiano x=0.013 y=0.011")
contour(density(homicidiospp),add=TRUE)

par(mfrow=c(1,2))
#Gráfico espacial de la función de intensidad usando el método scott y modelo epanechnikov
plot(intensidad_homicidios3,xlab="x",ylab="y",main="Epanechnikov sigma=0.007")
contour(density(homicidiospp),add=TRUE)

#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda y modelo quartic
plot(intensidad_homicidios4,xlab="x",ylab="y",main="Quartic sigma=0.007")
contour(density(homicidiospp),add=TRUE)

Se genera el mapa de la función de intensidad en donde se comprueba el patrón agregado de violencia en la zona oriente de la ciudad de Cali.

##gráfico con plotly (dinamico)
intensidad_raster=rasterFromXYZ(data.frame(intensidad_homicidios1)[,1:3])

par(mfrow=c(1,1))

##Adicionar el shape
require(rasterVis)
## Loading required package: rasterVis
## Loading required package: lattice
## 
## Attaching package: 'lattice'
## The following object is masked from 'package:spatstat':
## 
##     panel.histogram
## Loading required package: latticeExtra
levelplot(intensidad_raster,col.regions = colorRampPalette(c("gray","yellow","orange","red"))(1e3))

También se realizó la estimación de la intensidad a través del método kernel adaptativo, el cual plantea que el ancho de banda “bandwidth” es variable y se estima a través de la siguiente fórmula (6):

\[ \widehat{\lambda}_{h0} (x)=\sum_{y\in{Y}} h(y)^{-2} K(h(y)^{-1} (x-y)) =\sum_{y\in{Y}} K_{h(y)}(x-y) \]

en donde el ancho de banda variable se muestra en la siguiente expresión:

\[h(y)=h_0 n^{1/2}\lambda(y)^{-1/2}\]

En la siguiente gráfica se observan las diferencias entre el kernel clásico y adaptativo, en donde el patrón de intensidad es casi el mismo, sin embargo, el método adaptativo se focaliza más en los sitios donde realmente hay una agregación de los casos de homicidios, mientras el clásico sobreestima la intensidad.

intensidad_homicidios5=densityAdaptiveKernel(homicidiospp)
par(mfrow=c(1,3))
plot(intensidad,main="Kernel Clasico")
plot(intensidad_homicidios5,main="Kernel Adaptativo")
plot(intensidad-intensidad_homicidios5,main="Diferencia")

Métodos por máxima verosimilitud para estimar la Intensidad

Para un proceso poisson homogéneo con intensidad \(\lambda\) el logaritmo de la verosimilitud es:

\(log L(\lambda;x)=n(x)log(\lambda)-\lambda area(W)\)

Para un proceso poisson inhomogéneo con función de intensidad \(\lambda_\theta(x_i)x\) dependiendo del parámetro \(\theta\). El logaritmo de la verosimilitud para \(\theta\) es:

\(log L(\lambda;x)=\sum_{i=1}^nlog\lambda_{\theta}(x_i)-\int_{W}\lambda_{\theta}(\mu)d\mu\)

donde \(\int_{W}\lambda_{\theta}(\mu)d\mu\) es el número esperado de casos del proceso poisson inhomogénero con intensidad \(\lambda_{\theta}(\mu)\) en la región W.

Se ajustaron dos modelos paramétricos para estimar la intensidad, un modelo log-lineal (1) y un log-cuadrático (2). Denotados por las siguientes expresiones:

(1)\(\lambda_{\theta}(x,y)=exp(\theta_{0}+\theta_{1}x+\theta_{2}y)\)

(2)\(\lambda_{\theta}(x,y)=exp(\theta_{0}+\theta_{1}x+\theta_{2}y+\theta_{3}x^{2}+\theta_{4}y^{2})\)

#modelo log-lineal
loglin.loc <- ppm(homicidiospp, ~x+y)

#modelo log-cuadrático
logquad.loc <- ppm(homicidiospp, ~polynom(x, y, 2))

#predicción del modelo
par(mfrow=c(1,2),mar=c(0,0,1,2))

plot(predict(loglin.loc), main='log-linear')
points(homicidiospp, pch=16, cex=0.2)

plot(predict(logquad.loc), main='log-quadratic')
points(homicidiospp, pch=16, cex=0.2)

logquad=predict(logquad.loc)

Posterior a la estimación de los modelos se revisó la bondad de ajuste de los mismos a través del análisis de residuales suavizados usando un método kernel.

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 =  [-5.719e-07, 7.082e-07]
points(homicidiospp, pch=16, cex=0.2)
diagnose.ppm(logquad.loc, which = "smooth")
## Model diagnostics (raw residuals)
## Diagnostics available:
##  smoothed residual field
## range of smoothed field =  [-4.852e-07, 5.63e-07]
points(homicidiospp, pch=16, cex=0.2)

Las funciones de intensidad selecionadas fueron el kernel gaussiano usando el ancho de banda por máxima verosimilitud y el modleo logcuadrático. Los cuales se muestran a continuación, al igual que la diferencias en las funciones de intensidad estimadas:

par(mfrow=c(1,3))
plot(intensidad_homicidios1,main="Kernel Gaussiano")
plot(logquad,main="Modelo logquad")
plot(intensidad_homicidios1-logquad,main="Diferencia")

Estadísticos de Resumen

Para determinar si el patrón corresponde a un proceso aleatorio o agregado se realizó la estimación de los estadísticos de resumen de segundo orden como la función K, la correlación por pares y la función J. Para cada una de estas se utilizó la función de intensidad estimada a través del kernel gaussiano usando el método de máxima verosimilitud (ppl) para establecer el ancho de banda y el modelo logcuadrático.

A continuación, se muestran la funciones y los respectivos resultados:

Función K

Se estimó la función K, la cual se calcula a través de la siguiente formula en el caso homogéneo (1):

\(\widehat{K}(r)=\frac{\mid W\mid}{n(n-1)}\sum_i^n\sum_j^n1 \left\{d_{ij}\leq r\right\}e_{ij}(r)\)

Para el caso inhomogéneo su cálculo se realiza a través de la siguiente función (2):

\(\widehat{K}_{inhom}(t)=\frac{1}{\mid W\mid}\sum_i^n\sum_j^n\frac{w_{y_{i}y_{j}}1 \left\{\parallel y_{i}-y_{j}\parallel\leq t\right\}}{\lambda(y_{i})\lambda(y_{j})}\)

Donde \(w_{y_{i}y_{j}}\) es el factor de corrección de borde de Ripley´s.

Correctores de Borde

A través de la libreria spatstat se tienen diferentes opciones para realizar la corrección de borde en la estimación de los estadísticos, para este caso se usó la corrección por traslación la cual tiene la siguiente función:

\(\widehat{K}_{trans}(r)=\frac{1}{\overline{\lambda}n_{i}}\sum_{i=1}^n\sum_{j=1}^n 1\left\{d_{ij}\leq r\right\}\frac{|W|}{|W\cap(W-(x_{j}-x_{i}))|}\)

Se generó la gráfica para el caso Homogéneo en donde se observa un comportamiento agregado ya que la función observada se encuentra por encima del modelo poisson, sin embargo, en el caso inhomogéneo es al contrario, indicando aleatoriedad. Se utilizó el corrector de borde de traslación.

par(mfrow=c(1,3)) 
plot(Kest(homicidiospp, correction=c("translate")), main = "K-estimada (Homogénea)")
plot(Kinhom(homicidiospp, lambda = intensidad_homicidios1, correction=c("translate")), main = "K-estimada (Inhomogénea) kernel Gaussiano")
plot(Kinhom(homicidiospp, lambda = logquad, correction=c("translate")), main = "K-estimada (Inhomogénea) logquad")

Función de Correlación por pares

Algo similar ocurre con la función de correlación por pares, la cual se calcula a través de la siguiente fórmula (1):

\(g(\xi,\eta)=\frac{\rho^{(2)}(\xi,\eta)}{\rho(\xi) \rho(\eta)}\)

Donde \(g(\xi,\eta)>1\) indica agregación o aglomeración a distancias r, si \(g(\xi,\eta)<1\) indica regularidad (1).

Al igual que en la función K homogénea, se observa que el proceso es agregado dado que se encuentra por encima del cero, en el caso inhomogéneo se observan fluctuaciones, mostrando agreación para ciertas radios y aleatoriedad en otros.

par(mfrow=c(1,3)) 
plot(pcf(homicidiospp, correction=c("translate")), main = "PCF (Homogénea)")
plot(pcfinhom(homicidiospp, lambda = intensidad_homicidios1, correction=c("translate")), main = "PCF (Inhomogénea) Kernel")
plot(pcfinhom(homicidiospp, lambda = logquad, correction=c("translate")), main = "PCF (Inhomogénea) logquad")

Función J

También se estimó la función J, esta función en el caso homogéneo se calcula de la siguiente forma (1):

\(J(r)=(1-G(r))/(1-F(r))\)

Donde \(G(r)\) es la función de la distancia del vecino más cercano y \(F(r)\) es la función de espacio vacío. Para el caso inhomogéneo, su cálculo se realiza de la siguiente forma (3):

\(J(r)_{inhom}=(1-G(r)_{inhom})/(1-F(r)_{inhom})\)

Para la función J, se observa que la función homogénea hasta cierto radio está por debajo de 1, asi mismo, para el caso inhomogéneo.

par(mfrow=c(1,3)) 
plot(Jest(homicidiospp, correction=c("border")), main = "J (Homogénea)")
plot(Jinhom(homicidiospp, lambda = intensidad_homicidios1, correction=c("border")), main = "J (Inhomogénea) Kernel")
plot(Jinhom(homicidiospp, lambda = logquad, correction=c("border")), main = "J (Inhomogénea) logquad")

Se realizó el análisis utilizando una simulación de los envelops bajo la hipótesis de aleatoriedad, a través de las funciones homogeneas e inhomogéneas.

# Pointwise envelopes 

khomicidios <- envelope(homicidiospp, Kinhom, nsim = 99, fix.n = TRUE) 
## Generating 99 simulations of CSR with fixed number of points  ...
## 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.
khomicidiosh <- envelope(homicidiospp, Kest, nsim = 99, fix.n = TRUE) 
## Generating 99 simulations of CSR with fixed number of points  ...
## 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.
ghomicidios <- envelope(homicidiospp, pcfinhom, nsim = 99, fix.n = TRUE) 
## Generating 99 simulations of CSR with fixed number of points  ...
## 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.
ghomicidiosh <- envelope(homicidiospp, pcf, nsim = 99, fix.n = TRUE) 
## Generating 99 simulations of CSR with fixed number of points  ...
## 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.
Jhomicidios <- envelope(homicidiospp, Jinhom, nsim = 99, fix.n = TRUE) 
## Generating 99 simulations of CSR with fixed number of points  ...
## 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.
Jhomicidiosh <- envelope(homicidiospp, Jest, nsim = 99, fix.n = TRUE)
## Generating 99 simulations of CSR with fixed number of points  ...
## 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.
par(mfrow=c(1,2)) 
plot(khomicidiosh, main = "Envelopes- K homogénea")
plot(khomicidios, main = "Envelopes- K inhomogenea") 

par(mfrow=c(1,2)) 
plot(ghomicidiosh, main = "Envelopes- g homogénea")
plot(ghomicidios, main = "Envelopes- g inhomogenea")

par(mfrow=c(1,2)) 
plot(Jhomicidiosh, main = "Envelopes- J homogénea")
plot(Jhomicidios, main = "Envelopes- J inhomogenea")

Finalmente, se concluye que para la estimación de la intensidad el modelo que rpesenta mejor comportamiento es el gaussiano, teniendo en cuenta el método de diggle para establecer el ancho de banda (bw.diggle).

Referencias

  1. Moller, J., & Waagepetersen, R. P. (2003). Statistical inference and simulation for spatial point processes. CRC Press.
  2. Baddeley, A. J., Moller, J., & Waagepetersen, R. (2000). Non- and semi-parametric estimation of interaction in inhomogeneous point patterns. Statistica Neerlandica, 54(3), 329–350. doi:10.1111/1467-9574.00144.
  3. Van Lieshout, M. N. M. (2011). A J-function for inhomogeneous point processes. Statistica Neerlandica, 65(2), 183–201. doi:10.1111/j.1467-9574.2011.00482.
  4. Baddeley, A., Rubak, E., & Turner, R. (2015). Spatial point patterns: methodology and applications with R. CRC press. 5.Cronie, O., & Van Lieshout, M. N. M. (2018). A non-model-based approach to bandwidth selection for kernel estimators of spatial intensity functions. Biometrika, 105(2), 455-462.
  5. Davies, T. M., & Baddeley, A. (2018). Fast computation of spatially adaptive kernel estimates. Statistics and Computing, 28(4), 937-956.