Este documento presenta un análisis detallado de los accidentes automovilísticos registrados en los municipios del estado de Nuevo León. A través de técnicas de análisis espacial, modelos de regresión global y local, y herramientas visuales, se busca identificar patrones, clústeres y factores asociados a la incidencia de estos eventos.
El análisis se estructura en etapas: limpieza e importación de datos, análisis exploratorio descriptivo, evaluación de la autocorrelación espacial, ajuste de modelos predictivos globales (LM, SAR, SEM, SDM) y, finalmente, un modelo de regresión geográficamente ponderado (GWR). Todo el proceso se desarrolla desde un enfoque técnico riguroso, con interpretaciones que permiten traducir los hallazgos en recomendaciones para la formulación de políticas públicas.
Antes de comenzar con el procesamiento, se eliminan todos los objetos del entorno de R para garantizar una ejecución limpia del script.
Se cargan todas las librerías necesarias para el análisis espacial, asegurando que cada una esté instalada previamente.
packages <- c("sf", "dplyr", "ggplot2", "tmap", "spdep", "spatialreg", "spgwr", "GWmodel", "lmtest", "car", "tmaptools", "viridis")
lapply(packages, function(pkg) {
if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
library(pkg, character.only = TRUE)
})
## [[1]]
## [1] "sf" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "dplyr" "sf" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[3]]
## [1] "ggplot2" "dplyr" "sf" "stats" "graphics" "grDevices"
## [7] "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "tmap" "ggplot2" "dplyr" "sf" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "spdep" "spData" "tmap" "ggplot2" "dplyr" "sf"
## [7] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [13] "base"
##
## [[6]]
## [1] "spatialreg" "Matrix" "spdep" "spData" "tmap"
## [6] "ggplot2" "dplyr" "sf" "stats" "graphics"
## [11] "grDevices" "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "spgwr" "sp" "spatialreg" "Matrix" "spdep"
## [6] "spData" "tmap" "ggplot2" "dplyr" "sf"
## [11] "stats" "graphics" "grDevices" "utils" "datasets"
## [16] "methods" "base"
##
## [[8]]
## [1] "GWmodel" "Rcpp" "robustbase" "spgwr" "sp"
## [6] "spatialreg" "Matrix" "spdep" "spData" "tmap"
## [11] "ggplot2" "dplyr" "sf" "stats" "graphics"
## [16] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "lmtest" "zoo" "GWmodel" "Rcpp" "robustbase"
## [6] "spgwr" "sp" "spatialreg" "Matrix" "spdep"
## [11] "spData" "tmap" "ggplot2" "dplyr" "sf"
## [16] "stats" "graphics" "grDevices" "utils" "datasets"
## [21] "methods" "base"
##
## [[10]]
## [1] "car" "carData" "lmtest" "zoo" "GWmodel"
## [6] "Rcpp" "robustbase" "spgwr" "sp" "spatialreg"
## [11] "Matrix" "spdep" "spData" "tmap" "ggplot2"
## [16] "dplyr" "sf" "stats" "graphics" "grDevices"
## [21] "utils" "datasets" "methods" "base"
##
## [[11]]
## [1] "tmaptools" "car" "carData" "lmtest" "zoo"
## [6] "GWmodel" "Rcpp" "robustbase" "spgwr" "sp"
## [11] "spatialreg" "Matrix" "spdep" "spData" "tmap"
## [16] "ggplot2" "dplyr" "sf" "stats" "graphics"
## [21] "grDevices" "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "viridis" "viridisLite" "tmaptools" "car" "carData"
## [6] "lmtest" "zoo" "GWmodel" "Rcpp" "robustbase"
## [11] "spgwr" "sp" "spatialreg" "Matrix" "spdep"
## [16] "spData" "tmap" "ggplot2" "dplyr" "sf"
## [21] "stats" "graphics" "grDevices" "utils" "datasets"
## [26] "methods" "base"
A continuación, se importa un archivo shapefile que contiene información sobre los municipios de Nuevo León y su tasa de accidentes automovilísticos. Se imprime un resumen general para comprender la estructura del conjunto de datos.
data_sf <- st_read("/Users/gabrielmedina/Downloads/Quiz_3-2/nl_map/nl_mpios_auto_acc.shp", stringsAsFactors = FALSE)
## Reading layer `nl_mpios_auto_acc' from data source
## `/Users/gabrielmedina/Downloads/Quiz_3-2/nl_map/nl_mpios_auto_acc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 22 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -101.2068 ymin: 23.16268 xmax: -98.42158 ymax: 27.79914
## Geodetic CRS: WGS 84
summary(data_sf)
## OBJECTID CODELAG CVE_ENT IDUNICO Shape_Leng
## Min. : 1.0 Min. :2199 Min. :19 Min. :19001 Min. :0.2858
## 1st Qu.:13.5 1st Qu.:2212 1st Qu.:19 1st Qu.:19014 1st Qu.:0.8586
## Median :26.0 Median :2224 Median :19 Median :19026 Median :1.5393
## Mean :26.0 Mean :2224 Mean :19 Mean :19026 Mean :1.7867
## 3rd Qu.:38.5 3rd Qu.:2236 3rd Qu.:19 3rd Qu.:19038 3rd Qu.:2.4495
## Max. :51.0 Max. :2249 Max. :19 Max. :19051 Max. :4.7800
## Shape_Area OBJECTID_1 IDUNICO_1 mpio
## Min. :0.004224 Min. : 1.0 Min. :19001 Length:51
## 1st Qu.:0.021224 1st Qu.:13.5 1st Qu.:19014 Class :character
## Median :0.064960 Median :26.0 Median :19026 Mode :character
## Mean :0.113032 Mean :26.0 Mean :19026
## 3rd Qu.:0.145082 3rd Qu.:38.5 3rd Qu.:19038
## Max. :0.630891 Max. :51.0 Max. :19051
## auto_accid tasa_auto_ zona_urb zona_subur
## Min. : 1.0 Min. : 0.15 Length:51 Length:51
## 1st Qu.: 34.0 1st Qu.: 31.66 Class :character Class :character
## Median : 172.0 Median : 73.32 Mode :character Mode :character
## Mean : 2389.1 Mean : 2737.38
## 3rd Qu.: 624.5 3rd Qu.: 299.64
## Max. :42956.0 Max. :98483.48
## sexo aliento cinturon edad
## Min. :1.000 Length:51 Length:51 Min. :15.38
## 1st Qu.:2.000 Class :character Class :character 1st Qu.:32.10
## Median :2.000 Mode :character Mode :character Median :36.35
## Mean :1.804 Mean :37.61
## 3rd Qu.:2.000 3rd Qu.:41.46
## Max. :2.000 Max. :64.50
## pop densidad_p gini region
## Min. : 1071 Min. : 0.830 Min. :0.3000 Length:51
## 1st Qu.: 4182 1st Qu.: 4.245 1st Qu.:0.3450 Class :character
## Median : 15902 Median : 7.070 Median :0.3800 Mode :character
## Mean : 110003 Mean : 498.641 Mean :0.3922
## 3rd Qu.: 78173 3rd Qu.: 268.535 3rd Qu.:0.4400
## Max. :1124835 Max. :4748.840 Max. :0.4900
## grado_educ geometry
## Min. : 6.590 POLYGON :51
## 1st Qu.: 7.945 epsg:4326 : 0
## Median : 8.840 +proj=long...: 0
## Mean : 8.933
## 3rd Qu.: 9.645
## Max. :13.160
Se realiza una transformación inicial del dataset. Se renombran columnas para mayor claridad, se eliminan columnas técnicas que no aportan al análisis y se convierten algunas variables a tipo factor.
data_sf <- data_sf %>%
rename(
auto_accidentes = auto_accid,
tasa_auto_acc = tasa_auto_,
densidad_p = densidad_p
) %>%
select(-c(OBJECTID, OBJECTID_1, CODELAG, CVE_ENT, IDUNICO_1, Shape_Leng, Shape_Area)) %>%
mutate(across(c(region, zona_urb, zona_subur, aliento, cinturon, sexo), factor))
Se calculan estadísticas descriptivas (media y mediana) para el total estatal y por región, con el fin de explorar la distribución de los accidentes en el territorio.
state_stats <- data_sf %>% st_set_geometry(NULL) %>% summarise(media = mean(auto_accidentes), mediana = median(auto_accidentes))
region_stats <- data_sf %>% st_set_geometry(NULL) %>% group_by(region) %>% summarise(media = mean(auto_accidentes), mediana = median(auto_accidentes))
state_stats
## media mediana
## 1 2389.098 172
region_stats
## # A tibble: 5 × 3
## region media mediana
## <fct> <dbl> <dbl>
## 1 Este 198. 152.
## 2 Norte 965. 68.5
## 3 Oeste 2055. 340
## 4 Sur 79.9 70.5
## 5 ZMM 8551. 3825
Se visualiza la dispersión de accidentes a través de boxplots, permitiendo identificar la variabilidad y presencia de valores atípicos en cada región.
ggplot(data_sf, aes(x = region, y = auto_accidentes)) +
geom_boxplot(fill = viridis(1), alpha = 0.7) +
labs(title = "Dispersión de Accidentes por Región", x = "Región", y = "Accidentes por 10,000 hab.") +
theme_minimal()
Se comparan los niveles de accidentes en función del uso del cinturón y el sexo del conductor, explorando si existen diferencias importantes en la siniestralidad.
ggplot(data_sf, aes(x = cinturon, y = auto_accidentes)) +
geom_boxplot() +
labs(title = "Accidentes vs Uso de Cinturón", x = "Uso de cinturón", y = "Accidentes por 10,000 hab.") +
theme_minimal()
ggplot(data_sf, aes(x = sexo, y = auto_accidentes)) +
geom_boxplot() +
labs(title = "Accidentes vs Sexo", x = "Sexo (1=Hombre, 2=Mujer)", y = "Accidentes por 10,000 hab.") +
theme_minimal()
Se genera una matriz de contigüidad tipo Rook para detectar vecindades entre municipios. Esto será útil para evaluar dependencia espacial.
nb_rook <- poly2nb(data_sf, queen = FALSE)
listw_rook <- nb2listw(nb_rook, style = "W", zero.policy = TRUE)
plot(st_geometry(data_sf), col = "lightgrey", main = "Conectividad tipo Rook")
coords <- st_coordinates(st_centroid(data_sf))
## Warning: st_centroid assumes attributes are constant over geometries
plot(nb_rook, coords, add = TRUE, col = "red")
Se aplica la prueba de Moran Global para evaluar la presencia de autocorrelación espacial en la tasa de accidentes.
moran_global <- moran.test(data_sf$tasa_auto_acc, listw_rook, zero.policy = TRUE)
moran_global
##
## Moran I test under randomisation
##
## data: data_sf$tasa_auto_acc
## weights: listw_rook
##
## Moran I statistic standard deviate = 0.14327, p-value = 0.443
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.01516424 -0.02000000 0.00113921
Se realiza el análisis LISA (Moran Local) para detectar agrupamientos locales significativos de tasas altas o bajas de accidentes. Se clasifican los municipios según el tipo de clúster.
lisa <- localmoran(data_sf$tasa_auto_acc, listw_rook)
data_sf$lisa_I <- lisa[, 1]
data_sf$lisa_p <- lisa[, 5]
mean_rate <- mean(data_sf$tasa_auto_acc)
data_sf <- data_sf %>%
mutate(cluster = case_when(
tasa_auto_acc > mean_rate & lisa_I > 0 ~ "High-High",
tasa_auto_acc < mean_rate & lisa_I > 0 ~ "Low-Low",
tasa_auto_acc > mean_rate & lisa_I < 0 ~ "High-Low",
TRUE ~ "Low-High"
))
tm_shape(data_sf) +
tm_fill("cluster", palette = c("red", "blue", "orange", "lightblue"),
title = "Clústeres LISA") +
tm_borders() +
tm_layout(title = "Clústeres locales de accidentes (LISA)")
Se ajustan modelos de regresión global: un modelo lineal clásico (LM), un modelo espacial con rezago (SAR), uno con error espacial (SEM) y un modelo Durbin Espacial (SDM).
form <- log(tasa_auto_acc) ~ zona_urb + densidad_p + region + edad + I(edad^2) + sexo
lm_model <- lm(form, data = data_sf)
summary(lm_model)
##
## Call:
## lm(formula = form, data = data_sf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4440 -0.6831 -0.0673 0.5823 5.4782
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.1457434 4.5174960 -0.696 0.4902
## zona_urbno_interseccion 2.1542260 1.1543210 1.866 0.0694 .
## zona_urbsuburbana 1.1717921 1.4469273 0.810 0.4228
## densidad_p 0.0002395 0.0004104 0.584 0.5628
## regionNorte -0.1875490 1.0322058 -0.182 0.8567
## regionOeste -1.0841957 1.0442633 -1.038 0.3054
## regionSur -2.6271709 1.0063100 -2.611 0.0127 *
## regionZMM -1.7246836 1.4008889 -1.231 0.2255
## edad 0.3973415 0.2009280 1.978 0.0549 .
## I(edad^2) -0.0048484 0.0023373 -2.074 0.0445 *
## sexo2 1.0935001 0.7360987 1.486 0.1452
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.985 on 40 degrees of freedom
## Multiple R-squared: 0.3277, Adjusted R-squared: 0.1596
## F-statistic: 1.95 on 10 and 40 DF, p-value: 0.06633
sar_model <- lagsarlm(form, data = data_sf, listw = listw_rook, zero.policy = TRUE)
sem_model <- errorsarlm(form, data = data_sf, listw = listw_rook, zero.policy = TRUE)
sdm_model <- lagsarlm(form, data = data_sf, listw = listw_rook, type = "mixed", zero.policy = TRUE)
Se comparan los modelos con base en el criterio AIC. El modelo con menor AIC será considerado el mejor ajuste estadístico.
AIC(lm_model, sar_model, sem_model, sdm_model)
## df AIC
## lm_model 12 226.2857
## sar_model 13 228.0810
## sem_model 13 228.1449
## sdm_model 23 223.8674
Se evalúa si los residuos del modelo seleccionado (SDM) presentan autocorrelación espacial. La ausencia de autocorrelación sugiere un buen ajuste.
residuos_sdm <- residuals(sdm_model)
moran.test(residuos_sdm, listw_rook, zero.policy = TRUE)
##
## Moran I test under randomisation
##
## data: residuos_sdm
## weights: listw_rook
##
## Moran I statistic standard deviate = -0.70793, p-value = 0.7605
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.085332786 -0.020000000 0.008517018
Se visualiza el mapa de predicciones del modelo SDM para analizar su distribución geográfica.
data_sf$pred_sdm <- predict(sdm_model)
## This method assumes the response is known - see manual page
tm_shape(data_sf) +
tm_fill("pred_sdm", palette = "viridis", style = "quantile", title = "Predicciones SDM") +
tm_borders() +
tm_layout(title = "Modelo SDM – Predicciones por Municipio")
Se ajusta un modelo GWR con variables numéricas para evaluar cómo cambian los coeficientes según el espacio geográfico.
data_sp <- as(data_sf, "Spatial")
gwr_formula <- log(tasa_auto_acc) ~ densidad_p + edad
bw <- bw.gwr(gwr_formula, data = data_sp, adapt = TRUE, approach = "AIC")
## Adaptive bandwidth (number of nearest neighbours): 39 AICc value: 216.5233
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: 216.5141
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 219.1795
## Adaptive bandwidth (number of nearest neighbours): 34 AICc value: 216.3262
## Adaptive bandwidth (number of nearest neighbours): 37 AICc value: 216.1562
## Adaptive bandwidth (number of nearest neighbours): 37 AICc value: 216.1562
adapt_prop <- bw / nrow(data_sp)
gwr_model <- gwr(gwr_formula, data = data_sp, adapt = adapt_prop)
data_sf$pred_gwr <- gwr_model$SDF$pred
data_sf$local_R2 <- gwr_model$SDF$localR2
Se visualizan los resultados del GWR: predicciones y niveles de R² local para evaluar la variación espacial del ajuste.
tm_shape(data_sf) +
tm_fill("pred_gwr", palette = "viridis", style = "quantile", title = "Predicción GWR") +
tm_borders() +
tm_layout(title = "GWR – Predicción Local")
tm_shape(data_sf) +
tm_fill("local_R2", palette = "viridis", style = "quantile", title = "R² Local GWR") +
tm_borders() +
tm_layout(title = "GWR – Ajuste Local")
Este análisis permitió detectar que la siniestralidad vial en Nuevo León presenta una fuerte concentración en la zona metropolitana, revelando una distribución territorial altamente desigual. A través del índice de Moran y los clústeres LISA se identificaron zonas críticas que deben ser prioridad en la planeación de seguridad vial.
Los modelos espaciales, especialmente el SDM, mostraron ser superiores al modelo clásico, ya que permiten incorporar el efecto de los municipios vecinos. A través del SDM se reconoció que variables como la edad y la densidad poblacional no solo influyen de forma directa, sino también a través de los efectos espaciales indirectos.
El modelo GWR fue fundamental para evidenciar que la relación entre las variables explicativas y la tasa de accidentes varía geográficamente, lo que refuerza la necesidad de implementar políticas diferenciadas por región.
Finalmente, se recomienda extender este análisis en futuras investigaciones incorporando elementos temporales y de comportamiento, así como datos sobre infraestructura vial y fiscalización del tránsito. El presente estudio sienta las bases para intervenciones más focalizadas y efectivas.