Trabajo Final Geoestadística

Trabajo Final Geostadística

Introducción:

Este análisis forma parte de la evaluación final del módulo Geoestadística. La idea es llevar adelante un análisis de los siniestros sufridos por ciclistas en la Ciudad de Buenos Aires utilizando la técnica de Procesos y Patrones Puntuales. Este trabajo se encuentra disponible en mi blog personal

Objetivo:

Visualizar espacialmente los siniestros sufridos por ciclistas por parte de transporte público/taxi como contraparte

Analizar y modelar bajo la lógica de patrones espaciales de puntos estos fenómenos

Indagar acerca de la existencia de agrupamientos de siniestros.

Metodología:

Para llevar adelante los objetivos planteados, se utilizará el conjunto de técnicas específicas de los análisis de patrones puntuales. Estas técnicas permiten conocer e identificar altas zonas de densidad de un fenómeno particular así como también la distancia entre sucesos y su posible correlación con otras variables. Para ello, resultará de gran utilidad los resultados y el análisis de las estadísticas tanto de primer como de segundo orden. Siguiendo este orden de ideas, y dado el incremento de casos de siniestros sufridos por ciclistas, resulta de interés el abordaje de esta problemática a partir de la recolección de datos abiertos del portla de CABA. De esta manera, se arribará a un resumen espacial completo considerando la concentración y densidad de los siniestros acontecidos entre ciclistas y conductores de transporte público.

Fuentes de datos:

A los fines de los ideas planteadas, trabajaremos con los datos de siniestros ocurridos durante el año 2023, tomando como víctimas aquellos individuos que se transportan en bicicleta y, como contra-parte filtramos transporte público y taxi. Tambien utilizaremos las geometrías de comunas y barrios de la ciudad para visualizar nuestra información. Los datos serán proyectados en el CRS 5347.

Librerías y datos:

Procedemos a limpiar el ambiente de trabajo e instalar las librerías necesarias:

rm(list = ls())

# Instalamos librerías de trabajo
if (!require("pacman")) install.packages("pacman")
pacman::p_load("maptools",
               "lubridate",
               "tidyverse",
               "leaflet",
               "leaflet.extras",
               "sf",
               "units",
               "nngeo",
               "spatstat")

Limpieza y Procesamiento de datos:

Luego de activar las librerías de trabajo, procedemos a leer nuestros datos considerando las filtraciones planteadas en los osbjetivos para lo cual filtramos los datos que figuren “BICICLETA” como víctima y “TRANSPORTE PÚBLICO” y “TAXI” como contraparte. Tambien eliminamos aquellos registros que no tenga datos de latitud y longitud, tomamos el año 2023 y , finalmente convertimos nuestros datos en objeto espacial utilizando el sistema de coordenadas 5347.

Junto con los archivos espaciales auxiliares de comunas y barrios, ya podemos dar inicio a nuestro análisis.

siniestros <- read.csv("data/siniestros_viales_hechos.csv", sep = ";") %>% 
  filter(latitud != "SD" & latitud != "" & latitud != "sd",
         victima == "BICICLETA",
         contraparte == "TRANSPORTE PUBLICO" | contraparte == "TAXI") %>%
  mutate(longitud = as.numeric(longitud),
          latitud = as.numeric(latitud),
          year = year(fecha)) %>%
  filter( year == 2023) %>% 
  st_as_sf(coords =c("longitud", "latitud"), crs = 4326) %>% 
  st_transform(5347)

# Desde disco: data/comunas.csv
comunas <- st_read("data/comunas.geoson", 
                  crs = 4326) %>% 
  st_transform(5347)

## Reading layer `comunas' from data source 
##   `C:\Users\Guille\Documents\Proyectos\FLACSO\geo\TP FINAL\data\comunas.geoson' 
##   using driver `GeoJSON'
## Simple feature collection with 15 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53152 ymin: -34.70529 xmax: -58.33516 ymax: -34.52649
## Geodetic CRS:  WGS 84

# Desde disco: data/barrios.csv
barrios <- st_read("data/barrios.geoson", 
                  crs = 4326) %>%
  st_transform(5347)

## Reading layer `barrios' from data source 
##   `C:\Users\Guille\Documents\Proyectos\FLACSO\geo\TP FINAL\data\barrios.geoson' 
##   using driver `GeoJSON'
## Simple feature collection with 48 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -58.53152 ymin: -34.70529 xmax: -58.33516 ymax: -34.52649
## Geodetic CRS:  WGS 84

Visualizamos de forma general los datos en formato ggplot:

ggplot() + 
  geom_sf(data = barrios, color = "#afafaf", fill = NA)+
  geom_sf(data = comunas, color = "#585858", fill = NA) +
  geom_sf(data = siniestros, alpha = 0.7, color = "#660000") + 
  labs(title = "Distribución de siniestros sufridos por bicicletas",
       shape = "Tipo de comisaría") +
  theme_void()

Ya tenemos graficados los 94 siniestros registrados según los filtros aplicados. A continuación prepararemos los datos para el análisis de patrones de puntos En primer lugar, resulta necesario definir la ventana en la cual vamos a trabajar. Tal cual se estudió en el curso, la ventana elegida como región de estudio afecta la estimación de la intensidad.

CABA <- summarise(comunas)
CABA <- nngeo::st_remove_holes(CABA)

ventana <- as.owin(st_geometry(CABA))
unitname(ventana) <- "Meter"

plot(ventana)

Generamos el patrón de puntos (clase ppp) y generamos ruido para mover de forma aleatoria los puntos duplicados:

siniestros_ppp <- as.ppp(st_geometry(siniestros), W = ventana)

set.seed(200)
siniestros_ppp <- rjitter(siniestros_ppp, retry=TRUE, nsim = 1, 
                       radius = 2, drop = TRUE)

cat("\nExisten puntos duplicados:", any(duplicated.ppp(siniestros_ppp)))

## 
## Existen puntos duplicados: FALSE

plot(siniestros_ppp, cols = "#660000")

Ya tenemos nuestro objeto ppp en condiciones para comenzar a efectuar las correspondientes estadísticas de resumen. Por lo pronto, puede verse algún patrón en la zona céntrica de la ciudad.

Estadísticas de resumen de primer orden

Las estadísticas de primer orden permiten calcular intensidades proporcionando el numero esperado de puntos por unidad espacial sin considerar la estructura de dependencia de dichos puntos. Esto es útil para medir y observar las formas que adquieren las densidades del fenómeno en cuestión.

Intensidad del patrón de puntos

Veamos cuál es la intensidad del patrón de siniestros por unidad de área:

cat("Intensidad de delitos por metro cuadrado:", 
    intensity(siniestros_ppp) |> format(scientific = F), 
    "\n")

## Intensidad de delitos por metro cuadrado: 0.0000004601921

cat("Intensidad de delitos por kilómetro cuadrado:", 
    format( (intensity(siniestros_ppp) * (1000^2) ), 
            scientific = F, digits = 4, decimal.mark = ",") )

## Intensidad de delitos por kilómetro cuadrado: 0,4602

Intensidad por cuadrantes

Lo que haremos a continuación será dividir nuestra ventana en cuadrantes midiendo la intensidad en cada subdivisión:

QC_delitos <- quadratcount(siniestros_ppp, nx = 5, ny = 5)

QC_delitos

## tile
## Tile row 1, col 1 Tile row 1, col 2 Tile row 1, col 3 Tile row 1, col 4 
##                 0                 5                 1                 0 
## Tile row 2, col 1 Tile row 2, col 2 Tile row 2, col 3 Tile row 2, col 4 
##                 3                 7                10                 5 
## Tile row 2, col 5 Tile row 3, col 1 Tile row 3, col 2 Tile row 3, col 3 
##                 1                 2                12                16 
## Tile row 3, col 4 Tile row 3, col 5 Tile row 4, col 1 Tile row 4, col 2 
##                18                 2                 3                 3 
## Tile row 4, col 3 Tile row 4, col 4 Tile row 4, col 5 Tile row 5, col 1 
##                 0                 3                 2                 0 
## Tile row 5, col 2 Tile row 5, col 3 
##                 1                 0

plot(siniestros_ppp, 
     main = "Conteo por cuadrante",
     cols = "tomato4")
plot(QC_delitos, add = TRUE, cex = 1)

Tambien podemos ver esto a partir de hexágonos considerando 1300 metros como tamaño.

# Generamos hexágonos con lado 1300m
H <- hextess(ventana, 1300)

QC_delitos_HEX <- quadratcount(siniestros_ppp, tess = H)

plot(siniestros_ppp, 
     main = "Conteo por cuadrante (hexagonal)",
     cols = "tomato4")
plot(QC_delitos_HEX, add = TRUE, cex = 1)

intensity(QC_delitos_HEX, image = T) |> 
  plot(main = "Intensidad por cuadrante (hexagonal)")

plot(siniestros_ppp, add = T, cols = "white")

Los tres gráficos evidencian una cierta concentración de casos que aumenta notáblemente en el centro. Veremos si esto puede considerarse un patrón homogéneo de Poisson (CRS) tomando el test de hipótesis siguiente:

• Hipótesis nula (\(H_0\)): la intensidad es homogénea y la distribución de casos en los cuadrantes responde a una distribución de Poisson (CSR).
• Hipótesis alternativa (\(H_1\)) : la intensidad no es homogénea (en una forma no especificada).

El nivel de significación planteado será de \(0.05\)

quadrat.test(QC_delitos_HEX, method = "MonteCarlo", nsim = 1000)

## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  
## X2 = 117.72, p-value = 0.06194
## alternative hypothesis: two.sided
## 
## Quadrats: 68 tiles (irregular windows)

Es posible graficar los resultados del test (usamos pocos cuadrantes para facilitar la lectura), obteniendo el conteo real, el esperado y los residuos \(\chi^2\).

quadrat.test(siniestros_ppp, nx = 3, ny = 3, method = "MonteCarlo") |> 
  plot(cex = 1, main = "Test de hipótesis (valores por cuadrante)")

Dado que el p_value obtenido es mayor que el nivel de significancia, no resulta posible rechazar la hipótesis nula con lo cual hay evidencia que indica que el patrón de distribución espacial puede asimirse como homogeneo de Poisson.

Modelado de densidad según función Kernel

Mas allá de la homogeneidad del patrón, es posible la estimación de la función subyacente de la intensidad (densidad).

DD <- density(siniestros_ppp, edge = F)
plot(DD, main='densidad de siniestros Bicicleta - Transporte Publico o Taxi')

class(DD)

## [1] "im"

attr(DD, "sigma")

## [1] 2248.88

Al utilizar el valor True (T) en el parámetro edge se genera una corrección de borde que ajusta las estimaciones de densidades y evita posibles sesgos. Vamos a graficar los valores con diferentes métodos

par(mfrow = c(2, 3), mar = c(0.2, 0.2, 0.2, 0.2))

  density(siniestros_ppp, edge = T) |> plot()
  density(siniestros_ppp, diggle = T) |> plot()
  
  density(siniestros_ppp, diggle = T, sigma = bw.diggle) |> plot()
  density(siniestros_ppp, diggle = T, sigma = bw.ppl) |> plot()

  adaptive.density(siniestros_ppp, f = 0.1, method = "voronoi") |> plot()
  adaptive.density(siniestros_ppp, method = "kernel") |> plot()

par(mfrow = c(1, 1))

Vermos a continuación cómo varían los valores con los diferentes métodos:

cat("Cross Validated (Diggle):", bw.diggle(siniestros_ppp), "\n")

## Cross Validated (Diggle): 717.353

cat("Likelihood Cross Validation:", bw.ppl(siniestros_ppp), "\n")

## Likelihood Cross Validation: 1608.742

cat("Scott's Rule:", bw.scott(siniestros_ppp), "\n")

## Scott's Rule: 1743.177 1421.322

Estadísticas de resumen de segundo orden

Las funciones categorizadas como de segundo orden determinan la estructura de dependencia espacial inherente al patrón puntual. En otras palabras, estas estadísticas ponen el foco ya no tanto en cada punto individual sino en la posible correlación entre puntos.

Función K

Dentro de las ya definidas estadísticas de segundo orden se encuentra la función k que, básicamente, define la cantidad de eventos (en este caso siniestros sufridos por ciclistas) en un radio r alrededor de cualquier otro evento.

K_siniestros <- Kest(siniestros_ppp)

class(K_siniestros)

## [1] "fv"         "data.frame"

plot(K_siniestros)

El siguiente paso consiste en contrastar con la banda de confianza para la curva teórica con la función envelope que contendrá el patrón de puntos para determinar el intervalo y la función Kest ya graficada. Se generarán 50 simulaciones.

AcceptPoisson_K <- envelope(siniestros_ppp, Kest, nsim = 50)

## Generating 50 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.
## 
## Done.

plot(AcceptPoisson_K)

Una alternativa es observar el intervalo de confianza que presenta la distribución de la función en nuestro patrón de puntos (mediante bootstrap de Loh). Para eso usamos la función lohboot, pasando como parámetros el objeto de clase ppp y la función que queremos estimar.

IC_K <-  lohboot(siniestros_ppp, Kest, nsim = 20)

## 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.

plot(IC_K)

Los gráficos señalan que los sucesos tienden a formar agregaciones o clusters dado que por arriba de 1.000 la curva con los siniestros supera la curva teórica. Esto implica que pueden existir zonas puntuales peligrosas con mayor cantidad de siniestros.

Veamos los resultados de las siguientes funciones.

Función L de Besag

Los valores estimados superiores a los valores teóricos señalan la existencia de agregación o agrupamiento de siniestros. La existencia de agregación de puntos en ciertos lugares indica que la intensidad no es homogénea.

AcceptPoisson_L <- envelope(siniestros_ppp, Lest, nsim = 20)

## Generating 20 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
## 20.
## 
## Done.

plot(AcceptPoisson_L) 

Nuestra distribución vuelve a superar la distribución teórica.

Función g (correlación de a pares)

Si los puntos son independientes unos de otros (como es el caso del proceso Poisson), entonces \(g(r) = 1\). Por el contrario, \(g(u,v) > 1\) indica la presencia de atracción entre puntos (agregaciones) hasta una distancia \(r\), mientras que \(g (r) < 1\) indica repulsión entre puntos (regularidad).

AcceptPoisson_g <- envelope(siniestros_ppp, pcf, nsim = 20)

## Generating 20 simulations of CSR  ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
## 20.
## 
## Done.

plot(AcceptPoisson_g)

¿Que patrón sugiere la función g?

Función G

Mide la distribución de distancias desde un evento arbitrario hasta su evento más cercano.

Si la función empírica es mayor que la teórica, las distancias al vecino más cercano en el patrón son más cortas en comparación con un proceso de Poisson homogéneo. Si la función empírica es menor que la teórica, las distancias al vecino más cercano en el patrón son más largas en comparación con un proceso de Poisson homogéneo.

AcceptPoisson_G <- envelope(siniestros_ppp, Gest, nsim = 100)

## Generating 100 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, 
## 100.
## 
## Done.

plot(AcceptPoisson_G) 

Función F

Se llama también función de espacio vacío o distancia al próximo evento (siniestro) desde un punto fijo en el espacio. Resume información sobre huecos o vacíos (gaps) entre los puntos del dataset.

AcceptPoisson_F <- envelope(siniestros_ppp, Fest, nsim = 100)

## Generating 100 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, 
## 100.
## 
## Done.

plot(AcceptPoisson_F) 

El gráfico vuelve a mostrar patrones de agregación dado que nuestra función es menor que la teórica lo cual demuestra que hay un mayor porcentaje de distancias cortas al vecino mas cercano. En otras palabras, hay siniestros agrupados en algunas zonas lo cual evidencia zonas de peligro para ciclistas.

Función J

La función J tiene una poderosa ventaja en el hecho de considerar la función f y g en su cálculo tal cual lo indica la siguiente fórmula:

\[ J(r) = \frac{(1-G(r))}{1-F(r)}; \quad \text{para } F(r) < 1 \]

Para un proceso Poisson homogéneo (aleatorio completo) \(J(r) = 1\). Cuando \(J(r) < 1\) indica agregación o agrupamiento, mientras que \(J (r)> 1\) indica regularidad o dispersión

J_siniestros <- Jest(siniestros_ppp, 100)
plot(J_siniestros)

AcceptPoisson_J<- envelope(siniestros_ppp, Jest, nsim = 50)

## Generating 50 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.
## 
## Done.

plot(AcceptPoisson_J) 

Todas las funciones juntas

todas_siniestros <- allstats(siniestros_ppp)
plot(todas_siniestros)

plot(todas_siniestros, subset=c("r<=300", "r<=300", "r<=300", "r<=300"))

Conclusiones:

El análisis de patrón de puntos aplicado al análisis de siniestros sufridos por ciclistas en CABA permitió un abordaje más profundo del fenómeno. Las primeras pruebas de densidad arrojaron patrones de homogeneidad, mientras que las estadística de segundo orden evidenciaron posibles formas de agregación puntuales que muestran la existencia de zonas geográficas de importante peligro para los ciclistas. Como futuras ideas para continuar el trabajo puede pensarse el estudio de las avenidas y sus velocidades máximas, como así también el recorrido de colectivos calculando estadísticamente las zonas de mayor peligro para ciclistas.

Bibliografia

Fernandez Aviles, G., Montero, J. (2024). Fundamentos de ciencia de datos con R. Mc Graw Hill, Disponible: [https://cdr-book.github.io/index.html]

Baddeley, A., Rubak, E., & Turner, R. (2016). Spatial Point Patterns. Methodology and Applications with R. Taylor & Francis Group, LLC. https://doi.org/10.1201/b19708

Medina, J. & Solymosi, R. (2022). Crime Mapping and Spatial Data Analysis using R. Disponible: https://maczokni.github.io/crime_mapping/

Pebesma, E. J., & Bivand, R. S. (2020). Spatial Data Science. https://keen-swartz-3146c4.netlify.app/