Evidencia: Análisis Espacial y Predictivo de la Movilidad y Accidentes Viales en el Transporte Público de Nuevo León

Instalar paquetes y llamar librerías

library(sf)          
library(spdep)       
library(ggplot2)     
library(dplyr)      
library(leaflet)     
library(htmltools)
library(lubridate)
library(dplyr)
library(sf)
library(readxl)
library(foreign)
library(tigris)
library(rgeoda)
library(RColorBrewer)
library(viridis)
library(ggplot2)
library(tmap)
library(sf) 
library(sp)
library(spatialreg)
library(stargazer)
library(plm)
library(splm)
library(pspatreg)
library(regclass)
library(mctest)
library(lmtest)
library(spData)
library(mapview)
library(naniar)
library(dlookr)
library(caret)
library(e1071)
library(SparseM)
library(Metrics)
library(randomForest)
library(rpart.plot)
library(insight)
library(jtools)
library(xgboost)
library(DiagrammeR)
library(effects)
library(sf)
library(purrr)
library(geosphere)
library(gmapsdistance)
library(dplyr)
library(stringr)
library(knitr)
library(leaflet.extras)
library(tidyr)
library(prophet)
if (!require("pacman")) install.packages("pacman")
pacman::p_load(
  sf,          # Datos espaciales
  spdep,       # Autocorrelación espacial
  ggplot2,     # Visualización
  dbscan,      # Detección de clústers
  dplyr        # Manipulación de datos
)
if (!require("geosphere")) install.packages("geosphere")
if (!require("gmapsdistance"))install.packages("gmapsdistance")

Importar las bases de datos

datos <- read_excel("C:\\Users\\mari0\\OneDrive\\Documents\\R Studio\\8 semestre\\Evidencia\\BASE MUNICIPAL_ACCIDENTES DE TRANSITO GEORREFERENCIADOS_2023.xlsx")

geojson_data <- st_read("Transporte.geojson") %>% 
  st_transform(32614) %>%  # Convertir a UTM Zone 14N (Monterrey)
  filter(!st_is_empty(.))  # Eliminar geometrías vacías

## Reading layer `Transporte' from data source 
##   `C:\Users\mari0\OneDrive\Documents\R Studio\8 semestre\Evidencia\Transporte.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 74329 features and 4 fields
## Geometry type: LINESTRING
## Dimension:     XY
## Bounding box:  xmin: -100.7508 ymin: 25.3483 xmax: -99.7609 ymax: 26.13538
## Geodetic CRS:  WGS 84

Análisis Datos Tránsito en Nuevo León

La información de transporte proviene de la encuesta origen-destino del 2019 que el gobierno del estado de Nuevo León realizó como parte del Programa Integral de Movilidad Urbana Sustentable (PIMUS) del Área Metropolitana de Monterrey. Transconsult fue la empresa consultora responsable del levantamiento de la encuesta.

Hot Spots significativos de congestionamiento vehicular dentro del área metropolitana de Monterrey

# Desactivar geometría esférica (evita errores)
sf::sf_use_s2(FALSE)

## Spherical geometry (s2) switched off

# Crear variable numérica para análisis
geojson_data <- geojson_data %>%
  mutate(
    peso_transporte = case_when(
      Transporte %in% c("Autómovil", "TPUB", "Taxi") ~ 3,
      Transporte %in% c("Transporte de Personal", "Transporte escolar") ~ 2,
      TRUE ~ 1  # Para bicicletas, caminando y otros
    )
  )


# Obtener coordenadas de los centroides
coords <- st_coordinates(st_centroid(geojson_data))

## Warning: st_centroid assumes attributes are constant over geometries

# Calcular vecinos en un radio de 800 metros
nb <- dnearneigh(coords, d1 = 0, d2 = 800)

# Identificar y eliminar puntos aislados (sin vecinos)
aislados <- which(card(nb) == 0)
if(length(aislados) > 0){
  geojson_data <- geojson_data[-aislados, ]
  coords <- st_coordinates(st_centroid(geojson_data))
  nb <- dnearneigh(coords, d1 = 0, d2 = 800)
}

## Warning: st_centroid assumes attributes are constant over geometries

# Crear matriz de pesos espaciales
pesos <- nb2listw(nb, style = "B")  # Pesos binarios

# Calcular hotspots (Getis-Ord Gi*)
geojson_data$gi_score <- localG(
  x = geojson_data$peso_transporte,
  listw = pesos,
  zero.policy = FALSE
)

# Visualización de resultados

# Convertir LINESTRING a puntos (extremos de cada línea)
puntos_hotspots <- st_cast(geojson_data, "POINT") %>% 
  st_transform(4326)  # Convertir a WGS84

## Warning in st_cast.sf(geojson_data, "POINT"): repeating attributes for all
## sub-geometries for which they may not be constant

# Convertir los scores a numérico
geojson_data$gi_score_numeric <- as.numeric(geojson_data$gi_score)

# Convertir a puntos y filtrar
puntos_significativos <- geojson_data %>% 
  st_cast("POINT") %>% 
  st_transform(4326) %>% 
  filter(gi_score_numeric > quantile(gi_score_numeric, 0.97))  # Top 3%

## Warning in st_cast.sf(., "POINT"): repeating attributes for all sub-geometries
## for which they may not be constant

# Paleta de color para valores extremos
pal <- colorNumeric(
  palette = "Reds",
  domain = c(2.5, max(puntos_significativos$gi_score_numeric, na.rm = TRUE)),
  na.color = "transparent"
)

# Mapa final
leaflet(puntos_significativos) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    radius = ~sqrt(gi_score_numeric)*2,  # Tamaño proporcional al score
    color = ~pal(gi_score_numeric),
    fillOpacity = 0.7,
    stroke = FALSE,
    popup = ~paste("Score:", round(gi_score_numeric, 2))
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal,
    values = ~gi_score_numeric,
    title = "Hotspots<br>(Top 3%)"
  ) %>%
  setView(lng = -100.3, lat = 25.67, zoom = 12)

# Ver distribución de scores
ggplot(geojson_data, aes(x = gi_score_numeric)) + 
  geom_histogram(bins = 50) +
  geom_vline(xintercept = quantile(geojson_data$gi_score_numeric, 0.95), color = "red") +
  labs(title = "Distribución de Gi* Scores", subtitle = "Línea roja = percentil 95%")

La distribución es asimétrica con sesgo negativo: la mayor parte de los valores Gi* se concentran entre -10 y +8, con una larga cola hacia la izquierda (valores muy negativos).

Esto indica que la mayoría de los trayectos no presentan concentración significativa (Gi* cercano a 0), pero sí hay algunos con clara dispersión (Gi* muy negativo).

La línea roja marca el percentil 95, es decir, el umbral a partir del cual se considera que un Gi* score está en el 5% superior, lo cual indica un hotspot estadísticamente significativo.

# Convertir HoraOri a formato hora y luego a horas decimales
puntos_hora <- geojson_data %>%
  st_cast("POINT") %>%
  st_transform(4326) %>%
  mutate(
    hora_decimal = period_to_seconds(hms(HoraOri)) / 3600,  # Convertir a horas con decimales
    franja_horaria = cut(
      hora_decimal,
      breaks = c(0, 6, 9, 12, 15, 18, 21, 24),
      labels = c("Madrugada", "Mañana-Pico", "Mediodía", "Tarde", "Tarde-Pico", "Noche-Temprana", "Noche-Tardia"),
      right = FALSE,
      include.lowest = TRUE
    )
  ) %>%
  filter(!is.na(hora_decimal))  # Eliminar registros sin hora válida

## Warning in st_cast.sf(., "POINT"): repeating attributes for all sub-geometries
## for which they may not be constant

# Verificar
summary(puntos_hora$hora_decimal)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    5.00    8.00   12.50   12.26   16.00   22.83

table(puntos_hora$franja_horaria, useNA = "always")

## 
##      Madrugada    Mañana-Pico       Mediodía          Tarde     Tarde-Pico 
##           4098          43006          18450          35028          25250 
## Noche-Temprana   Noche-Tardia           <NA> 
##          19708           3034              0

# Paleta de colores para franjas
pal_hora <- colorFactor(
  palette = c("#1a2b5c", "#f7aa00", "#a4d4ae", "#ff6b6b", "#c81d25", "#5e17eb", "#2e2e2e"),
  domain = puntos_hora$franja_horaria
)

# Crear mapa
leaflet(puntos_hora) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    radius = ~ifelse(gi_score_numeric > 2, 6, 3),
    color = ~pal_hora(franja_horaria),
    fillOpacity = 0.7,
    popup = ~paste(
      "<strong>Hora:</strong> ", HoraOri, "<br>",
      "<strong>Franja:</strong> ", franja_horaria, "<br>",
      "<strong>Score:</strong> ", round(gi_score_numeric, 2)
    )
  ) %>%
  addLegend(
    position = "bottomright",
    pal = pal_hora,
    values = ~franja_horaria,
    title = "Franjas Horarias"
  ) %>%
  setView(lng = -100.3, lat = 25.67, zoom = 12)

# Resumen estadístico por franja horaria
puntos_hora %>%
  st_drop_geometry() %>%
  group_by(franja_horaria) %>%
  summarise(
    n_puntos = n(),
    avg_score = mean(gi_score_numeric),
    porcentaje = n() / nrow(puntos_hora) * 100
  ) %>%
  arrange(desc(avg_score)) %>%
  ggplot(aes(x = reorder(franja_horaria, -avg_score), y = avg_score)) +
  geom_col(aes(fill = franja_horaria)) +
  geom_text(aes(label = paste0(round(porcentaje,1), "%")), vjust = -0.5) +
  scale_fill_manual(values = c("#f7aa00", "#c81d25", "#ff6b6b", "#a4d4ae", "#1a2b5c", "#5e17eb", "#2e2e2e")) +
  labs(title = "Intensidad Promedio de Hotspots por Franja Horaria",
       x = "Franja Horaria",
       y = "Score Promedio (Gi*)") +
  theme_minimal()

Congestión significativa (hotspots positivos fuertes): - Madrugada y Noche-Temprana tienen valores de Gi* positivos → hay cierta concentración de congestión (aunque no muy intensa). - Pero su intensidad es ligera (Gi* < 1).

Congestión dispersa o sin concentración espacial (valores bajos/negativos de Gi*):

• Tarde, Mañana-Pico y Mediodía muestran valores negativos altos de Gi*, especialmente: - Tarde: -2.3 - Mañana-Pico: -1.7
• Esto indica que durante estas franjas horarias, los incidentes de congestión no se concentran espacialmente, sino que están más dispersos por la ciudad.

Transiciones interesantes:

La Noche-Tardia tiene un Gi* cercano a 0 → sin evidencia fuerte de concentración ni dispersión.

Tarde-Pico muestra apenas una ligera dispersión.

¿Qué significa esto ? En horas pico (mañana y tarde), aunque hay muchos viajes, los eventos de congestión están dispersos, lo cual sugiere que no hay zonas fijas de congestión, sino que se esparcen por varias partes de la ciudad.

En horas de baja actividad (madrugada y noche-temprana), los pocos viajes que ocurren tienden a concentrarse espacialmente, posiblemente en corredores viales clave o zonas centrales.

Planeación urbana o de transporte público debe enfocarse más en mitigar dispersión en horas pico, ya que es más difícil de controlar.

# Filtrar solo horas pico (mañana y tarde)
horas_pico <- puntos_hora %>%
  filter(franja_horaria %in% c("Mañana-Pico", "Tarde-Pico"))

# Preparar coordenadas ANTES de crear el mapa
coords <- horas_pico %>% 
  st_coordinates() %>% 
  as.data.frame()

# Crear mapa con heatmap
leaflet() %>%  # Nota: no pasamos el objeto sf directamente aquí
  addProviderTiles(providers$CartoDB.Positron) %>%
  addHeatmap(
    lng = coords$X,  # Usar las columnas ya extraídas
    lat = coords$Y,
    intensity = horas_pico$gi_score_numeric,
    radius = 15,
    blur = 20,
    max = 0.05,
    gradient = c("#f7aa00", "#c81d25")  # Gradiente de color directo
  ) %>%
  addLegend(
    position = "bottomright",
    colors = c("#f7aa00", "#c81d25"),
    labels = c("Baja intensidad", "Alta intensidad"),
    title = "Intensidad Horas Pico"
  ) %>%
  setView(lng = -100.3, lat = 25.67, zoom = 13)

Las áreas más intensamente coloreadas en rojo oscuro (alta intensidad de Gi*) corresponden a zonas con concentración estadísticamente significativa de eventos durante las horas pico. Estas incluyen:

• Zona centro de Monterrey
• San Nicolás de los Garza
• Guadalupe
• Santa Catarina
• Escobedo

También se observan focos en áreas más alejadas como Apodaca y algunas partes de Benito Juárez y Cadereyta Jiménez

Esto indica que, aunque la congestión se concentra principalmente en el corazón de la metrópoli, hay hotspots secundarios en zonas periféricas.

# Extraer coordenadas y datos
heat_data <- horas_pico %>%
  mutate(
    lon = st_coordinates(.)[,1],
    lat = st_coordinates(.)[,2]
  ) %>%
  st_drop_geometry()

# Mapa completo
leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  # Heatmap
  addHeatmap(
    data = heat_data,
    lng = ~lon,
    lat = ~lat,
    intensity = ~gi_score_numeric,
    radius = 20,
    blur = 15,
    max = 0.1,
    gradient = colorRampPalette(c("#f7aa00", "#c81d25"))(10)
  ) %>%
  
  # Puntos críticos (top 5%)
  addCircleMarkers(
    data = filter(heat_data, gi_score_numeric > quantile(gi_score_numeric, 0.95)),
    radius = ~sqrt(gi_score_numeric)*1.5,
    color = "#000000",
    fillColor = "#c81d25",
    fillOpacity = 0.8,
    popup = ~paste("Score:", round(gi_score_numeric, 2))
  ) %>%
  
  # Controles
  addLayersControl(
    overlayGroups = c("Heatmap", "Puntos Críticos")
  )

## Assuming "lon" and "lat" are longitude and latitude, respectively

Los puntos más oscuros y densos del heatmap, así como los círculos marcados, se concentran principalmente en:

• Centro de Monterrey
• San Nicolás de los Garza
• Santa Catarina
• Guadalupe
• Zonas del sur como Benito Juárez también presentan hotspots relevantes
• Algunas áreas periféricas como García o Pesquería también muestran aparición de hotspots aislados

---

Análisis Datos de Accidentes en Nuevo León

# Filtrar por EDO == 19
datos <- datos %>% 
  filter(EDO == 19) # 19 es Nuevo León

# Cambiar nombre de los municipios
datos <- datos %>%
  mutate(MPIO = recode(MPIO,
    "1" = "Abasolo",
    "6" = "Apodaca",
    "9" = "Cadereyta Jiménez",  
    "10" = "El Carmen",
    "12" = "Ciénega de Flores",
    "18" = "García",
    "19" = "San Pedro Garza García",
    "21" = "General Escobedo",
    "25" = "General Zuazua",
    "26" = "Guadalupe",
    "31" = "Juárez",
    "39" = "Monterrey",
    "41" = "Pesquería",
    "45" = "Salinas Victoria",
    "46" = "San Nicolás de los Garza",
    "47" = "Hidalgo",
    "48" = "Santa Catarina",
    "49" = "Santiago",
  ))

Hot Spots significativos de accidentes de transporte público dentro del área metropolitana de Monterrey

# Preparar los datos
datos <- datos %>% 
  filter(!is.na(LONGITUD) & !is.na(LATITUD))  # Eliminar NA

# Convertir a objeto espacial (CRS WGS84)
datos_sf <- st_as_sf(datos, coords = c("LONGITUD", "LATITUD"), crs = 4326) %>% 
  st_transform(32614)  # Proyectar a UTM zona 14N (Monterrey)

# Crear matriz de vecindad (radio de 1 km)
coords <- st_coordinates(datos_sf)
vecindad <- dnearneigh(coords, d1 = 0, d2 = 1000)  # Ajusta el radio (en metros)
pesos <- nb2listw(vecindad, style = "B", zero.policy = TRUE)

# Calcular hotspots
datos_sf$hotspot <- as.numeric(localG(datos_sf$CLASE, pesos) > 2)  # Gi* > 2 SD

# Aplicar DBSCAN (500m, mínimo 5 accidentes)
clusters <- dbscan(coords, eps = 500, minPts = 5)
datos_sf$cluster <- as.factor(clusters$cluster)

# Filtrar clústers significativos (excluir ruido, cluster = 0)
clusters_significativos <- datos_sf[datos_sf$cluster != 0, ]

# Clase como numérica
datos_sf$clase_numeric <- as.numeric(as.character(datos_sf$CLASE))

# Crear matriz de vecindad
pesos <- nb2listw(vecindad, style = "W", zero.policy = TRUE)

# Calcular Moran's I
moran.test(datos_sf$clase_numeric, pesos, zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  datos_sf$clase_numeric  
## weights: pesos  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 68.312, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.054679e-02     -1.431783e-05      2.001475e-07

# Visualización
ggplot() +
  # Hotspots (rojo = Gi* significativo)
  geom_sf(data = datos_sf, aes(color = as.factor(hotspot)), size = 2, alpha = 0.6) +
  scale_color_manual(
    values = c("0" = "gray", "1" = "red"),
    labels = c("No hotspot", "Hotspot")
  ) +
  # Clústers (círculos azules)
  geom_sf(data = clusters_significativos, aes(fill = cluster), shape = 21, size = 3) +
  scale_fill_viridis_d(option = "plasma") +
  labs(
    title = "Hotspots y Clústers de Accidentes en Monterrey",
    subtitle = "Rojo: Hotspots (Gi* > 2 SD) | Círculos: Clústers (DBSCAN)",
    caption = "Radio: 1 km (Hotspots) / 500 m (DBSCAN)"
  ) +
  theme_minimal()

ggplot() +
  geom_sf(
    data = filter(datos_sf, hotspot == 1),  # Solo hotspots
    color = "red",
    size = 2,
    alpha = 0.6
  ) +
  labs(title = "Hotspots Significativos (Gi* > 2 SD)") +
  theme_minimal()

# Resumen estadístico
cat("Número de hotspots:", sum(datos_sf$hotspot), "\n")

## Número de hotspots: NA

cat("Número de clústers:", max(clusters$cluster))

## Número de clústers: 142

# Verificar que 'CLASE' sea un factor
datos_sf$CLASE <- as.factor(datos_sf$CLASE)

# Mapa con ggplot2
ggplot() +
  geom_sf(
    data = datos_sf,
    aes(color = CLASE),  # Color por tipo de accidente
    size = 2,
    alpha = 0.7
  ) +
  scale_color_manual(
    values = c("1" = "red", "2" = "blue", "3" = "green"),  # Personaliza colores
    name = "Tipo de Accidente"  # Título de la leyenda
  ) +
  labs(
    title = "Distribución de Accidentes por Tipo en Monterrey",
    subtitle = "Color según la columna CLASE",
    caption = "Fuente: BASE MUNICIPAL 2023"
  ) +
  theme_minimal() +
  theme(
    legend.position = "bottom",
    plot.title = element_text(face = "bold", hjust = 0.5)
  )

Distribución por municipio y tipo de accidente

# Verificar conteo por municipio
tabla_mpio <- datos_sf %>%
  group_by(MPIO) %>%
  summarise(accidentes = n())

# Accidentes por tipo
tabla_clase <- datos_sf %>%
  group_by(CLASE) %>%
  summarise(cantidad = n())


ggplot(tabla_mpio, aes(x = reorder(MPIO, -accidentes), y = accidentes)) +
  geom_col(fill = "steelblue") +
  coord_flip() +
  labs(title = "Accidentes por Municipio", x = "Municipio", y = "Cantidad")

ggplot(tabla_clase, aes(x = as.factor(CLASE), y = cantidad, fill = as.factor(CLASE))) +
  geom_col() +
  labs(title = "Accidentes por Tipo (CLASE)", x = "Tipo de Accidente", y = "Cantidad") +
  theme_minimal()

Total de accidentes, muertos y heridos por municipio

resumen_mpio <- datos_sf %>%
  group_by(MPIO) %>%
  summarise(
    Total_Accidentes = n(),
    Total_Muertos = sum(TOTMUERTOS, na.rm = TRUE),
    Total_Heridos = sum(TOTHERIDOS, na.rm = TRUE)
  ) %>%
  arrange(desc(Total_Accidentes)) %>%
  slice_head(n = 10)  # Top 10 municipios


resumen_mpio_long <- resumen_mpio %>%
  pivot_longer(cols = c("Total_Muertos", "Total_Heridos"),
               names_to = "Tipo",
               values_to = "Cantidad")

ggplot(resumen_mpio_long, aes(x = reorder(MPIO, -Cantidad), y = Cantidad, fill = Tipo)) +
  geom_col(position = "dodge") +
  labs(title = "Top 10 Municipios por Severidad de Accidentes",
       x = "Municipio",
       y = "Cantidad",
       fill = "Tipo de Víctima") +
  scale_fill_manual(values = c("Total_Muertos" = "#c81d25", "Total_Heridos" = "#f7aa00")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Resumen de severidad total

datos_sf <- datos_sf %>%
  mutate(
    severidad = TOTMUERTOS + TOTHERIDOS
  )

ggplot(datos_sf, aes(x = severidad)) +
  geom_histogram(bins = 30, fill = "darkred") +
  labs(title = "Distribución de Severidad Total por Accidente", x = "Total Muertos + Heridos")

Autocorrelación Espacial (Moran’s I)

coords <- st_coordinates(datos_sf)
nb <- dnearneigh(coords, 0, 1000)  # 1 km
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

moran.test(datos_sf$severidad, lw, zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  datos_sf$severidad  
## weights: lw  
## n reduced by no-neighbour observations  
## 
## Moran I statistic standard deviate = 54.75, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      2.442388e-02     -1.431783e-05      1.992406e-07

Hotspots con Gi* en severidad

datos_sf$gi_score <- localG(datos_sf$severidad, listw = lw, zero.policy = TRUE)
datos_sf$gi_score_numeric <- as.numeric(datos_sf$gi_score)

# Mapa de Gi*
ggplot(datos_sf) +
  geom_sf(aes(color = gi_score_numeric), size = 1, alpha = 0.7) +
  scale_color_viridis(name = "Gi* Score") +
  labs(title = "Hotspots de Severidad (Gi*)", subtitle = "Valores altos indican concentración de alta severidad") +
  theme_minimal()

hotspots <- datos_sf %>%
  filter(!is.na(gi_score_numeric) & gi_score_numeric > quantile(gi_score_numeric, 0.95, na.rm = TRUE)) %>%
  st_transform(4326)

pal <- colorNumeric("Reds", domain = hotspots$gi_score_numeric)

leaflet(hotspots) %>%
  addProviderTiles(providers$CartoDB.Positron) %>%
  addCircleMarkers(
    radius = ~sqrt(gi_score_numeric)*2,
    fillColor = ~pal(gi_score_numeric),
    color = "#000000",
    fillOpacity = 0.7,
    stroke = FALSE,
    popup = ~paste("Score:", round(gi_score_numeric, 2))
  ) %>%
  addLegend(pal = pal, values = ~gi_score_numeric, title = "Gi* Score")

Distribución por hora del día

datos_sf <- datos_sf %>%
  mutate(
    hora_decimal = HORA + MINUTOS / 60,
    FRANJA = cut(
      hora_decimal,
      breaks = c(0, 6, 9, 12, 15, 18, 21, 24),
      labels = c("Madrugada", "Mañana-Pico", "Mediodía", "Tarde", "Tarde-Pico", "Noche-Temprana", "Noche-Tardía"),
      include.lowest = TRUE,
      right = FALSE
    )
  )

ggplot(datos_sf, aes(x = FRANJA)) +
  geom_bar(fill = "#0072B2") +
  labs(title = "Frecuencia de Accidentes por Franja Horaria", x = "Franja", y = "Cantidad")

Hotspots por franja horaria

# Función para mapa de hotspots por franja
mapa_hotspots_franja <- function(data, franja_label) {
  subset <- data %>% filter(FRANJA == franja_label)
  coords <- st_coordinates(subset)
  nb <- dnearneigh(coords, 0, 1000)
  lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

  subset$gi_score <- localG(subset$TOTMUERTOS + subset$TOTHERIDOS, lw, zero.policy = TRUE)

  subset_top <- subset %>% filter(gi_score > quantile(gi_score, 0.95, na.rm = TRUE))

  leaflet(subset_top %>% st_transform(4326)) %>%
    addProviderTiles(providers$CartoDB.Positron) %>%
    addCircleMarkers(
      radius = 5,
      color = "red",
      fillOpacity = 0.7,
      popup = ~paste("Score:", round(gi_score, 2))
    ) %>%
    addLegend("bottomright", colors = "red", labels = paste("Hotspots -", franja_label)) %>%
    setView(lng = -100.3, lat = 25.67, zoom = 11)
}

mapa_hotspots_franja(datos_sf, "Mañana-Pico")

mapa_hotspots_franja(datos_sf, "Mediodía")

mapa_hotspots_franja(datos_sf, "Tarde")

mapa_hotspots_franja(datos_sf, "Tarde-Pico")

mapa_hotspots_franja(datos_sf, "Noche-Tardía")

Distribución Espacial por Gravedad

datos_sf <- datos_sf %>%
  mutate(gravedad = case_when(
    TOTMUERTOS > 0 ~ "Grave (muertes)",
    TOTHERIDOS > 0 ~ "Media (heridos)",
    TRUE ~ "Leve (solo daños)"
  ))


ggplot(datos_sf) +
  geom_sf(aes(color = gravedad), alpha = 0.6) +
  scale_color_manual(values = c("Grave (muertes)" = "darkred", "Media (heridos)" = "orange", "Leve (solo daños)" = "gray")) +
  labs(title = "Distribución Espacial por Gravedad del Accidente")

Gráfico de líneas de accidentes por mes

datos_sf <- datos_sf %>%
  mutate(fecha = make_date(ANIO, MES, DIA)) %>%
  filter(!is.na(fecha))

# Agrupar por mes
mensual <- datos_sf %>%
  mutate(mes = floor_date(fecha, "month")) %>%
  group_by(mes) %>%
  summarise(n_accidentes = n())

# Gráfico
ggplot(mensual, aes(x = mes, y = n_accidentes)) +
  geom_line(color = "steelblue", size = 1.2) +
  geom_point(color = "darkred") +
  labs(
    title = "Tendencia mensual de accidentes",
    x = "Mes",
    y = "Número de accidentes"
  ) +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Pronósticos de Accidentes

library(dplyr)
library(lubridate)
library(stringr)

# Crear columna de fecha
datos <- datos %>%
  mutate(
    fecha = as.Date(paste(ANIO, MES, DIA, sep = "-")),
    diasemana_nombre = weekdays(fecha, abbreviate = FALSE),
    mes_nombre = months(fecha, abbreviate = FALSE)
  )

# Lista de columnas de transporte público
vehiculos_tp <- c("CAMPASAJ", "MICROBUS", "PASCAMION", "OMNIBUS", "TRANVIA")

# Agrupar por fecha y calcular variables
accidentes_por_fecha <- datos %>%
  group_by(fecha, diasemana_nombre, mes_nombre) %>%
  summarise(
    n_accidentes = n(),
    total_vehiculos = rowSums(across(c(AUTOMOVIL:OTROVEHIC)), na.rm = TRUE),
    vehiculos_promedio = total_vehiculos / n_accidentes,
    n_transporte_publico = rowSums(across(all_of(vehiculos_tp)), na.rm = TRUE),
    prop_transporte_publico = n_transporte_publico / total_vehiculos,
    total_muertos = sum(TOTMUERTOS, na.rm = TRUE),
    total_heridos = sum(TOTHERIDOS, na.rm = TRUE),
    tasa_mortalidad = total_muertos / n_accidentes,
    clase_accidente_frecuente = names(which.max(table(CLASE))),
    tipo_accidente_frecuente = names(which.max(table(TIPACCID))),
    .groups = "drop"
  )

## Warning: Returning more (or less) than 1 row per `summarise()` group was deprecated in
## dplyr 1.1.0.
## ℹ Please use `reframe()` instead.
## ℹ When switching from `summarise()` to `reframe()`, remember that `reframe()`
##   always returns an ungrouped data frame and adjust accordingly.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Convertir a factores para renombrar clases/tipos
accidentes_por_fecha <- accidentes_por_fecha %>%
  mutate(
    clase_accidente_frecuente = recode(clase_accidente_frecuente, `3` = "Solo daños"),
    tipo_accidente_frecuente = recode(tipo_accidente_frecuente, `1` = "Colisión con vehículo automotor")
  )

# Detectar outliers en número de accidentes
q1 <- quantile(accidentes_por_fecha$n_accidentes, 0.25)
q3 <- quantile(accidentes_por_fecha$n_accidentes, 0.75)
iqr <- q3 - q1
limite_sup <- q3 + 1.5 * iqr

accidentes_por_fecha <- accidentes_por_fecha %>%
  mutate(
    outlier = ifelse(n_accidentes > limite_sup, "Outlier", "Normal")
  )

glimpse(accidentes_por_fecha)

## Rows: 69,902
## Columns: 14
## $ fecha                     <date> 2023-01-01, 2023-01-01, 2023-01-01, 2023-01…
## $ diasemana_nombre          <chr> "domingo", "domingo", "domingo", "domingo", …
## $ mes_nombre                <chr> "enero", "enero", "enero", "enero", "enero",…
## $ n_accidentes              <int> 110, 110, 110, 110, 110, 110, 110, 110, 110,…
## $ total_vehiculos           <dbl> 1, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2, 2, 2,…
## $ vehiculos_promedio        <dbl> 0.009090909, 0.018181818, 0.018181818, 0.018…
## $ n_transporte_publico      <dbl> 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1,…
## $ prop_transporte_publico   <dbl> 0.0000000, 0.0000000, 0.5000000, 0.5000000, …
## $ total_muertos             <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ total_heridos             <dbl> 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, 23, …
## $ tasa_mortalidad           <dbl> 0.009090909, 0.009090909, 0.009090909, 0.009…
## $ clase_accidente_frecuente <chr> "Solo daños", "Solo daños", "Solo daños", "S…
## $ tipo_accidente_frecuente  <chr> "Colisión con vehículo automotor", "Colisión…
## $ outlier                   <chr> "Normal", "Normal", "Normal", "Normal", "Nor…

# 1. Agregación diaria 
accidentes_diarios <- accidentes_por_fecha %>%
  group_by(fecha) %>%
  summarise(
    n_accidentes = first(n_accidentes),
    dia_semana = first(diasemana_nombre),
    mes = first(mes_nombre),
    transporte_publico = sum(n_transporte_publico),
    muertos = sum(total_muertos),
    heridos = sum(total_heridos)
  ) %>%
  ungroup() %>%
  filter(fecha >= "2023-01-01") %>% 
  arrange(fecha)

# 2. Preparar datos para Prophet
df_prophet <- accidentes_diarios %>%
  select(ds = fecha, y = n_accidentes) %>%
  mutate(
    dia_semana = weekdays(ds),  # Extraer día de la semana como texto
    es_fin_de_semana = ifelse(dia_semana %in% c("sábado", "domingo", "Saturday", "Sunday"), 1, 0)
  )

# 3. Modelo Prophet con regresores adicionales
modelo_prophet <- prophet(
  growth = "flat",  # Asumimos estacionariedad en la tendencia
  yearly.seasonality = TRUE,
  weekly.seasonality = TRUE,
  daily.seasonality = FALSE
)

# Añadir regresores adicionales
modelo_prophet <- add_regressor(modelo_prophet, 'es_fin_de_semana')

# Ajustar modelo
modelo_prophet <- fit.prophet(modelo_prophet, df_prophet)

# 4. Crear dataframe futuro con características
futuro <- make_future_dataframe(modelo_prophet, periods = 730) %>%
  mutate(
    dia_semana = weekdays(ds),
    es_fin_de_semana = ifelse(dia_semana %in% c("sábado", "domingo", "Saturday", "Sunday"), 1, 0)
  )

# 5. Pronóstico
pronostico <- predict(modelo_prophet, futuro)

# Gráfico estático con componentes
plot(modelo_prophet, pronostico) +
  labs(title = "Pronóstico diario de accidentes 2023-2025",
       x = "Fecha", y = "Número de accidentes") +
  theme_minimal()

# Para ver los componentes estacionales
prophet_plot_components(modelo_prophet, pronostico)

Percepción del transporte

df <- read_excel("C:\\Users\\mari0\\OneDrive\\Documents\\R Studio\\8 semestre\\Evidencia\\Evidencia1_parte2\\PERCEPCION DE LA CIUDADANIA DEL TRANSPORTE PUBLICO DE LA ZONA METROPOLITANA DE MONTERREY  (Respuestas).xlsx",
                 sheet = "Respuestas de formulario 1")

colnames(df) <- c(
  "marca_temporal", "frecuencia_uso", "motivo_uso", "medio_transporte",
  "accesibilidad", "limpieza", "seguridad", "victima_delito",
  "problema_principal", "costo_percibido", "mejoras_importantes",
  "preparacion_mundial", "impacto_mundial", "edad", "genero",
  "ocupacion", "municipio_residencia", "aplicacion"
)

sum(is.na(df))

## [1] 3

moda <- function(x) {
  ux <- na.omit(unique(x))
  ux[which.max(tabulate(match(x, ux)))]
}
# Aplicar solo a columnas con NA
for (col in names(df)) {
  if (any(is.na(df[[col]]))) {
    df[[col]][is.na(df[[col]])] <- moda(df[[col]])
  }
}

Distribución por edad y género

ggplot(df, aes(x = edad)) +
  geom_bar(fill = "steelblue") +
  labs(title = "Distribución de edades", x = "Rango de edad", y = "Frecuencia") +
  theme_minimal()

ggplot(df, aes(x = genero)) +
  geom_bar(fill = "orchid") +
  labs(title = "Distribución por género", x = "Género", y = "Frecuencia") +
  theme_minimal()

Género vs Percepción de seguridad

ggplot(df, aes(x = genero, fill = seguridad)) +
  geom_bar(position = "fill") +
  labs(title = "Percepción de seguridad por género", x = "Género", y = "Proporción") +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

Percepción del costo por grupo de edad

ggplot(df, aes(x = edad, fill = costo_percibido)) +
  geom_bar(position = "fill") +
  labs(title = "Percepción del costo del transporte por edad", x = "Edad", y = "Proporción") +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

Problemas principales percibidos

ggplot(df, aes(x = problema_principal)) +
  geom_bar(fill = "darkorange") +
  coord_flip() +
  labs(title = "Problemas principales en el transporte", x = "Problema", y = "Frecuencia") +
  theme_minimal()

Pruebas Chi-Cuadrada

Prueba 1: Género vs percepción de seguridad

cat("1. Género vs Seguridad\n")

## 1. Género vs Seguridad

tabla1 <- table(df$genero, df$seguridad)
print(tabla1)

##            
##             No, me siento inseguro/a Sí, completamente
##   Femenino                        29                34
##   Masculino                       23                22
##            
##             Sí, pero con algunas preocupaciones
##   Femenino                                   46
##   Masculino                                  23

print(chisq.test(tabla1))

## 
##  Pearson's Chi-squared test
## 
## data:  tabla1
## X-squared = 1.5145, df = 2, p-value = 0.469

Prueba 2: Municipio vs preparación para el Mundial

cat("\n2. Municipio vs Preparación Mundial\n")

## 
## 2. Municipio vs Preparación Mundial

tabla2 <- table(df$municipio_residencia, df$preparacion_mundial)
print(tabla2)

##                 
##                  No, no está preparado Sí, completamente
##   Allende                            1                 0
##   Apodaca                            4                 0
##   Cadereyta                          1                 0
##   Escobedo                          17                 0
##   García                             4                 0
##   Guadalupe                         17                 0
##   Juárez                             4                 0
##   Monterrey                         72                 4
##   San Nicolás                        4                 0
##   Santa Catarina                    16                 0
##   Santiago                           7                 2
##                 
##                  Sí, pero requiere mejoras
##   Allende                                0
##   Apodaca                                2
##   Cadereyta                              0
##   Escobedo                               0
##   García                                 2
##   Guadalupe                              1
##   Juárez                                 1
##   Monterrey                             14
##   San Nicolás                            2
##   Santa Catarina                         2
##   Santiago                               0

print(chisq.test(tabla2))

## Warning in chisq.test(tabla2): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabla2
## X-squared = 24.421, df = 20, p-value = 0.2245

Prueba 3: Ocupación vs percepción del costo

cat("\n3. Ocupación vs Costo percibido\n")

## 
## 3. Ocupación vs Costo percibido

tabla3 <- table(df$ocupacion, df$costo_percibido)
print(tabla3)

##                 
##                  No, es demasiado barato No, es demasiado caro Sí, es adecuado
##   Desempleado(a)                       0                     2               0
##   Empleado(a)                          2                    77              31
##   Empresario(a)                        0                     0               1
##   Estudiante                           0                    31              33

print(chisq.test(tabla3))

## Warning in chisq.test(tabla3): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabla3
## X-squared = 13.176, df = 6, p-value = 0.04033

Prueba 4: Edad vs percepción de disponibilidad (usando problema_principal)

cat("\n4. Edad vs Problema principal (como proxy de disponibilidad)\n")

## 
## 4. Edad vs Problema principal (como proxy de disponibilidad)

tabla4 <- table(df$edad, df$problema_principal)
print(tabla4)

##                   
##                    Condiciones de las unidades Costo Puntualidad Seguridad
##   18-24 años                                17     7          48         2
##   25-34 años                                 6     5          12         1
##   35-44 años                                 3     4          19         0
##   45-54 años                                11     4          14         0
##   55 años o más                              6     0           8         2
##   Menos de 18 años                           3     0           3         2

print(chisq.test(tabla4))

## Warning in chisq.test(tabla4): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabla4
## X-squared = 28.021, df = 15, p-value = 0.02144

Prueba 5: Limpieza vs preparación para el Mundial

cat("\n5. Limpieza vs Preparación Mundial\n")

## 
## 5. Limpieza vs Preparación Mundial

tabla5 <- table(df$limpieza, df$preparacion_mundial)
print(tabla5)

##            
##             No, no está preparado Sí, completamente Sí, pero requiere mejoras
##   Bueno                        35                 4                         9
##   Excelente                     7                 0                         0
##   Mala                         61                 1                         2
##   Regular                      44                 1                        13

print(chisq.test(tabla5))

## Warning in chisq.test(tabla5): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabla5
## X-squared = 17.554, df = 6, p-value = 0.00745

Prueba 6: Puntualidad vs preparación para el Mundial

cat("\n6. Puntualidad (problema principal) vs Preparación Mundial\n")

## 
## 6. Puntualidad (problema principal) vs Preparación Mundial

tabla6 <- table(df$problema_principal, df$preparacion_mundial)
print(tabla6)

##                              
##                               No, no está preparado Sí, completamente
##   Condiciones de las unidades                    42                 0
##   Costo                                          18                 1
##   Puntualidad                                    85                 5
##   Seguridad                                       2                 0
##                              
##                               Sí, pero requiere mejoras
##   Condiciones de las unidades                         4
##   Costo                                               1
##   Puntualidad                                        14
##   Seguridad                                           5

print(chisq.test(tabla6))

## Warning in chisq.test(tabla6): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabla6
## X-squared = 24.756, df = 6, p-value = 0.0003789

Prueba 7: Disponibilidad negativa vs positiva y su relación con preparación Mundial

df$disponibilidad <- ifelse(df$problema_principal == "Disponibilidad de unidades", "Negativa", "Positiva")
cat("\n7. Disponibilidad vs Preparación Mundial\n")

## 
## 7. Disponibilidad vs Preparación Mundial

tabla7 <- table(df$disponibilidad, df$preparacion_mundial)
print(tabla7)

##           
##            No, no está preparado Sí, completamente Sí, pero requiere mejoras
##   Positiva                   147                 6                        24

print(chisq.test(tabla7))

## 
##  Chi-squared test for given probabilities
## 
## data:  tabla7
## X-squared = 199.63, df = 2, p-value < 2.2e-16

Prueba 8: Edad vs percepción del costo

cat("\n8. Edad vs Costo percibido\n")

## 
## 8. Edad vs Costo percibido

tabla8 <- table(df$edad, df$costo_percibido)
print(tabla8)

##                   
##                    No, es demasiado barato No, es demasiado caro
##   18-24 años                             2                    38
##   25-34 años                             0                    14
##   35-44 años                             0                    20
##   45-54 años                             0                    26
##   55 años o más                          0                     9
##   Menos de 18 años                       0                     3
##                   
##                    Sí, es adecuado
##   18-24 años                    34
##   25-34 años                    10
##   35-44 años                     6
##   45-54 años                     3
##   55 años o más                  7
##   Menos de 18 años               5

print(chisq.test(tabla8))

## Warning in chisq.test(tabla8): Chi-squared approximation may be incorrect

## 
##  Pearson's Chi-squared test
## 
## data:  tabla8
## X-squared = 19.893, df = 10, p-value = 0.03028