Este documento implementa un análisis completo de patrones puntuales espaciales aplicado a datos de siniestralidad vial en Bogotá para el año 2023, siguiendo la metodología del libro “Spatial Point Patterns: Methodology and Applications with R” de Adrian Baddeley, Ege Rubak y Rolf Turner.
En este conjunto de datos, cada siniestro se representa como un punto ubicado en el espacio mediante sus coordenadas de longitud y latitud. Además, a cada punto se le asocian marcas, que son atributos propios del evento y permiten caracterizarlo: por ejemplo, la clase de accidente (choque, atropello, volcamiento), la gravedad del siniestro, la fecha y hora de ocurrencia, la localidad y los indicadores de participación o condiciones (p. ej., presencia de bicicleta, moto, peatón, embriaguez, exceso de velocidad o huecos). Por su parte, las covariables son variables que describen el entorno donde ocurren los eventos y se usan para explicar variaciones sistemáticas en la intensidad espacial; en este estudio, Localidad puede interpretarse como una covariable espacial de tipo areal (asignada por zonas), mientras que variables como Hora_Acc o Dia_Semana_Acc pueden incorporarse como covariables temporales si se extiende el análisis a un enfoque espacio–temporal.
# Cargar paquetes necesarios
library(spatstat) # Análisis de patrones puntuales
library(sf) # Manejo de datos espaciales
library(sp) # Clases espaciales
library(tidyverse) # Manipulación de datos
library(lubridate) # Manejo de fechas
library(viridis) # Paletas de colores
library(knitr) # Tablas
library(kableExtra) # Tablas mejoradas
library(readxl)
# Configurar opciones
options(scipen = 999)# Directorio del proyecto
setwd("C:/Users/sarar/Documents/Maestría Estadística/Espacial/Patrones Puntuales")
# Cargar datos de siniestros
siniestros <- read_excel("SIGAT_ANUARIO_2023.xlsx")
# Exploración inicial
cat("Dimensiones del dataset:", dim(siniestros), "\n")Dimensiones del dataset: 14106 38 # Vista rápida
head(siniestros)# A tibble: 6 × 38
Codigo_Accidente Formulario Longitud Latitud Direccion Fecha_Acc
<dbl> <chr> <dbl> <dbl> <chr> <dttm>
1 10591851 A001571139 -74.2 4.60 AV AVENIDA C… 2023-05-26 00:00:00
2 10591854 A001571298 -74.2 4.57 KR 44 D - DG… 2023-05-31 00:00:00
3 10592149 A001570489 -74.1 4.62 CL 26 - KR 2… 2023-06-08 00:00:00
4 10592358 A001572190 -74.1 4.59 CL 6 - KR 5 … 2023-06-14 00:00:00
5 10592361 A001570128 -74.0 4.76 CL 187 A - K… 2023-06-14 00:00:00
6 10592370 A001571474 -74.1 4.74 CL 139 - KR … 2023-06-13 00:00:00
# ℹ 32 more variables: AA_Acc <dbl>, MM_Acc <chr>, DD_Mes_Acc <dbl>,
# Dia_Semana_Acc <chr>, Hora_Acc <dbl>, Min_Acc <dbl>, Localidad <chr>,
# Clase_Acc <chr>, Elemento_Choque <chr>, Tipo_Objeto_Fijo <chr>,
# Gravedad_Indicador_Tradicional <chr>, Gravedad_indicador_30d <chr>,
# Con_Bicicleta <chr>, Con_Carga <chr>, Con_Embriaguez <chr>,
# Con_Huecos <chr>, Con_Menores <chr>, Con_Moto <chr>, Con_Peaton <chr>,
# Con_Persona_Mayor <chr>, Con_Rutas <lgl>, Con_Tpi <chr>, Con_Tpp <chr>, …# Mostrar estructura
glimpse(siniestros)Rows: 14,106
Columns: 38
$ Codigo_Accidente <dbl> 10591851, 10591854, 10592149, 10592358,…
$ Formulario <chr> "A001571139", "A001571298", "A001570489…
$ Longitud <dbl> -74.16213, -74.15834, -74.07903, -74.07…
$ Latitud <dbl> 4.599679, 4.570567, 4.623901, 4.590091,…
$ Direccion <chr> "AV AVENIDA CIUDAD DE VILLAVICENCIO - C…
$ Fecha_Acc <dttm> 2023-05-26, 2023-05-31, 2023-06-08, 20…
$ AA_Acc <dbl> 2023, 2023, 2023, 2023, 2023, 2023, 202…
$ MM_Acc <chr> "Mayo", "Mayo", "Junio", "Junio", "Juni…
$ DD_Mes_Acc <dbl> 26, 31, 8, 14, 14, 13, 14, 14, 13, 23, …
$ Dia_Semana_Acc <chr> "viernes", "miércoles", "jueves", "miér…
$ Hora_Acc <dbl> 20, 9, 20, 6, 20, 10, 18, 5, 12, 8, 21,…
$ Min_Acc <dbl> 43, 40, 0, 41, 15, 35, 20, 48, 35, 26, …
$ Localidad <chr> "BOSA", "CIUDAD BOLIVAR", "TEUSAQUILLO"…
$ Clase_Acc <chr> "Choque", "Atropello", "Volcamiento", "…
$ Elemento_Choque <chr> "Vehículo", NA, NA, "Vehículo", "Vehícu…
$ Tipo_Objeto_Fijo <chr> NA, NA, NA, NA, NA, NA, NA, "Vehiculo e…
$ Gravedad_Indicador_Tradicional <chr> "Con Heridos", "Con Heridos", "Con Heri…
$ Gravedad_indicador_30d <chr> "Con Heridos", "Con Heridos", "Con Heri…
$ Con_Bicicleta <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Carga <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Embriaguez <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Huecos <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Menores <chr> NA, NA, NA, "SI", NA, NA, "SI", NA, NA,…
$ Con_Moto <chr> "SI", "SI", "SI", NA, "SI", "SI", "SI",…
$ Con_Peaton <chr> NA, "SI", NA, NA, NA, "SI", "SI", NA, N…
$ Con_Persona_Mayor <chr> NA, NA, NA, NA, NA, NA, NA, NA, "SI", N…
$ Con_Rutas <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Tpi <chr> "SI", NA, NA, "SI", "SI", NA, NA, NA, N…
$ Con_Tpp <chr> NA, NA, NA, NA, NA, NA, NA, NA, "SI", N…
$ Con_Velocidad <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Sitp <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Troncal <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Alimentador <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Zonal <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Provisional <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Articulado <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Biarticulado <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ Con_Padron_Dual <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,…# Limpiar datos: eliminar NAs en coordenadas
siniestros_clean <- siniestros %>%
filter(!is.na(Longitud) & !is.na(Latitud)) %>%
filter(Longitud != 0 & Latitud != 0) %>%
mutate(Fecha_Acc = as.Date(Fecha_Acc, format = "%Y-%m-%d"))
cat("\nDatos limpios:", nrow(siniestros_clean), "siniestros\n")
Datos limpios: 14106 siniestros# Resumen de coordenadas
summary(siniestros_clean[, c("Longitud", "Latitud")]) %>%
kable(caption = "Estadísticas de Coordenadas") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Longitud | Latitud | |
|---|---|---|
| Min. :-74.25 | Min. :4.044 | |
| 1st Qu.:-74.14 | 1st Qu.:4.596 | |
| Median :-74.11 | Median :4.628 | |
| Mean :-74.11 | Mean :4.636 | |
| 3rd Qu.:-74.08 | 3rd Qu.:4.679 | |
| Max. :-74.01 | Max. :4.823 |
Para un análisis espacial correcto, es fundamental trabajar con un sistema de coordenadas proyectado en metros en lugar de coordenadas geográficas (grados). Esto permite:
# (0) Preparar marcas (recomendado para análisis de marcas)
siniestros_clean <- siniestros_clean %>%
mutate(
Localidad = factor(Localidad),
Clase_Acc = factor(Clase_Acc),
Gravedad_Indicador_Tradicional = factor(Gravedad_Indicador_Tradicional),
Gravedad_indicador_30d = factor(Gravedad_indicador_30d),
Dia_Semana_Acc = factor(Dia_Semana_Acc,
levels = c("lunes","martes","miércoles","jueves",
"viernes","sábado","domingo"))
)
# 1) Convertir a objeto sf con coordenadas geográficas (WGS84)
pts_sf <- st_as_sf(
siniestros_clean,
coords = c("Longitud", "Latitud"),
crs = 4326, # WGS84
remove = FALSE
)
# 2) Proyectar a CRS en metros (MAGNA-SIRGAS / Colombia Bogota zone)
# EPSG:3116 es apropiado para Bogotá
pts_m <- st_transform(pts_sf, 3116)
# 3) Extraer coordenadas proyectadas en metros
xy <- st_coordinates(pts_m)
x <- xy[, 1]
y <- xy[, 2]
cat("✓ Coordenadas extraídas:\n")✓ Coordenadas extraídas:cat(" Rango X:", round(min(x)), "a", round(max(x)), "metros\n") Rango X: 981149 a 1007238 metroscat(" Rango Y:", round(min(y)), "a", round(max(y)), "metros\n") Rango Y: 938934 a 1025101 metros# 4) Crear ventana de observación (bounding box en metros)
W <- owin(xrange = range(x), yrange = range(y))
cat("✓ Ventana de observación creada\n")✓ Ventana de observación creadacat(" Área:", round(area(W) / 1e6, 2), "km²\n") Área: 2247.99 km²# 5) Crear objeto ppp (patrón puntual)
ppp_siniestros <- ppp(x = x, y = y, window = W)
cat("✓ Objeto ppp creado con", npoints(ppp_siniestros), "eventos\n")✓ Objeto ppp creado con 14106 eventoscat(" Intensidad global:", round(intensity(ppp_siniestros) * 1e6, 2), "siniestros/km²\n") Intensidad global: 6.27 siniestros/km²# 6) Añadir marcas (atributos de cada evento)
marks(ppp_siniestros) <- siniestros_clean %>%
select(
Localidad,
Clase_Acc,
Gravedad_Indicador_Tradicional,
Gravedad_indicador_30d,
Hora_Acc,
Dia_Semana_Acc
)
cat("✓ Marcas añadidas:", ncol(marks(ppp_siniestros)), "variables\n\n")✓ Marcas añadidas: 6 variables# 7) Diagnóstico: puntos duplicados (misma coordenada exacta)
n_dup <- sum(duplicated(cbind(x, y)))
cat("⚠ Puntos duplicados (misma coordenada):", n_dup, "\n")⚠ Puntos duplicados (misma coordenada): 307 cat(" Proporción duplicada:", round(100 * n_dup / npoints(ppp_siniestros), 2), "%\n") Proporción duplicada: 2.18 %if(n_dup > 0) {
cat(" Nota: Los duplicados pueden indicar múltiples siniestros en el mismo lugar.\n")
cat(" Se pueden tratar con jitter o análisis de puntos múltiples.\n")
} Nota: Los duplicados pueden indicar múltiples siniestros en el mismo lugar.
Se pueden tratar con jitter o análisis de puntos múltiples.# 8) Resumen del patrón puntual
cat("\n=== RESUMEN DEL PATRÓN PUNTUAL ===\n")
=== RESUMEN DEL PATRÓN PUNTUAL ===summary(ppp_siniestros)Marked planar point pattern: 14106 points
Average intensity 0.000006274925 points per square unit
*Pattern contains duplicated points*
Coordinates are given to 10 decimal places
Mark variables: Localidad, Clase_Acc, Gravedad_Indicador_Tradicional,
Gravedad_indicador_30d, Hora_Acc, Dia_Semana_Acc
Summary:
Localidad Clase_Acc Gravedad_Indicador_Tradicional
KENNEDY :1969 Atropello :2760 Con Heridos:12447
ENGATIVA :1295 Caida de ocupante: 798 Con Muertos: 544
SUBA :1230 Choque :9998 Solo Daños : 1115
PUENTE ARANDA:1068 Incendio : 2
BOSA :1044 Otro : 146
FONTIBON : 874 Volcamiento : 402
(Other) :6626
Gravedad_indicador_30d Hora_Acc Dia_Semana_Acc
Con Heridos:12466 Min. : 0.00 lunes :1865
Con Muertos: 525 1st Qu.: 8.00 martes :2018
Solo Daños : 1115 Median :13.00 miércoles:2048
Mean :12.76 jueves :1975
3rd Qu.:18.00 viernes :2266
Max. :23.00 sábado :2228
domingo :1706
Window: rectangle = [981148.8, 1007237.8] x [938934.1, 1025100.6] units
(26090 x 86170 units)
Window area = 2247990000 square units# Reescalar a km para que distancias y sigma estén en km
ppp_siniestros <- rescale(ppp_siniestros, s = 1000, unitname = c("km","kilometres"))
cat("\n✓ Patrón reescalado: unidades en km\n")
✓ Patrón reescalado: unidades en kmcat(" Área ventana:", round(area(ppp_siniestros$window), 2), "km²\n") Área ventana: 2247.99 km²cat(" Intensidad global:", round(intensity(ppp_siniestros), 2), "siniestros/km²\n") Intensidad global: 6.27 siniestros/km²# Información adicional del CRS
cat("\nInformación del sistema de coordenadas:\n")
Información del sistema de coordenadas:cat("CRS original (WGS84):", st_crs(pts_sf)$input, "\n")CRS original (WGS84): EPSG:4326 cat("CRS proyectado:", st_crs(pts_m)$input, "\n")CRS proyectado: EPSG:3116 cat("Unidades:", st_crs(pts_m)$units_gdal, "\n")Unidades: metre En esta sección se construye el patrón puntual que servirá como objeto central para el análisis con spatstat. Primero, las coordenadas originales en longitud/latitud (WGS84) se convierten a un objeto espacial sf y luego se proyectan al sistema EPSG:3116, cuya unidad es el metro. Tras la proyección, las coordenadas resultantes presentan rangos de aproximadamente 981,149 a 1,007,238 metros en el eje X y 938,934 a 1,025,101 metros en el eje Y, lo cual confirma que el patrón se encuentra en un marco consistente en unidades métricas.
A continuación, se define una ventana de observación 𝑊como un rectángulo envolvente que contiene todos los puntos proyectados. Esta ventana se usa como aproximación inicial del área donde se observa el proceso espacial. La ventana obtenida tiene un área de 2247.99 km², valor que corresponde al rectángulo mínimo que cubre todas las ubicaciones registradas; por ello, puede ser mayor que el área administrativa real de Bogotá y se interpreta como una ventana provisional, útil para comenzar el análisis exploratorio (más adelante puede reemplazarse por un polígono oficial del límite urbano o del distrito).
Con la ventana definida, se construye el objeto ppp denominado ppp_siniestros, que contiene 14106 eventos (siniestros) y sus ubicaciones. El resumen reporta un aviso de puntos duplicados: existen 307 siniestros con coordenadas idénticas, equivalente al 2.18% del total. Esto es esperable en datos de siniestralidad, ya que varios eventos pueden ocurrir en el mismo cruce o porque el geocodificador asigna la misma coordenada a una dirección repetida.
Finalmente, se calcula la intensidad global (densidad promedio de eventos por unidad de área). En este caso, la intensidad es 6.27 siniestros/km², que se obtiene como ( = n/|W| ), donde (n=14106) es el número de eventos y (|W|=2247.99) km² es el área de la ventana. Este valor es un promedio espacial: no implica que los siniestros estén uniformemente distribuidos, sino que sirve como referencia inicial antes de estudiar la inhomogeneidad (variación espacial de la intensidad) y las interacciones entre eventos.
Además, se incorporan marcas al patrón (6 variables: Localidad, Clase_Acc, dos indicadores de gravedad, Hora_Acc y Dia_Semana_Acc), lo cual permite análisis posteriores diferenciados por tipo de siniestro, gravedad y componentes temporales.
plot(ppp_siniestros,
main = "Patrón de Siniestros Viales - Bogotá 2023",
pch = 16, cex = 0.3, col = "white",
xlab = "Coordenada X (km)",
ylab = "Coordenada Y (km)")
Para mejorar la visualización sin afectar el análisis estadístico, crearemos un patrón auxiliar eliminando outliers extremos solo para gráficos.
# Detectar outliers extremos (más allá de 3 desviaciones estándar)
outliers_x <- abs(scale(x)) > 3
outliers_y <- abs(scale(y)) > 3
outliers <- outliers_x | outliers_y
cat("Outliers detectados:", sum(outliers), "de", length(x), "puntos\n")Outliers detectados: 14 de 14106 puntoscat("Porcentaje de outliers:", round(100 * sum(outliers) / length(x), 2), "%\n")Porcentaje de outliers: 0.1 %# IMPORTANTE: Crear patrón SOLO para visualización (no para análisis)
if(sum(outliers) > 0) {
x_plot <- x[!outliers]
y_plot <- y[!outliers]
siniestros_plot <- siniestros_clean[!outliers, ]
# Nueva ventana para visualización
W_plot <- owin(xrange = range(x_plot), yrange = range(y_plot))
# Patrón puntual solo para gráficos
ppp_plot <- ppp(x = x_plot, y = y_plot, window = W_plot)
marks(ppp_plot) <- siniestros_plot %>%
select(Localidad, Clase_Acc, Gravedad_Indicador_Tradicional,
Gravedad_indicador_30d, Hora_Acc, Dia_Semana_Acc)
# Reescalar a km
ppp_plot <- rescale(ppp_plot, s = 1000, unitname = c("km", "kilometres"))
cat("\n✓ Patrón para visualización creado:\n")
cat(" Puntos:", npoints(ppp_plot), "\n")
cat(" Área:", round(area(ppp_plot$window), 2), "km²\n")
cat(" Intensidad aparente:", round(intensity(ppp_plot), 2), "siniestros/km²\n")
cat("\n✓ Patrón original (para análisis) conservado:\n")
cat(" Puntos:", npoints(ppp_siniestros), "\n")
cat(" Área:", round(area(ppp_siniestros$window), 2), "km²\n")
cat(" Intensidad real:", round(intensity(ppp_siniestros), 2), "siniestros/km²\n")
} else {
cat("\nNo se detectaron outliers significativos\n")
ppp_plot <- ppp_siniestros
}
✓ Patrón para visualización creado:
Puntos: 14092
Área: 857.23 km²
Intensidad aparente: 16.44 siniestros/km²
✓ Patrón original (para análisis) conservado:
Puntos: 14106
Área: 2247.99 km²
Intensidad real: 6.27 siniestros/km²# Mapa simple y limpio
plot(ppp_plot,
pch = 20, # Punto pequeño
cex = 0.3, # Tamaño pequeño
cols = "#E31A1C", # Rojo vibrante
main = "Distribución de Siniestros Viales - Bogotá 2023",
xlab = "Coordenada X (km)",
ylab = "Coordenada Y (km)",
lwd = 0.5)
# Añadir recuadro
box(lwd = 1.5)
El gráfico del patrón puntual usando la ventana rectangular (W) puede verse poco informativo, porque unos pocos puntos extremos “estiran” el rectángulo de observación y comprimen visualmente la nube principal de eventos. Para mejorar exclusivamente la visualización, se detectan coordenadas extremas mediante un criterio simple (valores más allá de 3 desviaciones estándar en (x) o (y)).
En este caso se identificaron 14 puntos (aprox. 0.1% del total). Al excluirlos se obtiene una ventana rectangular más ajustada, lo que permite apreciar mejor la forma general del patrón y las zonas de concentración.
# Usar ppp_plot para visualización
dens_general <- density(ppp_plot, sigma = 1)
plot(dens_general,
main = "Densidad de Siniestros (Kernel, σ = 1 km)",
col = hcl.colors(100, "YlOrRd", rev = TRUE),
las = 1)
contour(dens_general,
add = TRUE,
col = "black",
lwd = 0.5,
drawlabels = FALSE)
plot(ppp_plot,
add = TRUE,
pch = 20,
cex = 0.2,
col = rgb(0, 0, 0, 0.1))
Con el objetivo de explorar la distribución espacial de la siniestralidad vial, se estimó la intensidad espacial ((s)) del patrón puntual mediante un estimador kernel. En el contexto de procesos puntuales, ((s)) representa la tasa local esperada de eventos por unidad de área alrededor de la ubicación (s); por tanto, su estimación permite identificar zonas con mayor o menor concentración de siniestros de manera descriptiva.
Sea ({x_1,,x_n}) el conjunto de ubicaciones observadas. El estimador kernel de la intensidad puede escribirse como:
[ (s)={i=1}^{n} K{}(s-x_i), K_{}(u)=K!(), ]
donde (K()) es una función kernel (típicamente gaussiana en spatstat) y () es el parámetro de suavizado. En este análisis se utilizó (= 1) km (las coordenadas del patrón están expresadas en kilómetros, i.e., Unit of length: 1 km), lo cual corresponde a una escala urbana razonable para sintetizar la variación espacial de los siniestros sin amplificar excesivamente la micro-variabilidad a nivel de intersecciones.
El resultado es una superficie ((s)) representada como mapa de calor, sobre la cual se superponen curvas de nivel (isolíneas) para resaltar gradientes y “crestas” de intensidad, y los puntos originales del patrón para contrastar la estimación suavizada con las observaciones. Visualmente, la superficie de intensidad evidencia una marcada inhomogeneidad espacial, con zonas de mayor intensidad y áreas con intensidad menor dentro de la ventana de estudio.
Es importante notar que este mapa de densidad kernel constituye una herramienta exploratoria: describe la estructura espacial promedio a la escala (), pero no establece por sí mismo significancia estadística de agrupamiento. La evaluación inferencial del patrón (p. ej., bajo CSR o modelos inhomogéneos) se realiza posteriormente mediante funciones resumen y envolventes Monte Carlo.
# Usar ppp_siniestros (patrón completo) para estadísticas
freq_localidad <- table(marks(ppp_siniestros)$Localidad) %>%
as.data.frame() %>%
rename(Localidad = Var1, Frecuencia = Freq) %>%
arrange(desc(Frecuencia)) %>%
mutate(
Porcentaje = round(100 * Frecuencia / sum(Frecuencia), 2),
Acumulado = round(cumsum(Porcentaje), 2)
)
freq_localidad %>%
head(10) %>%
kable(caption = "Top 10 Localidades con Mayor Siniestralidad") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Localidad | Frecuencia | Porcentaje | Acumulado |
|---|---|---|---|
| KENNEDY | 1969 | 13.96 | 13.96 |
| ENGATIVA | 1295 | 9.18 | 23.14 |
| SUBA | 1230 | 8.72 | 31.86 |
| PUENTE ARANDA | 1068 | 7.57 | 39.43 |
| BOSA | 1044 | 7.40 | 46.83 |
| FONTIBON | 874 | 6.20 | 53.03 |
| USAQUEN | 755 | 5.35 | 58.38 |
| CIUDAD BOLIVAR | 687 | 4.87 | 63.25 |
| TEUSAQUILLO | 632 | 4.48 | 67.73 |
| SAN CRISTOBAL | 601 | 4.26 | 71.99 |
freq_localidad %>%
head(10) %>%
ggplot(aes(x = reorder(Localidad, Frecuencia), y = Frecuencia)) +
geom_col(fill = "#E31A1C", alpha = 0.8) +
geom_text(aes(label = Frecuencia), hjust = -0.2, size = 3) +
coord_flip() +
labs(
title = "Siniestros por Localidad - Top 10",
x = "Localidad",
y = "Número de Siniestros"
) +
theme_minimal(base_size = 12) +
theme(
plot.title = element_text(face = "bold", size = 14),
panel.grid.major.y = element_blank(),
panel.grid.minor = element_blank()
) +
ylim(0, max(freq_localidad$Frecuencia) * 1.1)
Dado que el patrón puntual ppp_siniestros incluye la variable Localidad como marca (mark), se realizó un análisis descriptivo de frecuencia por categoría para identificar en qué localidades se concentra el mayor número de siniestros. Este análisis no incorpora distancias ni dependencia espacial, su propósito es resumir la composición del patrón por una variable administrativa y aportar contexto a los hallazgos espaciales.
En primer lugar, se contabilizó el número de siniestros por localidad mediante:
Los resultados muestran que la siniestralidad se concentra fuertemente en un subconjunto de localidades. En particular, el Top 10 está encabezado por Kennedy (1,969; 13.96%), seguida por Engativá (1,295; 9.18%) y Suba (1,230; 8.72%). En conjunto, las tres localidades acumulan aproximadamente 31.86% de los siniestros registrados. Al ampliar al Top 10, se alcanza cerca de 71.99% del total, lo cual evidencia una concentración marcada del conteo en estas áreas.
Este patrón de concentración por localidades debe interpretarse con cautela: las frecuencias dependen tanto del tamaño del área como de la exposición (población, volumen de tráfico, cantidad de vías, centralidades, etc.).
# Usar patrón completo para estadísticas
freq_clase <- table(marks(ppp_siniestros)$Clase_Acc) %>%
as.data.frame() %>%
rename(Clase = Var1, Frecuencia = Freq) %>%
arrange(desc(Frecuencia)) %>%
mutate(Porcentaje = round(100 * Frecuencia / sum(Frecuencia), 2))
freq_clase %>%
kable(caption = "Distribución por Clase de Accidente") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Clase | Frecuencia | Porcentaje |
|---|---|---|
| Choque | 9998 | 70.88 |
| Atropello | 2760 | 19.57 |
| Caida de ocupante | 798 | 5.66 |
| Volcamiento | 402 | 2.85 |
| Otro | 146 | 1.04 |
| Incendio | 2 | 0.01 |
# Usar ppp_plot para visualización
ppp_por_clase <- split(ppp_plot, f = marks(ppp_plot)$Clase_Acc)
plot(ppp_por_clase,
pch = 20,
cex = 0.4,
main = "",
cols = hcl.colors(length(ppp_por_clase), "Set 2")
)
La distribución por clase evidencia un fuerte desbalance en el tipo de siniestros registrados. En particular, Choque concentra la gran mayoría de los eventos con 9,998 casos (70.88%), seguido por Atropello con 2,760 casos (19.57%). A bastante distancia aparecen Caída de ocupante con 798 (5.66%) y Volcamiento con 402 (2.85%). Las categorías Otro (146; 1.04%) e Incendio (2; 0.01%) son muy poco frecuentes.
En el mapa por paneles (patrón separado por Clase_Acc), esta diferencia de tamaños se refleja directamente en la densidad visual de puntos: el panel de Choque aparece mucho más “cargado”, mientras que clases raras como Incendio muestran muy pocos eventos, por lo que cualquier aparente concentración debe interpretarse con cautela (puede deberse únicamente al bajo tamaño muestral).
Aun así, el gráfico sugiere dos conclusiones exploratorias importantes:
Dado el desbalance de frecuencias, en secciones posteriores es razonable centrar la inferencia espacial principalmente en las clases Choque y Atropello, y tratar las clases muy raras (Incendio, y en menor medida Otro) como exploratorias.
# Usar patrón completo para estadísticas
freq_gravedad <- table(marks(ppp_siniestros)$Gravedad_Indicador_Tradicional) %>%
as.data.frame() %>%
rename(Gravedad = Var1, Frecuencia = Freq) %>%
arrange(desc(Frecuencia)) %>%
mutate(Porcentaje = round(100 * Frecuencia / sum(Frecuencia), 2))
freq_gravedad %>%
kable(caption = "Distribución por Gravedad del Siniestro") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Gravedad | Frecuencia | Porcentaje |
|---|---|---|
| Con Heridos | 12447 | 88.24 |
| Solo Daños | 1115 | 7.90 |
| Con Muertos | 544 | 3.86 |
# Usar ppp_plot para visualización
ppp_por_gravedad <- split(ppp_plot,
f = marks(ppp_plot)$Gravedad_Indicador_Tradicional)
plot(ppp_por_gravedad,
pch = 20,
cex = 0.5,
main = "",
cols = c("#2ECC71", "#E74C3C", "#F39C12"))
A partir de la marca Gravedad_Indicador_Tradicional se construyó una tabla de frecuencias y porcentajes. Los resultados muestran una clara dominancia de eventos Con Heridos (12,447; 88.24%), seguidos por Solo Daños (1,115; 7.90%) y Con Muertos (544; 3.86%). Esta asimetría indica que la mayor parte de los registros corresponde a siniestros con afectación a personas, mientras que los eventos fatales representan una fracción menor del total, aunque de alto interés para priorización de intervenciones.
Para evaluar si los distintos niveles de gravedad exhiben patrones espaciales diferenciados, el patrón puntual se particionó por gravedad mediante split(ppp_plot, f = marks(...)$Gravedad_Indicador_Tradicional). Esto permite visualizar tres subpatrones (Con Heridos, Con Muertos, Solo Daños) sobre la misma ventana de observación. La visualización sugiere que los eventos de menor frecuencia (en particular, Con Muertos) presentan una distribución más dispersa y con menor densidad aparente, lo cual es consistente con su tamaño muestral reducido, mientras que Con Heridos reproduce más claramente las zonas de alta concentración observadas en el patrón total.
# Usar patrón completo para estadísticas temporales
hora_data <- marks(ppp_siniestros)$Hora_Acc
hist(hora_data,
breaks = 24,
col = viridis(24),
border = "white",
main = "Distribución de Siniestros por Hora del Día",
xlab = "Hora (0-23)",
ylab = "Frecuencia",
las = 1)
abline(h = mean(table(hora_data)), lty = 2, col = "red", lwd = 2)
hora_tabla <- table(hora_data)
hora_max <- as.numeric(names(which.max(hora_tabla)))
cat("\nHora pico:", hora_max, "horas con", max(hora_tabla), "siniestros\n")
Hora pico: 7 horas con 909 siniestrosSe estudió la marca Hora_Acc mediante un histograma con 24 barras (0–23). Adicionalmente, se trazó una línea horizontal en la frecuencia promedio por hora, ({n} = _{h=0}^{23} n_h), para identificar horas con siniestralidad por encima del nivel medio. El patrón horario muestra un comportamiento claramente no uniforme, con un pico alrededor de las 7:00 (909 siniestros), consistente con un aumento de exposición durante la hora de mayor movilidad (inicio de jornada). En general, se observa un incremento en franjas asociadas a desplazamientos laborales/estudiantiles y una disminución en horas de menor circulación.
# Usar patrón completo
dia_data <- marks(ppp_siniestros)$Dia_Semana_Acc
dia_tabla <- table(dia_data)
barplot(dia_tabla,
col = viridis(7),
border = "white",
main = "Distribución de Siniestros por Día de la Semana",
las = 2,
ylab = "Número de Siniestros")
text(x = seq_along(dia_tabla) * 1.2 - 0.5,
y = dia_tabla + 50,
labels = dia_tabla,
pos = 3,
cex = 0.9)
Finalmente, se evaluó Dia_Semana_Acc mediante un gráfico de barras con frecuencias por día. Los resultados evidencian variación sistemática a lo largo de la semana, con mayor número de siniestros hacia finales de semana y durante el fin de semana, mientras que el domingo presenta el menor conteo. Este patrón es compatible con cambios en la exposición (volumen de viajes), actividades recreativas y diferencias en comportamiento vial entre días hábiles y no hábiles.
El análisis de primer orden estudia la intensidad del proceso puntual, es decir, el número esperado de eventos por unidad de área.
# =========================================================
# 4) Análisis de Primer Orden: Intensidad
# Nota: ppp_siniestros = patrón completo (en km)
# ppp_plot = patrón para visualización (en km)
# =========================================================
# --- Intensidad global ---
intensidad_global <- intensity(ppp_plot)
cat("═══ INTENSIDAD GLOBAL ═══\n\n")═══ INTENSIDAD GLOBAL ═══cat("Intensidad:", round(intensidad_global, 2), "siniestros/km²\n")Intensidad: 16.44 siniestros/km²cat("Total de eventos:", npoints(ppp_plot), "\n")Total de eventos: 14092 cat("Área de estudio:", round(area(ppp_plot$window), 2), "km²\n\n")Área de estudio: 857.23 km²# --- Conteo por cuadrantes (10x10) + diagnóstico V/M ---
Q <- quadratcount(ppp_plot, nx = 10, ny = 10)
conteos <- as.vector(Q)
conteos <- conteos[!is.na(conteos)]
vm <- var(conteos) / mean(conteos)
cat("═══ CUADRANTES (10×10) ═══\n\n")═══ CUADRANTES (10×10) ═══cat("Conteo mínimo:", min(conteos), "\n")Conteo mínimo: 0 cat("Conteo máximo:", max(conteos), "\n")Conteo máximo: 677 cat("Conteo medio:", round(mean(conteos), 2), "\n")Conteo medio: 140.92 cat("Varianza:", round(var(conteos), 2), "\n")Varianza: 33342.74 cat("Relación Varianza/Media (V/M):", round(vm, 2), "\n\n")Relación Varianza/Media (V/M): 236.61 cat("Interpretación (primer orden):\n")Interpretación (primer orden):cat("- Bajo CSR homogéneo (Poisson con intensidad constante), se esperaría V/M ≈ 1.\n")- Bajo CSR homogéneo (Poisson con intensidad constante), se esperaría V/M ≈ 1.cat("- Un V/M muy grande sugiere fuerte INHOMOGENEIDAD (intensidad no constante).\n\n")- Un V/M muy grande sugiere fuerte INHOMOGENEIDAD (intensidad no constante).# --- Test CSR (Chi-cuadrado sobre cuadrantes) ---
test_csr <- quadrat.test(ppp_plot, nx = 10, ny = 10)
cat("═══ TEST CSR (CHI-CUADRADO) ═══\n\n")═══ TEST CSR (CHI-CUADRADO) ═══cat("X²:", round(test_csr$statistic, 2), "\n")X²: 23424.15 cat("gl:", test_csr$parameter, "\n")gl: 99 cat("p-valor:", format.pval(test_csr$p.value, digits = 4), "\n\n")p-valor: < 0.00000000000000022 if(test_csr$p.value < 0.05) {
cat("Conclusión: Se RECHAZA CSR homogéneo (p < 0.05).\n")
cat("Interpretación: la intensidad espacial NO es constante (patrón inhomogéneo).\n\n")
} else {
cat("Conclusión: No se rechaza CSR homogéneo.\n")
cat("Interpretación: no hay evidencia suficiente contra intensidad constante\n")
cat("(a esta escala de cuadrantes).\n\n")
}Conclusión: Se RECHAZA CSR homogéneo (p < 0.05).
Interpretación: la intensidad espacial NO es constante (patrón inhomogéneo).# --- Intensidad por cuadrante (mismo Q) ---
intensidades_Q <- intensity(Q)
cat("═══ INTENSIDAD POR CUADRANTE ═══\n\n")═══ INTENSIDAD POR CUADRANTE ═══cat("Intensidad mínima:", round(min(intensidades_Q, na.rm = TRUE), 2), "siniestros/km²\n")Intensidad mínima: 0 siniestros/km²cat("Intensidad máxima:", round(max(intensidades_Q, na.rm = TRUE), 2), "siniestros/km²\n")Intensidad máxima: 78.98 siniestros/km²cat("Rango:", round(max(intensidades_Q, na.rm = TRUE) - min(intensidades_Q, na.rm = TRUE), 2),
"siniestros/km²\n")Rango: 78.98 siniestros/km²A partir del patrón puntual observado en la zona urbana de Bogotá (ventana ajustada sin outliers extremos), se estimó la intensidad global como:
\[\hat{\lambda} = \frac{n}{|W|} = \frac{14{,}092}{857.23} = 16.44 \text{ siniestros/km}^2\]
donde \(n = 14{,}092\) eventos fueron observados en una ventana de área \(|W| = 857.23\) km².
Esto significa que, en promedio, en el área urbana de estudio ocurren aproximadamente 16.44 siniestros por cada km². Sin embargo, este valor es solo un promedio espacial y no implica que la intensidad sea uniforme en toda la zona urbana.
Al dividir la ventana en una rejilla de 10×10 cuadrantes, los conteos por celda muestran una variabilidad muy alta:
La razón Varianza/Media resulta:
\[V/M = \frac{\text{Var}(N_i)}{\mathbb{E}[N_i]} = \frac{33{,}342.74}{140.92} = 236.61\]
Interpretación del ratio V/M:
Bajo un proceso de Poisson homogéneo (CSR con intensidad constante \(\lambda\)), los conteos por cuadrante siguen una distribución Poisson y se esperaría aproximadamente \(V/M \approx 1\). El valor obtenido (\(V/M = 236.61\)) es enormemente mayor que 1, lo cual indica que los conteos por cuadrante son mucho más variables de lo que se esperaría bajo homogeneidad.
En términos de primer orden, esto sugiere con fuerza que la intensidad no es constante en el espacio. El patrón es claramente inhomogéneo: existen zonas con alta concentración de siniestros y otras con muy pocos o ningún evento.
Los 34 cuadrantes sin eventos (34%) corresponden principalmente a áreas con baja o nula cobertura vial dentro de la ventana (parques urbanos, zonas industriales, cerros orientales), lo cual contribuye a inflar el ratio \(V/M\).
Para evaluar si esta inhomogeneidad es un artefacto de la definición de la ventana o refleja variación real en el patrón urbano, se realizó una prueba adicional utilizando un convex hull (envolvente convexa) como ventana alternativa. El convex hull elimina automáticamente áreas periféricas y genera una ventana más ajustada (538 km²) que reduce los cuadrantes vacíos de 34 a 12 (-65%).
Resultado: Incluso con convex hull, el ratio \(V/M\) solo se redujo a 197.02 (-17%), permaneciendo muy por encima de 1. Esto confirma que la inhomogeneidad es genuina y no un artefacto de la ventana: refleja diferencias reales en la intensidad de siniestralidad entre zonas urbanas de Bogotá (por ejemplo, Kennedy con λ > 40 siniestros/km² versus zonas residenciales con λ < 10 siniestros/km²).
El test \(\chi^2\) de bondad de ajuste aplicado a los cuadrantes evalúa formalmente la hipótesis nula:
\[H_0: \lambda(s) = \lambda \text{ (constante) para todo } s \in W\]
Resultados:
Conclusión estadística:
Se rechaza contundentemente la hipótesis de CSR homogéneo (\(p < 0.05\)). La evidencia es abrumadora contra la hipótesis de intensidad constante.
Esta evidencia es consistente con el diagnóstico \(V/M\): existe un componente fuerte de primer orden (inhomogeneidad), es decir, la intensidad \(\lambda(s)\) varía espacialmente de forma sistemática.
Al convertir los conteos por cuadrante en intensidades locales \(\hat{\lambda}_i = N_i / |A_i|\), donde \(|A_i|\) es el área del cuadrante \(i\), se observa un rango muy amplio:
Interpretación:
Esta diferencia muestra que existen cuadrantes con casi ausencia de eventos y otros con una intensidad muy superior al promedio global (16.44 siniestros/km²). En particular, algunos cuadrantes alcanzan intensidades de hasta 4.8 veces el promedio urbano, lo cual confirma la existencia de hot spots claramente diferenciados.
Esta heterogeneidad espacial marcada refuerza la conclusión de que el proceso generador de siniestros viales no es uniforme en el espacio urbano, sino que está fuertemente modulado por factores espaciales como la densidad de tráfico, el tipo de vías, la centralidad urbana y las características locales de cada zona.
En conjunto, los resultados indican que:
El rechazo de CSR homogéneo tiene importantes consecuencias metodológicas para el análisis subsecuente:
Pregunta pendiente:
Dado que se confirmó inhomogeneidad (\(\lambda(s)\) no constante), la pregunta de segundo orden es: ¿Existe interacción espacial entre eventos más allá de la variación en intensidad?
Es decir, una vez “descontada” la inhomogeneidad de primer orden, ¿los eventos muestran:
El análisis de segundo orden estudia la dependencia espacial entre eventos, es decir, si los siniestros tienden a ocurrir cerca unos de otros (agregación) o se evitan mutuamente (regularidad), más allá de la variación en intensidad ya detectada.
Dado que confirmamos inhomogeneidad (\(\lambda(s)\) no constante), usaremos funciones de segundo orden con corrección por inhomogeneidad.
Primero estimamos la intensidad local mediante densidad kernel, que será usada para normalizar las funciones de segundo orden:
# =========================================================
# Estimación de intensidad espacial λ(s) con kernel
# =========================================================
cat("═══ ESTIMACIÓN DE INTENSIDAD λ(s) ═══\n\n")═══ ESTIMACIÓN DE INTENSIDAD λ(s) ═══# Ancho de banda óptimo (regla de Scott para datos espaciales)
sigma_bw <- bw.scott(ppp_plot)
cat("Ancho de banda (σ):\n")Ancho de banda (σ):cat(" Método: Scott (isotropic)\n") Método: Scott (isotropic)cat(" σ_x =", round(sigma_bw[1], 3), "km\n") σ_x = 0.894 kmcat(" σ_y =", round(sigma_bw[2], 3), "km\n") σ_y = 1.327 kmcat(" σ promedio =", round(mean(sigma_bw), 3), "km\n\n") σ promedio = 1.11 km# Calcular densidad kernel (intensidad estimada)
lambda_hat <- density(ppp_plot, sigma = sigma_bw)
cat("✓ Intensidad λ(s) estimada\n")✓ Intensidad λ(s) estimadacat(" Intensidad mínima:", round(min(lambda_hat), 2), "siniestros/km²\n") Intensidad mínima: 0 siniestros/km²cat(" Intensidad máxima:", round(max(lambda_hat), 2), "siniestros/km²\n") Intensidad máxima: 81.06 siniestros/km²cat(" Intensidad media:", round(mean(lambda_hat), 2), "siniestros/km²\n\n") Intensidad media: 16.48 siniestros/km²# Visualizar intensidad estimada
plot(lambda_hat,
main = expression(paste("Intensidad Espacial Estimada ", lambda, "(s)")),
col = hcl.colors(100, "YlOrRd", rev = TRUE),
las = 1)
# Añadir contornos
contour(lambda_hat, add = TRUE, col = "black", lwd = 0.8, drawlabels = FALSE)
# Añadir puntos (semi-transparentes)
plot(ppp_plot, add = TRUE, pch = 20, cex = 0.2, col = rgb(0,0,0,0.15))
Este mapa muestra la intensidad local estimada mediante suavizado kernel. Las zonas rojas (alta λ) son áreas con mayor tasa esperada de siniestros, mientras que las zonas amarillas/blancas tienen menor intensidad. Esta superficie será usada para normalizar las funciones de segundo orden.
La función \(K_{\text{inhom}}(r)\) evalúa si los eventos están más agregados o más regulares de lo esperado bajo un proceso de Poisson inhomogéneo (independencia espacial con λ(s) variable).
# =========================================================
# Función K inhomogénea
# =========================================================
cat("═══ FUNCIÓN K INHOMOGÉNEA ═══\n\n")═══ FUNCIÓN K INHOMOGÉNEA ═══# Calcular K inhomogénea con corrección de borde
K_inhom <- Kinhom(ppp_plot,
lambda = lambda_hat,
correction = "border",
r = seq(0, 5, by = 0.1)) # Hasta 5 km
cat("✓ K_inhom calculada\n\n")✓ K_inhom calculada# Transformación L para estabilizar varianza
L_inhom <- Linhom(ppp_plot,
lambda = lambda_hat,
correction = "border",
r = seq(0, 5, by = 0.1))
cat("✓ L_inhom calculada\n\n")✓ L_inhom calculada# =========================================================
# Solución: Forzar lambda_hat a ser estrictamente positivo
# =========================================================
cat("═══ CORRECCIÓN DE λ(s) PARA ENVOLVENTES ═══\n\n")═══ CORRECCIÓN DE λ(s) PARA ENVOLVENTES ═══# Clonar lambda_hat y forzar valores positivos
lambda_hat_pos <- lambda_hat
lambda_hat_pos[lambda_hat_pos <= 0] <- min(lambda_hat[lambda_hat > 0])
cat("✓ λ(s) corregida para ser estrictamente positiva\n")✓ λ(s) corregida para ser estrictamente positivacat(" Valores negativos o cero reemplazados por:",
round(min(lambda_hat[lambda_hat > 0]), 4), "\n\n") Valores negativos o cero reemplazados por: 0 # Envolventes Monte Carlo con lambda corregida
cat("Generando envolventes Monte Carlo... (puede tardar 5-10 minutos)\n")Generando envolventes Monte Carlo... (puede tardar 5-10 minutos)cat(" Simulaciones: 99 réplicas bajo Poisson inhomogéneo\n\n") Simulaciones: 99 réplicas bajo Poisson inhomogéneoset.seed(42)
env_K <- envelope(ppp_plot,
fun = Kinhom,
lambda = lambda_hat_pos, # Usar lambda corregida
nsim = 99,
correction = "border",
r = seq(0, 5, by = 0.1),
verbose = FALSE,
savefuns = FALSE) # No guardar funciones individuales
cat("✓ Envolventes K calculadas\n\n")✓ Envolventes K calculadasenv_L <- envelope(ppp_plot,
fun = Linhom,
lambda = lambda_hat_pos, # Usar lambda corregida
nsim = 99,
correction = "border",
r = seq(0, 5, by = 0.1),
verbose = FALSE,
savefuns = FALSE)
cat("✓ Envolventes L calculadas\n\n")✓ Envolventes L calculadaspar(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
# Panel 1: K inhomogénea
plot(env_K,
main = expression(paste("Función ", K[inhom], "(r)")),
xlab = "Distancia r (km)",
ylab = expression(K[inhom](r)),
legend = FALSE,
lwd = 2)
# Línea teórica (Poisson inhomogéneo)
lines(env_K$r, pi * env_K$r^2, lty = 2, col = "red", lwd = 2)
legend("topleft",
legend = c("Observado", "Teórico (Poisson inhom)", "Envolvente 99%"),
lty = c(1, 2, 1),
col = c("black", "red", "grey"),
lwd = c(2, 2, 8),
bty = "n")
# Panel 2: L - r
plot(env_L,
. - r ~ r,
main = expression(paste("Función ", L[inhom], "(r) - r")),
xlab = "Distancia r (km)",
ylab = expression(L[inhom](r) - r),
legend = FALSE,
lwd = 2)
abline(h = 0, lty = 2, col = "red", lwd = 2)
legend("topleft",
legend = c("Observado", "Independencia", "Envolvente 99%"),
lty = c(1, 2, 1),
col = c("black", "red", "grey"),
lwd = c(2, 2, 8),
bty = "n")
par(mfrow = c(1, 1))En esta sección se evalúa la interacción espacial de segundo orden controlando la inhomogeneidad de primer orden (es decir, permitiendo que la intensidad ((s)) varíe en el espacio). Para ello se usan las funciones (K_{}(r)) y (L_{}(r)-r), junto con envolventes de Monte Carlo (99%) bajo un Poisson inhomogéneo.
Observaciones:
Conclusión (K):
No se detecta evidencia estadísticamente significativa de agregación ni de regularidad. El patrón es compatible con un proceso de Poisson inhomogéneo, es decir, eventos independientes con intensidad espacial variable.
Observaciones:
Conclusión (L):
No hay evidencia significativa de agregación a ninguna escala espacial considerada (0–5 km). Una vez controlada la variación de ((s)), el patrón resulta consistente con independencia espacial.
Resultado principal:
Una vez controlada la inhomogeneidad de la intensidad ((s)), no se observa interacción espacial significativa entre los siniestros.
Esto implica que:
Los hot spots no sugieren un “contagio espacial” entre siniestros ni clusters auto-organizados. Más bien reflejan condiciones del entorno que elevan la intensidad local ((s)), por ejemplo: