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")

6. Conclusiones Finales

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.