Análisis de resultados y frente de Pareto

Con los 6,000 resultados agregados por el proceso de Python y R, pasamos a la fase final de análisis.

Carga y limpieza de datos

# 1. Cargar librerías
suppressPackageStartupMessages({
  library(tidyverse)
  library(sf)
  library(skimr)
  library(GGally)
  library(plotly)
  library(kableExtra)
})

# 2. Definir rutas a los archivos con tus rutas específicas
ruta_resultados_csv <- "D:/Users/Dayling Paez/Downloads/FIA FOUNDATION/Modelo FIA Enforcement/batch_resultados_r.csv"
ruta_shapefile_vias <- "D:/Users/Dayling Paez/Downloads/FIA FOUNDATION/shp/capa_modelo_opti.shp"

# 3. Cargar y limpiar datos
resultados_limpios <- read_csv(ruta_resultados_csv, show_col_types = FALSE) %>%
  filter(!is.na(seed))

# Cargamos los datos espaciales una sola vez aquí para usarlos después.
vias_sf <- st_read(ruta_shapefile_vias, quiet = TRUE) %>%
  janitor::clean_names()

Análisis estadístico y gráficos

Resumen estadístico

resultados_limpios %>%
  select(-seed) %>%
  skim()

Data summary
Name	Piped data
Number of rows	5500
Number of columns	4
_______________________
Column type frequency:
numeric	4
________________________
Group variables	None

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
velocidad_promedio_red	1	40.43	0.03	40.32	40.41	40.43	40.45	40.55	▁▃▇▂▁
indice_exceso_red	1	0.18	0.00	0.18	0.18	0.18	0.18	0.18	▁▇▂▁▁
indice_riesgo_red	1	4.09	0.01	4.06	4.08	4.09	4.09	4.11	▁▂▇▅▁
probabilidad_siniestro_red	1	0.57	0.00	0.57	0.57	0.57	0.57	0.57	▁▃▇▅▁

Gráficos de distribución y relaciones

# Gráfico 1: Distribuciones (Violin Plots)
resultados_limpios %>%
  pivot_longer(cols = -seed, names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = variable, y = valor, fill = variable)) +
  geom_violin(trim = FALSE, alpha = 0.8) +
  geom_boxplot(width = 0.1, fill = "white") +
  facet_wrap(~variable, scales = "free") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "none", axis.text.x = element_blank()) +
  labs(title = paste("Distribución de resultados en", nrow(resultados_limpios), "Escenarios"),
       x = "Variable objetivo", y = "Valor")

# Gráfico 2: Relaciones Cruzadas (Pair Plot)
GGally::ggpairs(
  resultados_limpios,
  columns = c("velocidad_promedio_red", "indice_exceso_red", "indice_riesgo_red", "probabilidad_siniestro_red"),
  upper = list(continuous = GGally::wrap("cor", size = 4)),
  lower = list(continuous = GGally::wrap("points", alpha = 0.1, size = 0.5)),
  title = "Relaciones cruzadas entre variables objetivo"
)

Identificación de frentes de Pareto

is_dominated <- function(punto, todos_los_puntos) {
   any(
    colSums(t(todos_los_puntos) <= punto) == ncol(todos_los_puntos) &
    colSums(t(todos_los_puntos) < punto) > 0
  )
}
objetivos_min <- c("indice_exceso_red", "probabilidad_siniestro_red", "indice_riesgo_red")
valores_objetivos <- resultados_limpios %>% select(all_of(objetivos_min))
es_pareto <- !sapply(1:nrow(valores_objetivos), function(i) {
  is_dominated(as.numeric(valores_objetivos[i, ]), valores_objetivos[-i, ])
})
frente_pareto <- resultados_limpios[es_pareto, ]
cat(paste("Se encontraron", nrow(frente_pareto), "escenarios en el Frente de Pareto.\n"))

## Se encontraron 18 escenarios en el Frente de Pareto.

Visualización y Análisis del Escenario 4417

En esta sección, en lugar de calcular dinámicamente una solución de compromiso, forzamos el análisis para que se centre específicamente en el Escenario 4417 para fines de reproducibilidad y comparación.

# --- 1. DEFINIR EL PUNTO BASE (CONDICIÓN ACTUAL) ---
punto_base <- tibble(
  seed = -1,
  velocidad_promedio_red = 40.5916,
  indice_exceso_red = 0.1832,
  indice_riesgo_red = 4.0238,
  probabilidad_siniestro_red = 0.5691
)

# --- 2. FORZAR LA SELECCIÓN DE LA SEMILLA 4417 ---
# En lugar de buscar, seleccionamos directamente la fila del escenario deseado.
punto_compromiso <- resultados_limpios %>%
  filter(seed == 4417)

cat(paste("Análisis enfocado en la solución del Seed:", punto_compromiso$seed, "\n"))

## Análisis enfocado en la solución del Seed: 4417

# --- 3. PREPARAR DATOS PARA GRÁFICOS ---
datos_grafico <- bind_rows(
  resultados_limpios %>% mutate(tipo = "6000 escenarios simulados"),
  frente_pareto %>% mutate(tipo = "Frente de Pareto"),
  punto_base %>% mutate(tipo = "Situación actual"),
  punto_compromiso %>% mutate(tipo = "Solución de compromiso (Seed 4417)") # Etiqueta actualizada
) %>%
  mutate(tipo = factor(tipo, levels = c("6000 escenarios simulados", "Frente de Pareto", 
                                      "Situación actual", "Solución de compromiso (Seed 4417)")))

# --- 4. DEFINIR ESCALAS MANUALES ---
colores_manuales <- c("6000 escenarios simulados" = "grey80", "Frente de Pareto" = "#0072B2", 
                      "Situación actual" = "black", "Solución de compromiso (Seed 4417)" = "firebrick")
formas_manuales <- c("6000 escenarios simulados" = 16, "Frente de Pareto" = 16, 
                     "Situación actual" = 4, "Solución de compromiso (Seed 4417)" = 17)
tamanos_manuales <- c("6000 escenarios simulados" = 1.5, "Frente de Pareto" = 3, 
                      "Situación actual" = 5, "Solución de compromiso (Seed 4417)" = 5)

# --- 5. CREAR GRÁFICOS 2D ---
# Proyección 1: Velocidad vs Exceso
ggplot(datos_grafico, aes(x = velocidad_promedio_red, y = indice_exceso_red, color = tipo, shape = tipo, size = tipo)) +
  geom_point(alpha = 0.7, stroke = 1.2) + scale_color_manual(values = colores_manuales, name = "Leyenda") +
  scale_shape_manual(values = formas_manuales, name = "Leyenda") + scale_size_manual(values = tamanos_manuales, name = "Leyenda") +
  guides(color = guide_legend(override.aes = list(size = 4))) + labs(title = "Velocidad vs. índice de exceso", x = "Velocidad promedio (km/h)", y = "Índice de exceso") + theme_bw(base_size = 14)

# Proyección 2: Velocidad vs Siniestralidad
ggplot(datos_grafico, aes(x = velocidad_promedio_red, y = probabilidad_siniestro_red, color = tipo, shape = tipo, size = tipo)) +
  geom_point(alpha = 0.7, stroke = 1.2) + scale_color_manual(values = colores_manuales, name = "Leyenda") +
  scale_shape_manual(values = formas_manuales, name = "Leyenda") + scale_size_manual(values = tamanos_manuales, name = "Leyenda") +
  guides(color = guide_legend(override.aes = list(size = 4))) + labs(title = "Velocidad vs. probabilidad de siniestro", x = "Velocidad promedio (km/h)", y = "Probabilidad de siniestro") + theme_bw(base_size = 14)

# Proyección 3: Exceso vs Siniestralidad
ggplot(datos_grafico, aes(x = indice_exceso_red, y = probabilidad_siniestro_red, color = tipo, shape = tipo, size = tipo)) +
  geom_point(alpha = 0.7, stroke = 1.2) + scale_color_manual(values = colores_manuales, name = "Leyenda") +
  scale_shape_manual(values = formas_manuales, name = "Leyenda") + scale_size_manual(values = tamanos_manuales, name = "Leyenda") +
  guides(color = guide_legend(override.aes = list(size = 4))) + labs(title = "Análisis del Escenario 4417", x = "Índice de exceso", y = "Probabilidad de siniestro") + theme_bw(base_size = 14)

Resultados Clave del Escenario 4417

resultados_optimos <- punto_compromiso %>%
  select(
    `Semilla del Escenario` = seed,
    `Velocidad promedio (km/h)` = velocidad_promedio_red,
    `Índice de exceso` = indice_exceso_red,
    `Índice de riesgo` = indice_riesgo_red,
    `Probabilidad de siniestro` = probabilidad_siniestro_red
  )
knitr::kable(
  resultados_optimos,
  caption = "Valores de las Funciones Objetivo para el Escenario 4417",
  digits = 4
) %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE, position = "center")

Valores de las Funciones Objetivo para el Escenario 4417
Semilla del Escenario	Velocidad promedio (km/h)	Índice de exceso	Índice de riesgo	Probabilidad de siniestro
4417	40.3194	0.1769	4.0725	0.5674

Identificación de ubicaciones óptimas

obtener_ubicaciones_por_seed <- function(seed_id, datos_base_sf, n_camaras_nuevas) {
  set.seed(seed_id)
  elegibles <- datos_base_sf %>% 
    filter(tiene_cama == 0) %>%
    pull(pk_id_calz)
  n_a_tomar <- min(n_camaras_nuevas, length(elegibles))
  ubicaciones_elegidas <- sample(elegibles, size = n_a_tomar, replace = FALSE)
  return(ubicaciones_elegidas)
}

Mapa 1: Ubicaciones para el Escenario 4417

seed_optimo <- punto_compromiso$seed
ubicaciones_optimas_final <- obtener_ubicaciones_por_seed(seed_id = seed_optimo, datos_base_sf = vias_sf, n_camaras_nuevas = 100)

vias_mapa_compromiso <- vias_sf %>%
  mutate(
    tipo_camara = case_when(
      pk_id_calz %in% ubicaciones_optimas_final ~ "Nueva (solución óptima)",
      tiene_cama == 1 ~ "Existente",
      TRUE ~ "Sin cámara"
    ),
    tipo_camara = factor(tipo_camara, levels = c("Nueva (solución óptima)", "Existente", "Sin cámara"))
  )
puntos_camaras <- vias_mapa_compromiso %>%
  filter(tipo_camara %in% c("Nueva (solución óptima)", "Existente")) %>%
  st_centroid(warn = FALSE)

ggplot() +
  geom_sf(data = vias_mapa_compromiso, aes(color = tipo_camara, linewidth = tipo_camara), alpha = 0.9) +
  geom_sf(data = puntos_camaras, aes(color = tipo_camara), size = 2.5, shape = 16) +
  scale_color_manual(name = "Tipo de cámara", values = c("Nueva (solución óptima)" = "#0072B2", "Existente" = "firebrick", "Sin cámara" = "grey80")) +
  scale_linewidth_manual(values = c("Nueva (solución óptima)" = 0.8, "Existente" = 0.8, "Sin cámara" = 0.15), guide = "none") +
  labs(title = "Ubicación de cámaras para la solución de compromiso", subtitle = paste("Escenario:", seed_optimo)) +
  theme_void() +
  theme(plot.title = element_text(size = 16, face = "bold", hjust = 0.5), plot.subtitle = element_text(size = 12, hjust = 0.5), legend.position = "bottom")

Mapa 2: Mapa de Calor de Ubicaciones Prioritarias

top_seeds <- list(
  arrange(resultados_limpios, desc(velocidad_promedio_red)),
  arrange(resultados_limpios, indice_exceso_red),
  arrange(resultados_limpios, indice_riesgo_red),
  arrange(resultados_limpios, probabilidad_siniestro_red)
) %>%
  map(~ head(., 10) %>% pull(seed)) %>%
  unlist()
seeds_optimos <- unique(c(top_seeds, frente_pareto$seed))
seeds_optimos <- na.omit(seeds_optimos)

lista_completa_tramos <- unlist(lapply(seeds_optimos, function(s) {
  obtener_ubicaciones_por_seed(seed_id = s, datos_base_sf = vias_sf, n_camaras_nuevas = 100)
}))
frecuencia_tramos <- as.data.frame(table(lista_completa_tramos)) %>%
  rename(pk_id_calz = lista_completa_tramos, frecuencia = Freq) %>%
  mutate(pk_id_calz = as.numeric(as.character(pk_id_calz)))

vias_mapeo <- vias_sf %>%
  left_join(frecuencia_tramos, by = "pk_id_calz") %>%
  mutate(frecuencia = replace_na(frecuencia, 0))

ggplot() +
  geom_sf(data = vias_mapeo, color = "grey80", size = 0.2, alpha = 0.5) +
  geom_sf(data = vias_mapeo %>% filter(frecuencia > 0), aes(color = frecuencia, alpha = frecuencia, size = frecuencia)) +
  geom_sf(data = vias_mapeo %>% filter(tiene_cama == 1) %>% st_centroid(warn = FALSE), color = "firebrick", size = 1, alpha = 0.7, shape = 16) +
  scale_color_viridis_c(name = "Frecuencia de selección en escenarios óptimos") +
  scale_alpha_continuous(range = c(0.4, 1.0), guide = "none") +
  scale_size_continuous(range = c(0.4, 2), guide = "none") +
  labs(title = "Mapa de calor de ubicaciones óptimas para 100 nuevas cámaras", subtitle = "Basado en la frecuencia de aparición en los mejores escenarios de simulación", caption = "Puntos rojos representan las cámaras existentes.") +
  theme_void() +
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5), plot.caption = element_text(hjust = 0.5))

Tablas de Tramos y Corredores Prioritarios

# Tabla 1: TOP 20 TRAMOS más frecuentes
top_tramos <- vias_mapeo %>%
  st_drop_geometry() %>%
  filter(frecuencia > 0) %>%
  select(pk_id_calz, tipos_via, frecuencia) %>%
  arrange(desc(frecuencia)) %>%
  head(20)
knitr::kable(top_tramos, col.names = c("ID del tramo (pk_id_calz)", "Corredor vial (tipos_via)", "Frecuencia"),
             caption = "Top 20 tramos prioritarios para la ubicación de 100 nuevas cámaras", align = "lcr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, position = "center")

Top 20 tramos prioritarios para la ubicación de 100 nuevas cámaras
ID del tramo (pk_id_calz)	Corredor vial (tipos_via)	Frecuencia
503605	AVENIDA JORGE URIBE BOTERO ( KR 35 )	4
603810	AVENIDA LUIS CARLOS GALAN SARMIENTO ( AC 39 )	4
385342	AVENIDA FERROCARRIL DE OCCIDENTE ( CL 28 )	4
298272	AVENIDA CARACAS ( )	4
508314	AVENIDA PEPE SIERRA ( CL 116 )	4
184483	AVENIDA FUCHA ( CL 8 SUR Y CL 11 SUR )	3
91026442	AVENIDA JORGE ELIECER GAITAN ( CL 26 )	3
24122993	AVENIDA SANTA BARBARA ( KR 19 )	3
410287	AVENIDA PRIMERO DE MAYO ( )	3
523749	AVENIDA DEL CONGRESO EUCARISTICO ( AK 68 )	3
24123116	AVENIDA BOYACA ( )	3
24119890	AVENIDA CIUDAD DE QUITO ( AK 30 )	3
24120850	AVENIDA CIUDAD DE QUITO ( AK 30 )	3
91018066	AVENIDA CIUDAD DE CALI ( )	3
91019449	AVENIDA LAUREANO GOMEZ ( AK 9 )	3
504137	AVENIDA CARACAS ( )	3
24120557	AVENIDA ALFREDO D BATEMAN ( AV SUBA )	3
213434	AVENIDA DE LOS CERROS ( AV CIRCUNVALAR )	3
514327	AVENIDA CHILE ( CL 72 )	3
507140	AVENIDA CIUDAD DE LIMA ( CL 19 )	3

# Tabla 2: TOP 15 CORREDORES con mayor frecuencia acumulada
top_corredores <- vias_mapeo %>%
  st_drop_geometry() %>%
  group_by(tipos_via) %>%
  summarise(
    longitud_total_km = first(longitud_k),
    numero_tramos_totales = n(),
    frecuencia_absoluta = sum(frecuencia, na.rm = TRUE),
    numero_tramos_sugeridos = sum(frecuencia > 0)
  ) %>%
  ungroup() %>%
  mutate(
    freq_rel_tramos = frecuencia_absoluta / numero_tramos_totales,
    pct_tramos_sugeridos = numero_tramos_sugeridos / numero_tramos_totales,
    freq_por_km = if_else(longitud_total_km > 0, frecuencia_absoluta / longitud_total_km, 0),
    `% Tramos sugeridos` = scales::percent(pct_tramos_sugeridos, accuracy = 0.1)
  ) %>%
  arrange(desc(frecuencia_absoluta)) %>%
  head(15) %>%
  select( Corredor = tipos_via, `Frecuencia absoluta` = frecuencia_absoluta, `Frec. / Tramos totales` = freq_rel_tramos, `% Tramos sugeridos`, `Frec. / km` = freq_por_km)
knitr::kable(top_corredores, caption = "Top 15 corredores prioritarios", digits = 2, align = "lcrrr") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, position = "center")

Top 15 corredores prioritarios
Corredor	Frecuencia absoluta	Frec. / Tramos totales	% Tramos sugeridos	Frec. / km
AVENIDA BOYACA ( )	321	0.35	30.1%	3.23
AVENIDA CARACAS ( )	244	0.41	32.9%	3.11
AVENIDA PASEO DE LOS LIBERTADORES ( AUTONORTE )	226	0.37	32.0%	2.66
AVENIDA CIUDAD DE QUITO ( AK 30 )	191	0.42	34.9%	2.72
AVENIDA ALBERTO LLERAS CAMARGO ( AK 7 )	170	0.38	33.6%	3.48
AVENIDA JORGE ELIECER GAITAN ( CL 26 )	166	0.40	33.3%	2.20
AVENIDA CIUDAD DE CALI ( )	148	0.37	30.8%	2.81
AVENIDA CIUDAD DE VILLAVICENCIO ( )	141	0.43	35.2%	4.31
AVENIDA DE LOS CERROS ( AV CIRCUNVALAR )	132	0.41	32.6%	3.40
AVENIDA PRIMERO DE MAYO ( )	118	0.35	29.9%	3.98
AVENIDA MEDELLIN ( CL 80 )	114	0.46	39.4%	2.65
AVENIDA CHILE ( CL 72 )	102	0.45	37.0%	5.30
AVENIDA DEL CONGRESO EUCARISTICO ( AK 68 )	102	0.39	33.0%	2.06
AVENIDA DE LAS AMERICAS ( )	98	0.40	33.1%	3.06
AUTOPISTA AL LLANO ( )	90	0.38	31.4%	3.87

Distribución de Cámaras por Corredor para el Escenario 4417 )

if (!exists("ubicaciones_optimas_final")) {
  ubicaciones_optimas_final <- obtener_ubicaciones_por_seed(seed_id = seed_optimo, datos_base_sf = vias_sf, n_camaras_nuevas = 100)
}
tramos_solucion_optima <- vias_sf %>%
  filter(pk_id_calz %in% ubicaciones_optimas_final)
distribucion_por_corredor <- tramos_solucion_optima %>%
  st_drop_geometry() %>%
  group_by(Corredor = tipos_via) %>%
  summarise(`Número de Nuevas Cámaras` = n()) %>%
  arrange(desc(`Número de Nuevas Cámaras`))

knitr::kable(
  distribucion_por_corredor,
  caption = paste("Asignación de las 100 nuevas cámaras por corredor vial para el Escenario", seed_optimo,")")
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE, position = "center")

Asignación de las 100 nuevas cámaras por corredor vial para el Escenario 4417 )
Corredor	Número de Nuevas Cámaras
AVENIDA PASEO DE LOS LIBERTADORES ( AUTONORTE )	9
AVENIDA BOYACA ( )	7
AVENIDA CIUDAD DE QUITO ( AK 30 )	7
AVENIDA CARACAS ( )	6
AVENIDA ALBERTO LLERAS CAMARGO ( AK 7 )	5
AVENIDA DEL CONGRESO EUCARISTICO ( AK 68 )	5
AUTOPISTA AL LLANO ( )	4
AVENIDA DE LAS AMERICAS ( )	4
AVENIDA BATALLON CALDAS ( KR 50 Y TV 42 )	3
AVENIDA IBERIA ( CL 133 )	3
AVENIDA LAUREANO GOMEZ ( AK 9 )	3
AVENIDA AGOBERTO MEJIA CIFUENTES ( )	2
AVENIDA CAMINO DEL PRADO ( CL 138 )	2
AVENIDA CIUDAD DE VILLAVICENCIO ( )	2
AVENIDA DE LOS CERROS ( AV CIRCUNVALAR )	2
AVENIDA JORGE ELIECER GAITAN ( CL 26 )	2
AVENIDA JOSE CELESTINO MUTIS ( CL 60 )	2
AVENIDA MEDELLIN ( CL 80 )	2
AVENIDA PASEO DEL COUNTRY ( KR 15 )	2
AVENIDA PEPE SIERRA ( CL 116 )	2
AVENIDA SANTA BARBARA ( KR 19 )	2
AUTOPISTA DEL SUR Y AVENIDA FERROCARRIL ( AC 57R SUR )	1
AVENIDA ALFREDO D BATEMAN ( AV SUBA )	1
AVENIDA ALSACIA ( CL 12 )	1
AVENIDA BOSA ( )	1
AVENIDA CIUDAD DE CALI ( )	1
AVENIDA CIUDAD DE LIMA ( CL 19 )	1
AVENIDA DE LA CONSTITUCION ( KR 68D Y AV ROJAS Y TV 49 )	1
AVENIDA DE LOS COMUNEROS ( CL 6 )	1
AVENIDA DEL FERROCARRIL DEL SUR ( TV 49 )	1
AVENIDA EL RINCON ( )	1
AVENIDA FRANCISCO DE MIRANDA ( CL 45 )	1
AVENIDA FUCHA ( CL 8 SUR Y CL 11 SUR )	1
AVENIDA JIMENEZ DE QUESADA ( CL 13 )	1
AVENIDA JOSE ASUNCION SILVA ( CL 27 SUR Y CL 28 SUR )	1
AVENIDA LA ESMERALDA ( AK 60 Y AK 48 )	1
AVENIDA LA SIRENA ( CL 153 )	1
AVENIDA LUIS CARLOS GALAN SARMIENTO ( AC 39 )	1
AVENIDA MANUEL CEPEDA VARGAS ( CL 6 SUR )	1
AVENIDA MARISCAL SUCRE ( KR 19 Y KR 24 )	1
AVENIDA MONTES ( CL 3 )	1
AVENIDA MORISCA ( AC 90 )	1
AVENIDA PRIMERO DE MAYO ( )	1
AVENIDA QUIROGA ( CL 36 Y DG 39 )	1
AVENIDA SANTA LUCIA ( CL 44 SUR )	1

Modelo FIA Optimización - Análisis de Resultados

Tomás Villescas & Dayling Paez

2025-10-20