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 queden sólo 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)

plot(homicidiospp)

Se generó la gráfica de la función K para el caso Homogéneo e Inhomogéneo en donde se observa un comportamiento agregado ya que se encuentra por encima del poisson.

plot(Kest(homicidiospp), main = "K-estimada (Homogénea)")

plot(Kinhom(homicidiospp), main = "K-estimada (Inhomogénea)")
## number of data points exceeds 1000 - computing border correction estimate only

Algo similar ocurre con la función g en la cual se concluye que el proceso es agregado porque se encuentra por encima del cero.

plot(pcf(homicidiospp), main = "pair-correlation estimada (Homogénea)")

plot(pcfinhom(homicidiospp), main = "pair-correlation 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,2)) 
#20 simulaciones
khomicidios <- envelope(homicidiospp, Kest, 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, pcf, 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 usando K") #agregado
plot(ghomicidios, main = "Envelopes usando g")

Se estima la intensidad por medio del método de diggle.

##intensidad proceso aleatorio##
#intensidad espacial
intensidad_homicidios <- density(homicidiospp,diggle=TRUE)
#Gráfico espacial de la función de intensidad
par(mfrow=c(1,1)) 
plot(intensidad_homicidios,xlab="x",ylab="y",main="Intensidad espacial")

summary(intensidad_homicidios)
## real-valued pixel image
## 128 x 128 pixel array (ny, nx)
## enclosing rectangle: [-76.59076, -76.46125] x [3.331819, 3.505871] units
## dimensions of each pixel: 0.00101 x 0.001359781 units
## Image is defined on the full rectangular grid
## Frame area = 0.0225417997208965 square units
## Pixel values
##  range = [1.638757, 245287.8]
##  integral = 1162
##  mean = 51548.68
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_homicidios)[,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'))

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.

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