Para este análisis se usaron datos de homicidios para la ciudad de Cali y contiene la georeferenciación de los casos de ocurridos durante el año 2017 en el area metropolitana de la Ciudad.

#cargando librerias
require(spatstat)
require(rgdal)
require(maptools)
require(sp)
require(raster)
require(stpp)
require(KernSmooth)

#leyendo shape
comunas=shapefile("C:/Users/delia/Desktop/Doctorado/Estadistica_Espacial/Ejercicio_NoParamétrico/Comunas/Comunas/bcs_lim_comunas_WGS84.shp")

homicidios=shapefile("C:/Users/delia/Desktop/Doctorado/Estadistica_Espacial/Ejercicio_NoParamétrico/Homicidios2017/Homicidios2017/2017_0_WGS84.shp")

#creando la ventana del patron puntual espacial
comunas_win=owin(xrange = extent(comunas)[1:2],yrange =extent(comunas)[3:4])

Se realizó un proceso de depuración de los datos garantizando que sólo queden los puntos al interior de la ciudad por medio de las operaciones espaciales crop y mask.

homicidios=crop(homicidios,comunas)

#Gráfico espacial
par(mfrow=c(1,1))
plot(homicidios,window=comunas)
## Warning in title(...): "window" is not a graphical parameter

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.

#Creando un objeto de clase ppp (patron puntual)
homicidios2=coordinates(homicidios)
homicidiospp <- ppp(x=homicidios2[,1],y=homicidios2[,2],window=comunas_win)

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

Se aplica la prueba Chi-cuadrado para evalaur si la intensidad tiene un comportamiento homogeneo o no, bajo la hipótesis nula de aleatoriedad.

#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 = 2126.9, df = 19, p-value < 2.2e-16
## alternative hypothesis: two.sided
## 
## Quadrats: 4 by 5 grid of tiles

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

Por otro lado, se estimó la función K, la cual se calcula a través de la siguiente formula:

\(\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)\)

Se generó la gráfica para el caso Homogéneo e Inhomogéneoen donde se observa un comportamiento agregado ya que se encuentra por encima del modelo poisson en ambos casos.

par(mfrow=c(1,2)) 
plot(Kest(homicidiospp), main = "K-estimada (Homogénea)")
plot(Kinhom(homicidiospp), main = "K-estimada (Inhomogénea)")

Algo similar ocurre con la función g en la cual se concluye que el proceso es agregado dado se encuentra por encima del cero, tanto para el caso homogéneo como inhomogéneo.

par(mfrow=c(1,2)) 
plot(pcf(homicidiospp), main = "G-estimada (Homogénea)")
plot(pcfinhom(homicidiospp), main = "G-estimada (Inhomogénea)")

Asi mismo, para la función J se observa que la función homogénea e inhomogénea están por debajo de 1, indicando agregación.

par(mfrow=c(1,2)) 
plot(Jest(homicidiospp), main = "J-estimada (Homogénea)")
plot(Jinhom(homicidiospp), main = "J-estimada (Inhomogénea)")

Se realiza el análisis utilizando una simulación de los envelops y nuevamente se concluye que el patrón es agregado.

# Pointwise envelopes 
par(mfrow=c(1,3)) 
#20 simulaciones
khomicidios <- envelope(homicidiospp, Kinhom, nsim = 20, fix.n = TRUE) 
## Generating 20 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.
## 
## Done.
ghomicidios <- envelope(homicidiospp, pcfinhom, nsim = 20, fix.n = TRUE) 
## Generating 20 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.
## 
## Done.
Jhomicidios <- envelope(homicidiospp, Jinhom, nsim = 20, fix.n = TRUE) 
## Generating 20 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.
## 
## Done.
plot(khomicidios, main = "Envelopes- K inhomogenea") 
plot(ghomicidios, main = "Envelopes- g inhomogenea")
plot(Jhomicidios, main = "Envelopes- J inhomogenea")

Se estima la intensidad usando el kernel “gaussiano” usando diferentes anchos de banda, además usando el método “epanechnikov” y “quartic”. En todos se incluyó la corrección de borde de Diggle. El kernel corregido de Diggle se estima a través de la siguiente fórmula:

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

##intensidad proceso aleatorio##
#intensidad espacial
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.ppl(homicidiospp), kernel="epanechnikov")
intensidad_homicidios3 <- density.ppp(homicidiospp,diggle=TRUE,sigma=bw.ppl(homicidiospp), kernel="quartic")
#Gráfico espacial de la función de intensidad usando la validación cruzada para la selección del ancho de banda (método Diggle)
par(mfrow=c(2,2))
plot(intensidad_homicidios,xlab="x",ylab="y",main="Modelo Gaussiano sigma=0.0013")
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 (método verosimilitud)
plot(intensidad_homicidios1,xlab="x",ylab="y",main="Modelo Gaussiano sigma=0.005")
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 método epanechnikov
plot(intensidad_homicidios2,xlab="x",ylab="y",main="Método epanechnikov sigma=0.005")
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 método quartic
plot(intensidad_homicidios3,xlab="x",ylab="y",main="Método quartic sigma=0.005")
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
intensidad_raster=mask(intensidad_raster,comunas)
levelplot(intensidad_raster,col.regions = colorRampPalette(c("gray","yellow","orange","red"))(1e3))+layer(sp.polygons(comunas, lwd=0.5, col='black'))

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)

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 =  [-77120, 104900]
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 =  [-27270, 36440]
points(homicidiospp, pch=16, cex=0.2)