Modelo FIA Enforcement

Modelo FIA Enforcement

Este documento presenta el desarrollo de un modelo estadístico para entender los factores que influyen en la operación y seguridad vial de Bogotá. El objetivo principal de este análisis, enmarcado en el proyecto de la Fundación FIA, la Universidad de Los Andes, ITDP y el ACC, es identificar las características de las vías que se asocian con mayores velocidades de operación, índices de exceso de velocidad, índices ponderados de riesgo por siniestralidad y, fundamentalmente, con la probabilidad de que ocurra un siniestro vial. Para ello, se utiliza una base de datos detallada con información geométrica, de infraestructura, operacional y socioeconómica de 11,561 tramos viales y sus manzanas adyacentes.

Cargar librerías

library(tidyverse)
library(readxl)
library(janitor)     # limpiar nombres de variables
library(caret)
library(corrplot)    # Matriz de correlación
library(MASS)        # stepAIC (usado por caret::train con *StepAIC)
library(car)         # Diagnóstico de modelos de regresión lineal
library(pROC)        # Diagnóstico LOGIT
library(ranger)      # Random Forest
library(xgboost)     # XGBoost
library(ggplot2)
library(gridExtra)
library(pdp)         # Partial Dependence Plots para Random Forest
set.seed(123)

Aceleración de CPU para usar multiples núcleos

library(doParallel)

## Warning: package 'doParallel' was built under R version 4.4.3

## Loading required package: foreach

## 
## Attaching package: 'foreach'

## The following objects are masked from 'package:purrr':
## 
##     accumulate, when

## Loading required package: iterators

## Loading required package: parallel

cl <- makePSOCKcluster(detectCores() - 1)
registerDoParallel(cl)

Cargar datos

Con la función glimpse() se puede ver la estructura de la base de datos con 11,561 tramos viales (filas) y 45 características (columnas) para analizar. La función clean_names() de la librería janitor estandariza los nombres de las variables, dejándolos en minúscula y con guiones bajos en lugar de espacios o caracteres especiales.

# Estandarizar nombres y dejar todo en minúscula con guión bajo
vias <- read_xlsx("C:/Users/tomas/Documents/UniAndes/FIA/Modelo/Modelo FIA Enforcement/Poligonos_final_tabla.xlsx") %>% janitor::clean_names()

glimpse(vias)

## Rows: 11,561
## Columns: 45
## $ objectid      <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1…
## $ pk_id_calz    <dbl> 506785, 178752, 180590, 603999, 601252, 91020357, 159387…
## $ aashto_bue    <dbl> 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ ancho         <dbl> 11.37, 10.13, 9.81, 11.02, 6.88, 12.13, 9.12, 6.65, 5.87…
## $ camara_1      <dbl> 151.35180, 817.62667, 1368.80344, 1373.69035, 788.90999,…
## $ camara_100    <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ camara_2      <dbl> 192.16309, 949.57494, 1377.11683, 1398.42713, 794.67181,…
## $ camara_200    <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ camara_3      <dbl> 2270.3148, 1938.3025, 1438.9690, 2441.9172, 1502.5749, 7…
## $ camara_300    <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ camara_4      <dbl> 2274.8975, 2020.2705, 1512.7726, 2463.7115, 1502.5749, 1…
## $ camara_5      <dbl> 3496.7556, 2230.9255, 1554.7296, 2490.8699, 1538.1069, 2…
## $ camara_500    <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,…
## $ cantidad_c    <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2,…
## $ carriles      <dbl> 3, 3, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 3, 2, 2, 2, 2, 2, 2,…
## $ ciclorruta    <dbl> 0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1,…
## $ civ           <dbl> 3000582, 12001970, 13001270, 11010285, 2000161, 13002801…
## $ densidad_k    <dbl> 0.00000000, 0.00000000, 0.00000000, 0.00000000, 0.000000…
## $ estaashto     <chr> "ACEPTABLE", "POBRE", "GRAVE", "ACEPTABLE", "POBRE", "FA…
## $ estrato_fi    <dbl> 3, 3, 4, 5, 6, 4, 3, 3, 0, 0, 3, 0, 5, 0, 0, 0, 3, 3, 2,…
## $ flujo_1       <chr> "MEDIO", "MEDIO", "MEDIO", "ALTO", "MEDIO", "MUY ALTO", …
## $ hora_mas_5    <dbl> 0, 6, 0, 0, 0, 0, 0, 3, 0, 0, 0, 0, 21, 6, 0, 6, 8, 0, 0…
## $ longitud_k    <dbl> 6.713203, 8.689926, 8.689926, 7.090326, 6.254239, 21.069…
## $ longitudtr    <dbl> 106.65589, 60.96325, 63.94370, 132.94178, 10.03453, 31.6…
## $ mean_denem    <dbl> 670.62158, 243.02288, 110.80748, 15.43021, 109.59107, 21…
## $ mean_denpo    <dbl> 166.5713824, 446.3473554, 274.8294960, 55.0010019, 101.2…
## $ mean_estra    <dbl> 2.66694441, 2.94095363, 3.04038592, 4.02889950, 5.408990…
## $ mean_v_ref    <dbl> 2326059.8, 1623369.1, 1546083.4, 909783.2, 4176573.8, 19…
## $ n_arboles     <dbl> 4, 10, 0, 3, 5, 3, 4, 3, 47, 3, 0, 8, 0, 1, 0, 1, 0, 15,…
## $ n_muer        <chr> "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "…
## $ n_semaforo    <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,…
## $ object        <dbl> 2, 3, 5, 6, 7, 8, 9, 10, 13, 14, 15, 19, 20, 21, 22, 23,…
## $ pres_anden    <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ rutas_sitp    <dbl> 46, 4, 15, 19, 11, 11, 4, 26, 13, 25, 23, 4, 10, 1, 21, …
## $ shape_leng    <dbl> 237.99538, 143.19366, 147.04768, 287.90036, 33.38148, 85…
## $ sitp          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, 1,…
## $ tipo          <chr> "Zona Histórica", "Residential", "Residential", "Residen…
## $ tipomallae    <chr> "Arterial", "Arterial", "Arterial", "Arterial", "Arteria…
## $ tipos_via     <chr> "AVENIDA CIUDAD DE LIMA ( CL 19 )", "AVENIDA COLOMBIA ( …
## $ vel_promed    <dbl> 18.28583, 45.75000, 29.54750, 14.90000, 28.92042, 36.205…
## $ zat           <dbl> 993, 234, 287, 133, 260, 292, 303, 187, 381, 230, 125, 8…
## $ shape_length  <dbl> 237.99538, 143.19366, 147.04768, 287.90036, 33.38148, 85…
## $ shape_area    <dbl> -1209.73507, -606.46401, -627.05544, -1454.93216, -66.84…
## $ n_her         <dbl> 7, 0, 2, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 3, 0, 3, 1, 2,…
## $ indice_riesgo <dbl> 7, 0, 2, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 2, 3, 0, 3, 1, 2,…

Variables

Se crean nuevas variables derivadas y se convierten algunas variables categóricas en factores con niveles de referencia específicos. Además, se asegura que la variable dependiente de siniestros_dummy para el modelo logit sea un factor binario que marca la existencia de un siniestro con heridos o muertos. De igual manera, la variable dependiente indice_exceso se define como la proporción de horas en las que la velocidad promedio excede los 50 km/h, limitada entre 0 y 1.

vias <- vias %>%
  mutate(
    # Derivadas
    n_muer = as.integer(n_muer),
    n_her = as.integer(n_her),
    indice_exceso = pmin(pmax(hora_mas_5 / 24, 0), 1),
    siniestros_dummy_num = as.integer(replace_na(n_muer,0) + replace_na(n_her,0) > 0),
    
    # Ingeniería de cámaras (evita multicolinealidad)
    dist_min_camara = suppressWarnings(pmin(camara_1, camara_2, camara_3, camara_4, camara_5, na.rm = TRUE)),
    
    # Eliminar "SIN DATO" de estaashto
    estaashto = if_else(estaashto == "SIN DATO", "BUENO", estaashto),
    
    # Factores + definición de niveles de referencia
    tipomallae = forcats::fct_relevel(as.factor(tipomallae), "Arterial"),
    estaashto = forcats::fct_relevel(as.factor(estaashto), "FALLADO"),
    aashto_bue = forcats::fct_relevel(as.factor(aashto_bue), "0"),
    tipo = forcats::fct_relevel(as.factor(tipo), "Sin Construir"),
    estrato_fi = forcats::fct_relevel(as.factor(estrato_fi), "1"),
    flujo_1 = forcats::fct_relevel(as.factor(flujo_1), "BAJO"),
    pres_anden = forcats::fct_relevel(as.factor(pres_anden), "0"),
    ciclorruta = forcats::fct_relevel(as.factor(ciclorruta), "0"),
    sitp = forcats::fct_relevel(as.factor(sitp), "0"),
    
    # Variable objetivo para el Logit: la creo como factor con niveles "Sí"/"No"
    siniestros_dummy = factor(if_else(siniestros_dummy_num == 1, "Sí", "No"), levels = c("Sí", "No"))
  )

vias %>% dplyr::select(vel_promed, indice_exceso, indice_riesgo, siniestros_dummy) %>% summary()

##    vel_promed    indice_exceso   indice_riesgo    siniestros_dummy
##  Min.   : 2.25   Min.   :0.000   Min.   :  0.00   Sí:6361         
##  1st Qu.:27.66   1st Qu.:0.000   1st Qu.:  0.00   No:5200         
##  Median :36.14   Median :0.000   Median :  1.00                   
##  Mean   :36.01   Mean   :0.152   Mean   :  5.84                   
##  3rd Qu.:43.76   3rd Qu.:0.250   3rd Qu.:  4.00                   
##  Max.   :80.78   Max.   :1.000   Max.   :146.00

Test de colinealidad

Antes de construir los modelos, realizo un test de multicolinealidad. El gráfico de calor que se genera aquí me permite detectar visualmente variables que están fuertemente correlacionadas (que miden lo mismo).

Los cuadrados de color rojo intenso indican una fuerte correlación positiva (ej. ancho y carriles), mientras que los azules intensos indican una fuerte correlación negativa. Los colores pálidos o blancos señalan una relación débil.

num_vars <- vias %>%
  dplyr::select(where(is.numeric)) %>%
  dplyr::select(-objectid, -pk_id_calz, -civ, -shape_area, -shape_leng, -shape_length)

cor_mat <- cor(num_vars, use = "pairwise.complete.obs")
corrplot(cor_mat, type = "upper", method = "color", tl.cex = .6, order = "hclust")

nzv_tbl <- caret::nearZeroVar(num_vars, saveMetrics = TRUE) %>%
  tibble::rownames_to_column("var") %>%
  filter(nzv) %>%
  arrange(freqRatio)
nzv_tbl

##          var freqRatio percentUnique zeroVar  nzv
## 1 camara_100  24.29759    0.01729954   FALSE TRUE

Limpieza de variables con alta colinealidad

Tras revisar la matriz de correlación y el análisis de varianza, selecciono una serie de variabels predictoras en el objeto x_comunes. Así, se plantean las formulas de los modelos lineales y del logit. Además, en el logit se incluye vel_promed e indice_exceso como predictor, ya que puede influir en la probabilidad de siniestros.

# Núcleo común de predictores
x_comunes <- ~ 
  # Variables de Jerarquía y Diseño Vial
  tipomallae + ancho + carriles + ciclorruta + 
  # Variables Socioeconómicas y de Entorno
  mean_v_ref + mean_denpo + mean_denem + n_arboles + 
  # Variables Operacionales
  n_semaforo + rutas_sitp + estaashto +
  # Variables de Cámaras y Corredor
  dist_min_camara + longitud_k

# 3 modelos lineales
form_vel   <- as.formula(update(vel_promed    ~ 1, x_comunes))
form_exc   <- as.formula(update(indice_exceso ~ 1, x_comunes))
form_risk  <- as.formula(update(indice_riesgo ~ 1, update(x_comunes, . ~ . + vel_promed + indice_exceso)))

# Modelo lineal para medir la eficiencia de cada rango de distancia de la cámara (camara_1, camara_2, camara_3, camara_4 y camara_5) y eliminando dist_min_camara
form_vel_buffer <- as.formula(update(vel_promed ~ 1, update(x_comunes, . ~ . - dist_min_camara + camara_100 + camara_200 + camara_300 + camara_500)))

# Logit (presencia de siniestros) con las variables vel_promed e indice_exceso añadidas

form_logit <- as.formula(update(siniestros_dummy ~ 1, update(x_comunes, . ~ . + vel_promed + indice_exceso)))

Modelos de regresión lineal

Para construir los modelos de regresión lineal y logit, primero divido los datos en conjuntos de entrenamiento (70%) y prueba (30%). Luego, la función train() de la librería caret permite ajustar los modelos con validación cruzada de 5 pliegues y preprocesamiento que incluye centrado y escalado de las variables predictoras.

Con el método "lmStepAIC" se realiza una regresión por pasos (stepwise) que selecciona automáticamente el mejor subconjunto de variables predictoras de la lista x_comunes. El algoritmo busca la combinación de variables que minimice el Criterio de Información de Akaike (AIC), un indicador que balancea la precisión del modelo con su simplicidad. Se realiza esta regresión por pasos para cada uno de los tres modelos lineales: velocidad promedio (vel_promed), índice de exceso de velocidad (indice_exceso) e índice ponderado de riesgo por siniestralidad (indice_riesgo).

Adicionalmente, se construye un modelo de regresión lineal para evaluar la eficiencia de diferentes rangos de distancia a cámaras de fotodetección (camara_100, camara_200, camara_300, camara_500), excluyendo la variable dist_min_camara para evitar multicolinealidad.

Resultados stepwise

En este bloque, presento los resultados finales de los tres modelos lineales seleccionados por el algoritmo stepAIC.

La función stepAIC() se aplica a cada uno de los modelos lineales ajustados previamente, utilizando los datos de entrenamiento correspondientes. El argumento direction = "both" indica que el algoritmo puede agregar o eliminar variables en cada paso para encontrar el modelo con el AIC más bajo. El parámetro trace = FALSE suprime la salida detallada del proceso de selección.

La función summary() muestra los coeficientes estimados, errores estándar, valores t y significancia estadística de cada variable incluida en el modelo final.

audit_vel   <- MASS::stepAIC(lm(form_vel,   data = train_vel),  direction = "both", trace = FALSE)
audit_exc   <- MASS::stepAIC(lm(form_exc,   data = train_exc),  direction = "both", trace = FALSE)
audit_risk  <- MASS::stepAIC(lm(form_risk,  data = train_risk), direction = "both", trace = FALSE)
audit_vel_buffer <- MASS::stepAIC(lm(form_vel_buffer, data = train_vel_buffer), direction = "both", trace = FALSE)

summary(audit_vel)

## 
## Call:
## lm(formula = vel_promed ~ tipomallae + ancho + carriles + ciclorruta + 
##     mean_v_ref + mean_denpo + mean_denem + n_arboles + n_semaforo + 
##     rutas_sitp + estaashto + dist_min_camara + longitud_k, data = train_vel)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.304  -6.617   0.308   6.755  43.189 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             2.295e+01  2.434e+00   9.427  < 2e-16 ***
## tipomallaeTroncal       2.467e+00  3.210e-01   7.687 1.68e-14 ***
## ancho                   4.596e-01  6.459e-02   7.116 1.20e-12 ***
## carriles                1.330e+00  2.332e-01   5.704 1.21e-08 ***
## ciclorruta1             8.474e-01  2.311e-01   3.666 0.000248 ***
## mean_v_ref              1.202e-06  1.500e-07   8.013 1.27e-15 ***
## mean_denpo             -7.109e-03  6.677e-04 -10.648  < 2e-16 ***
## mean_denem             -2.075e-02  1.173e-03 -17.700  < 2e-16 ***
## n_arboles               7.398e-02  1.853e-02   3.992 6.60e-05 ***
## n_semaforo             -3.894e+00  3.469e-01 -11.226  < 2e-16 ***
## rutas_sitp             -2.755e-02  8.030e-03  -3.431 0.000604 ***
## estaashtoACEPTABLE      2.197e+00  2.401e+00   0.915 0.360208    
## estaashtoBUENO          5.195e+00  2.385e+00   2.178 0.029443 *  
## estaashtoGRAVE          3.466e+00  2.574e+00   1.347 0.178175    
## estaashtoMUY POBRE      2.007e+00  2.446e+00   0.821 0.411908    
## estaashtoPOBRE          2.279e+00  2.421e+00   0.941 0.346533    
## estaashtoSATISFACTORIO  2.805e+00  2.389e+00   1.174 0.240485    
## dist_min_camara         5.625e-04  7.854e-05   7.162 8.64e-13 ***
## longitud_k              8.129e-02  4.301e-03  18.900  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.08 on 8074 degrees of freedom
## Multiple R-squared:  0.2075, Adjusted R-squared:  0.2057 
## F-statistic: 117.4 on 18 and 8074 DF,  p-value: < 2.2e-16

summary(audit_exc)

## 
## Call:
## lm(formula = indice_exceso ~ tipomallae + ancho + carriles + 
##     ciclorruta + mean_v_ref + mean_denpo + mean_denem + n_arboles + 
##     n_semaforo + estaashto + dist_min_camara + longitud_k, data = train_exc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.53038 -0.13618 -0.05133  0.06631  0.99759 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)            -1.643e-01  5.354e-02  -3.068  0.00216 ** 
## tipomallaeTroncal       9.230e-02  6.591e-03  14.004  < 2e-16 ***
## ancho                   1.130e-02  1.421e-03   7.953 2.07e-15 ***
## carriles                2.984e-02  5.128e-03   5.818 6.19e-09 ***
## ciclorruta1             3.335e-02  5.081e-03   6.564 5.56e-11 ***
## mean_v_ref              3.318e-08  3.287e-09  10.095  < 2e-16 ***
## mean_denpo             -2.094e-04  1.470e-05 -14.248  < 2e-16 ***
## mean_denem             -3.470e-04  2.572e-05 -13.493  < 2e-16 ***
## n_arboles               6.327e-04  4.075e-04   1.553  0.12048    
## n_semaforo             -5.545e-02  7.635e-03  -7.263 4.14e-13 ***
## estaashtoACEPTABLE      3.969e-02  5.285e-02   0.751  0.45266    
## estaashtoBUENO          1.028e-01  5.250e-02   1.959  0.05019 .  
## estaashtoGRAVE          6.173e-02  5.666e-02   1.090  0.27596    
## estaashtoMUY POBRE      3.329e-02  5.383e-02   0.618  0.53635    
## estaashtoPOBRE          3.333e-02  5.328e-02   0.626  0.53158    
## estaashtoSATISFACTORIO  4.838e-02  5.259e-02   0.920  0.35762    
## dist_min_camara         1.797e-05  1.694e-06  10.606  < 2e-16 ***
## longitud_k              1.452e-03  8.905e-05  16.306  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2219 on 8075 degrees of freedom
## Multiple R-squared:  0.2163, Adjusted R-squared:  0.2146 
## F-statistic: 131.1 on 17 and 8075 DF,  p-value: < 2.2e-16

summary(audit_risk)

## 
## Call:
## lm(formula = indice_riesgo ~ tipomallae + carriles + ciclorruta + 
##     mean_v_ref + mean_denpo + mean_denem + n_arboles + n_semaforo + 
##     rutas_sitp + estaashto + dist_min_camara + longitud_k, data = train_risk)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.528  -6.093  -3.411   0.042 139.928 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             6.485e-01  3.152e+00   0.206  0.83698    
## tipomallaeTroncal       2.424e+00  4.166e-01   5.818 6.19e-09 ***
## carriles                8.961e-01  2.234e-01   4.011 6.11e-05 ***
## ciclorruta1             2.113e+00  2.983e-01   7.085 1.50e-12 ***
## mean_v_ref             -1.952e-06  1.948e-07 -10.018  < 2e-16 ***
## mean_denpo              1.892e-03  8.658e-04   2.185  0.02892 *  
## mean_denem              4.739e-03  1.521e-03   3.116  0.00184 ** 
## n_arboles              -5.742e-02  2.393e-02  -2.399  0.01646 *  
## n_semaforo              2.466e+00  4.504e-01   5.475 4.52e-08 ***
## rutas_sitp              1.058e-01  1.042e-02  10.156  < 2e-16 ***
## estaashtoACEPTABLE      5.882e-01  3.118e+00   0.189  0.85037    
## estaashtoBUENO          2.328e+00  3.097e+00   0.752  0.45216    
## estaashtoGRAVE          5.378e-01  3.342e+00   0.161  0.87216    
## estaashtoMUY POBRE      1.096e-01  3.176e+00   0.035  0.97248    
## estaashtoPOBRE          8.324e-01  3.143e+00   0.265  0.79112    
## estaashtoSATISFACTORIO  1.082e+00  3.102e+00   0.349  0.72721    
## dist_min_camara        -1.577e-04  1.020e-04  -1.546  0.12209    
## longitud_k              7.966e-03  5.575e-03   1.429  0.15310    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.09 on 8075 degrees of freedom
## Multiple R-squared:  0.05729,    Adjusted R-squared:  0.05531 
## F-statistic: 28.87 on 17 and 8075 DF,  p-value: < 2.2e-16

summary(audit_vel_buffer)

## 
## Call:
## lm(formula = vel_promed ~ tipomallae + ancho + carriles + ciclorruta + 
##     mean_v_ref + mean_denpo + mean_denem + n_arboles + n_semaforo + 
##     rutas_sitp + estaashto + longitud_k + camara_100 + camara_200, 
##     data = train_vel_buffer)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.011  -6.678   0.217   6.782  43.723 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             2.425e+01  2.434e+00   9.963  < 2e-16 ***
## tipomallaeTroncal       1.999e+00  3.152e-01   6.341 2.40e-10 ***
## ancho                   4.638e-01  6.490e-02   7.147 9.64e-13 ***
## carriles                1.388e+00  2.340e-01   5.932 3.12e-09 ***
## ciclorruta1             6.256e-01  2.297e-01   2.723  0.00648 ** 
## mean_v_ref              8.999e-07  1.441e-07   6.244 4.47e-10 ***
## mean_denpo             -7.731e-03  6.642e-04 -11.640  < 2e-16 ***
## mean_denem             -2.064e-02  1.176e-03 -17.545  < 2e-16 ***
## n_arboles               7.120e-02  1.860e-02   3.828  0.00013 ***
## n_semaforo             -3.927e+00  3.482e-01 -11.278  < 2e-16 ***
## rutas_sitp             -3.849e-02  7.904e-03  -4.870 1.14e-06 ***
## estaashtoACEPTABLE      2.456e+00  2.408e+00   1.020  0.30793    
## estaashtoBUENO          5.505e+00  2.392e+00   2.301  0.02141 *  
## estaashtoGRAVE          3.673e+00  2.582e+00   1.423  0.15484    
## estaashtoMUY POBRE      2.297e+00  2.453e+00   0.936  0.34913    
## estaashtoPOBRE          2.533e+00  2.428e+00   1.043  0.29683    
## estaashtoSATISFACTORIO  3.177e+00  2.396e+00   1.326  0.18486    
## longitud_k              8.416e-02  4.332e-03  19.427  < 2e-16 ***
## camara_100              1.211e+00  8.364e-01   1.448  0.14764    
## camara_200             -1.273e+00  6.313e-01  -2.017  0.04375 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.11 on 8073 degrees of freedom
## Multiple R-squared:  0.2029, Adjusted R-squared:  0.201 
## F-statistic: 108.1 on 19 and 8073 DF,  p-value: < 2.2e-16

LOGIT

Finalmente, construyo un modelo de regresión logística (logit) para predecir la probabilidad de que ocurra un siniestro vial con heridos o muertos en un tramo vial.

# Ver distribución de la variable dependiente
table(df_train$siniestros_dummy)

## 
##   Sí   No 
## 4453 3640

# Logit simple
modelo_logit <- glm(form_logit,
  data = df_train, family = binomial(link = "logit")
)

summary(modelo_logit)

## 
## Call:
## glm(formula = form_logit, family = binomial(link = "logit"), 
##     data = df_train)
## 
## Coefficients:
##                          Estimate Std. Error z value Pr(>|z|)    
## (Intercept)             9.048e-01  5.097e-01   1.775 0.075841 .  
## tipomallaeTroncal      -3.819e-01  6.783e-02  -5.630 1.80e-08 ***
## ancho                  -4.687e-02  1.360e-02  -3.448 0.000565 ***
## carriles               -1.053e-01  4.859e-02  -2.168 0.030175 *  
## ciclorruta1            -2.088e-01  4.795e-02  -4.355 1.33e-05 ***
## mean_v_ref              3.620e-07  3.198e-08  11.319  < 2e-16 ***
## mean_denpo             -5.303e-04  1.401e-04  -3.785 0.000154 ***
## mean_denem             -1.800e-03  2.561e-04  -7.027 2.10e-12 ***
## n_arboles               1.446e-02  3.961e-03   3.650 0.000262 ***
## n_semaforo             -5.325e-01  7.593e-02  -7.013 2.33e-12 ***
## rutas_sitp             -1.970e-02  1.733e-03 -11.362  < 2e-16 ***
## estaashtoACEPTABLE     -4.606e-01  4.937e-01  -0.933 0.350851    
## estaashtoBUENO         -5.363e-01  4.905e-01  -1.093 0.274287    
## estaashtoGRAVE         -6.986e-01  5.296e-01  -1.319 0.187079    
## estaashtoMUY POBRE     -5.644e-01  5.029e-01  -1.122 0.261739    
## estaashtoPOBRE         -4.573e-01  4.977e-01  -0.919 0.358134    
## estaashtoSATISFACTORIO -4.681e-01  4.913e-01  -0.953 0.340702    
## dist_min_camara         5.404e-05  1.611e-05   3.354 0.000795 ***
## longitud_k              9.619e-04  9.130e-04   1.054 0.292065    
## vel_promed              6.236e-03  3.322e-03   1.877 0.060482 .  
## indice_exceso          -9.753e-02  1.500e-01  -0.650 0.515539    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 11137  on 8092  degrees of freedom
## Residual deviance: 10554  on 8072  degrees of freedom
## AIC: 10596
## 
## Number of Fisher Scoring iterations: 4

Modelo Random Forest

El modelo Random Forest es un método de aprendizaje de máquina basado en árboles de decisión. Es útil cuando hay relaciones no lineales en los datos y cuando se busca una alta precisión predictiva. El modelo Random Forest es un algoritmo de aprendizaje supervisado que se basa en la construcción de múltiples árboles de decisión y la combinación de sus predicciones para mejorar la precisión y reducir el sobreajuste.

En este bloque, entreno cuatro modelos de Random Forest, uno para cada una de mis variables objetivo, utilizando la implementación ranger que es significativamente más rápida. Para evitar el sobreajuste, en lugar de dejar que el modelo crezca sin control, defino una grilla de parámetros (tuneGrid_rf) que fuerza a los árboles a ser más simples, probando con un número limitado de variables en cada división (mtry) y exigiendo un número mínimo de observaciones en cada nodo final (min.node.size).

# Random Forest
set.seed(123)
ctrl_clf$allowParallel <- TRUE # Para que sea más rápido
ctrl_reg$allowParallel <- TRUE

# Define una grilla de parámetros para evitar sobreajuste
tuneGrid_rf <- expand.grid(
  .mtry = c(2, 5, 10),  # Prueba con menos variables en cada árbol
  .splitrule = "variance",
  .min.node.size = c(10, 20) # Fuerza a que cada "hoja" del árbol tenga al menos 10 o 20 observaciones
)

# Se construyen los 4 modelos de Random Forest
rf_vel <- train(
  form_vel, data = train_vel,
  method = "ranger",
  trControl = ctrl_reg,
  preProcess = c("center","scale"),
  tuneGrid = tuneGrid_rf,
  metric = "RMSE",
  importance = 'permutation'
)

rf_exc <- train(
  form_exc, data = train_exc,
  method = "ranger",
  trControl = ctrl_reg,
  preProcess = c("center","scale"),
  tuneGrid = tuneGrid_rf,
  metric = "RMSE",
  importance = 'permutation'
)

rf_risk <- train(
  form_risk, data = train_risk,
  method = "ranger",
  trControl = ctrl_reg,
  preProcess = c("center","scale"),
  tuneGrid = tuneGrid_rf,
  metric = "RMSE",
  importance = 'permutation'
)

rf_log <- train(
  form_logit, data = train_log,
  method = "ranger",
  trControl = ctrl_clf,
  preProcess = c("center","scale"),
  tuneLength = 5,
  metric = "ROC",
  importance = 'permutation'
)

rf_vel

## Random Forest 
## 
## 8093 samples
##   13 predictor
## 
## Pre-processing: centered (18), scaled (18) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6474, 6473, 6474, 6476, 6475 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  RMSE      Rsquared   MAE     
##    2    10             8.872265  0.4321157  6.906363
##    2    20             8.998369  0.4159980  7.032423
##    5    10             8.073841  0.4994363  5.964992
##    5    20             8.237318  0.4829361  6.180401
##   10    10             7.964102  0.5083760  5.797593
##   10    20             8.086914  0.4953538  5.974505
## 
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 10, splitrule = variance
##  and min.node.size = 10.

plot(varImp(rf_vel), top= 20)

rf_exc

## Random Forest 
## 
## 8093 samples
##   13 predictor
## 
## Pre-processing: centered (18), scaled (18) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6474, 6475, 6475, 6473, 6475 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  RMSE       Rsquared   MAE      
##    2    10             0.1898146  0.4618782  0.1314205
##    2    20             0.1924523  0.4466231  0.1337566
##    5    10             0.1753296  0.5159030  0.1126765
##    5    20             0.1787736  0.4997706  0.1170068
##   10    10             0.1746032  0.5170763  0.1099883
##   10    20             0.1771262  0.5046334  0.1137163
## 
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 10, splitrule = variance
##  and min.node.size = 10.

plot(varImp(rf_exc), top= 20)

rf_risk

## Random Forest 
## 
## 8093 samples
##   15 predictor
## 
## Pre-processing: centered (20), scaled (20) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6476, 6475, 6474, 6473, 6474 
## Resampling results across tuning parameters:
## 
##   mtry  min.node.size  RMSE      Rsquared    MAE     
##    2    10             12.87282  0.08559415  7.113752
##    2    20             12.86618  0.08783738  7.105279
##    5    10             12.95366  0.08337872  7.189869
##    5    20             12.88523  0.08719986  7.150558
##   10    10             13.03338  0.08233801  7.262006
##   10    20             12.97931  0.08246554  7.217442
## 
## Tuning parameter 'splitrule' was held constant at a value of variance
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were mtry = 2, splitrule = variance
##  and min.node.size = 20.

plot(varImp(rf_risk), top= 20)

rf_log

## Random Forest 
## 
## 8093 samples
##   15 predictor
##    2 classes: 'Sí', 'No' 
## 
## Pre-processing: centered (20), scaled (20) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6475, 6474, 6475, 6474, 6474 
## Resampling results across tuning parameters:
## 
##   mtry  splitrule   ROC        Sens       Spec     
##    2    gini        0.6765952  0.7421947  0.5093407
##    2    extratrees  0.6718378  0.7792517  0.4500000
##    6    gini        0.6632170  0.7040192  0.5500000
##    6    extratrees  0.6583215  0.6939154  0.5354396
##   11    gini        0.6606212  0.6993019  0.5505495
##   11    extratrees  0.6576179  0.6909994  0.5392857
##   15    gini        0.6593593  0.6984017  0.5472527
##   15    extratrees  0.6576364  0.6916710  0.5425824
##   20    gini        0.6580655  0.6977286  0.5475275
##   20    extratrees  0.6582884  0.6880795  0.5417582
## 
## Tuning parameter 'min.node.size' was held constant at a value of 1
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 2, splitrule = gini
##  and min.node.size = 1.

plot(varImp(rf_log), top= 20)

Modelo XGBoost

Para llevar el análisis predictivo al siguiente nivel, implemento un modelo de eXtreme Gradient Boosting (XGBoost). Este es considerado actualmente uno de los algoritmos de mejor rendimiento para datos tabulares. A diferencia del Random Forest que construye árboles de forma independiente, XGBoost los construye de manera secuencial, donde cada nuevo árbol se enfoca en corregir los errores cometidos por el anterior.

Nuevamente, para controlar el riesgo de sobreajuste, defino una grilla de hiperparámetros (tuneGrid_xgb) que restringe la complejidad del modelo, limitando la profundidad de los árboles (max_depth) e introduciendo penalizaciones (gamma). Entrenaré un modelo XGBoost para cada una de mis cuatro variables objetivo para comparar su rendimiento con los modelos anteriores.

# Define una grilla con parámetros más restrictivos para evitar sobreajuste
tuneGrid_xgb <- expand.grid(
  nrounds = 100,
  eta = 0.05,
  max_depth = c(3, 4), # Profundidad máxima del árbol más pequeña
  gamma = 1, # Mayor regularización
  colsample_bytree = 0.7,
  min_child_weight = 5, # Nodos más grandes
  subsample = 0.7
)

# Modelo XGBoost para vel_promed
xgb_vel <- train(
  form_vel, data = train_vel,
  method = "xgbTree", # <-- El método para XGBoost
  trControl = ctrl_reg,
  preProcess = c("center", "scale"),
  tuneGrid = tuneGrid_xgb,
  metric = "RMSE"
)

# Modelo XGBoost para indice_exceso
xgb_exc <- train(
  form_exc, data = train_exc,
  method = "xgbTree", # <-- El método para XGBoost
  trControl = ctrl_reg,
  preProcess = c("center", "scale"),
  tuneGrid = tuneGrid_xgb,
  metric = "RMSE"
)

# Modelo XGBoost para indice_riesgo
xgb_risk <- train(
  form_risk, data = train_risk,
  method = "xgbTree", # <-- El método para XGBoost
  trControl = ctrl_reg,
  preProcess = c("center", "scale"),
  tuneGrid = tuneGrid_xgb,
  metric = "RMSE"
)

# Modelo XGBoost para el logit (clasificación)
xgb_logit <- train(
  form_logit, data = train_log,
  method = "xgbTree", # <-- El método para XGBoost
  trControl = ctrl_clf,
  preProcess = c("center", "scale"),
  tuneGrid = expand.grid(nrounds = c(100, 200), # Un ejemplo de ajuste de hiperparámetros
                           eta = c(0.05, 0.1),
                           max_depth = c(3, 5),
                           gamma = 0,
                           colsample_bytree = 0.8,
                           min_child_weight = 1,
                           subsample = 0.8),
  metric = "ROC"
)

xgb_vel

## eXtreme Gradient Boosting 
## 
## 8093 samples
##   13 predictor
## 
## Pre-processing: centered (18), scaled (18) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6475, 6474, 6475, 6473, 6475 
## Resampling results across tuning parameters:
## 
##   max_depth  RMSE      Rsquared   MAE     
##   3          9.417173  0.3213195  7.423847
##   4          9.147025  0.3584736  7.166943
## 
## Tuning parameter 'nrounds' was held constant at a value of 100
## Tuning
##  parameter 'min_child_weight' was held constant at a value of 5
## 
## Tuning parameter 'subsample' was held constant at a value of 0.7
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nrounds = 100, max_depth = 4, eta
##  = 0.05, gamma = 1, colsample_bytree = 0.7, min_child_weight = 5 and
##  subsample = 0.7.

plot(varImp(xgb_vel), top= 20)

xgb_exc

## eXtreme Gradient Boosting 
## 
## 8093 samples
##   13 predictor
## 
## Pre-processing: centered (18), scaled (18) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6475, 6474, 6475, 6475, 6473 
## Resampling results across tuning parameters:
## 
##   max_depth  RMSE       Rsquared   MAE      
##   3          0.2036976  0.3558674  0.1429910
##   4          0.2005517  0.3764835  0.1402436
## 
## Tuning parameter 'nrounds' was held constant at a value of 100
## Tuning
##  parameter 'min_child_weight' was held constant at a value of 5
## 
## Tuning parameter 'subsample' was held constant at a value of 0.7
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nrounds = 100, max_depth = 4, eta
##  = 0.05, gamma = 1, colsample_bytree = 0.7, min_child_weight = 5 and
##  subsample = 0.7.

plot(varImp(xgb_exc), top= 20)

xgb_risk

## eXtreme Gradient Boosting 
## 
## 8093 samples
##   15 predictor
## 
## Pre-processing: centered (20), scaled (20) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6473, 6474, 6475, 6475, 6475 
## Resampling results across tuning parameters:
## 
##   max_depth  RMSE      Rsquared    MAE     
##   3          12.93823  0.07718419  7.052258
##   4          12.90191  0.08131200  7.011898
## 
## Tuning parameter 'nrounds' was held constant at a value of 100
## Tuning
##  parameter 'min_child_weight' was held constant at a value of 5
## 
## Tuning parameter 'subsample' was held constant at a value of 0.7
## RMSE was used to select the optimal model using the smallest value.
## The final values used for the model were nrounds = 100, max_depth = 4, eta
##  = 0.05, gamma = 1, colsample_bytree = 0.7, min_child_weight = 5 and
##  subsample = 0.7.

plot(varImp(xgb_risk), top= 20)

xgb_logit

## eXtreme Gradient Boosting 
## 
## 8093 samples
##   15 predictor
##    2 classes: 'Sí', 'No' 
## 
## Pre-processing: centered (20), scaled (20) 
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 6475, 6474, 6475, 6474, 6474 
## Resampling results across tuning parameters:
## 
##   eta   max_depth  nrounds  ROC        Sens       Spec     
##   0.05  3          100      0.6755658  0.7738602  0.4637363
##   0.05  3          200      0.6798934  0.7558945  0.4947802
##   0.05  5          100      0.6816291  0.7594868  0.4903846
##   0.05  5          200      0.6852564  0.7386070  0.5271978
##   0.10  3          100      0.6794372  0.7624066  0.4983516
##   0.10  3          200      0.6836734  0.7442190  0.5159341
##   0.10  5          100      0.6815769  0.7379296  0.5195055
##   0.10  5          200      0.6781645  0.7186204  0.5348901
## 
## Tuning parameter 'gamma' was held constant at a value of 0
## Tuning
## 
## Tuning parameter 'min_child_weight' was held constant at a value of 1
## 
## Tuning parameter 'subsample' was held constant at a value of 0.8
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were nrounds = 200, max_depth = 5, eta
##  = 0.05, gamma = 0, colsample_bytree = 0.8, min_child_weight = 1 and
##  subsample = 0.8.

plot(varImp(xgb_logit), top= 20)

Validación de los modelos (predicción)

Ahora que los modelos están entrenados, llega el paso más crucial: evaluar su rendimiento en un campo de juego neutral. Para ello, utilizo el conjunto de datos de prueba (df_test), que contiene información que los modelos nunca han visto durante su entrenamiento. Esto me permite obtener una medida honesta y realista de qué tan bien pueden generalizar sus predicciones a nuevos escenarios viales en Bogotá.

Para los modelos de regresión, calcularé el RMSE (Error Cuadrático Medio) y el R-cuadrado. Para los modelos de clasificación, generaré matrices de confusión para evaluar la precisión y la sensibilidad, y dibujaré una curva ROC comparativa para determinar qué modelo es el mejor para distinguir entre tramos con y sin siniestros.

# Función para evaluar modelos de regresión
evaluar_regresion <- function(modelo, datos_prueba, variable_objetivo) {
  predicciones <- predict(modelo, newdata = datos_prueba)
  postResample(pred = predicciones, obs = datos_prueba[[variable_objetivo]])
}

# --- 1. Evaluación de modelos para Velocidad Promedio ---
resultados_vel_lineal <- evaluar_regresion(m_vel, df_test, "vel_promed")
resultados_vel_rf    <- evaluar_regresion(rf_vel, df_test, "vel_promed")
resultados_vel_xgb   <- evaluar_regresion(xgb_vel, df_test, "vel_promed")

# --- 2. Evaluación de modelos para Índice de Exceso ---
resultados_exc_lineal <- evaluar_regresion(m_exc, df_test, "indice_exceso")
resultados_exc_rf    <- evaluar_regresion(rf_exc, df_test, "indice_exceso")
resultados_exc_xgb   <- evaluar_regresion(xgb_exc, df_test, "indice_exceso")

# --- 3. Evaluación de modelos para Índice de Riesgo ---
resultados_risk_lineal <- evaluar_regresion(m_risk, df_test, "indice_riesgo")
resultados_risk_rf    <- evaluar_regresion(rf_risk, df_test, "indice_riesgo")
resultados_risk_xgb   <- evaluar_regresion(xgb_risk, df_test, "indice_riesgo")

resultados_vel_lineal

##       RMSE   Rsquared        MAE 
## 10.1897855  0.2002195  8.1816445

resultados_vel_rf

##      RMSE  Rsquared       MAE 
## 8.0279455 0.5050522 5.8080060

resultados_vel_xgb

##      RMSE  Rsquared       MAE 
## 9.2143836 0.3572529 7.3046444

resultados_exc_lineal

##      RMSE  Rsquared       MAE 
## 0.2228191 0.2023049 0.1582864

resultados_exc_rf

##      RMSE  Rsquared       MAE 
## 0.1737999 0.5151767 0.1065803

resultados_exc_xgb

##      RMSE  Rsquared       MAE 
## 0.1998678 0.3688992 0.1394732

resultados_risk_lineal

##        RMSE    Rsquared         MAE 
## 13.76171944  0.04533893  7.42720812

resultados_risk_rf

##        RMSE    Rsquared         MAE 
## 13.55363900  0.07543141  7.19924564

resultados_risk_xgb

##        RMSE    Rsquared         MAE 
## 13.61296336  0.06569377  7.17262113

# --- Crear una tabla comparativa de resultados ---

tabla_resultados <- data.frame(
  Modelo = c("Velocidad promedio - Lineal", "Velocidad promedio - Random Forest", "Velocidad promedio - XGBoost",
             "Índice de exceso - Lineal", "Índice de exceso - Random Forest", "Índice de exceso - XGBoost",
             "Índice de riesgo - Lineal", "Índice de riesgo - Random Forest", "Índice de riesgo - XGBoost"),
  RMSE = c(resultados_vel_lineal["RMSE"], resultados_vel_rf["RMSE"], resultados_vel_xgb["RMSE"],
           resultados_exc_lineal["RMSE"], resultados_exc_rf["RMSE"], resultados_exc_xgb["RMSE"],
           resultados_risk_lineal["RMSE"], resultados_risk_rf["RMSE"], resultados_risk_xgb["RMSE"]),
  Rsquared = c(resultados_vel_lineal["Rsquared"], resultados_vel_rf["Rsquared"], resultados_vel_xgb["Rsquared"],
               resultados_exc_lineal["Rsquared"], resultados_exc_rf["Rsquared"], resultados_exc_xgb["Rsquared"],
               resultados_risk_lineal["Rsquared"], resultados_risk_rf["Rsquared"], resultados_risk_xgb["Rsquared"]),
  MAE = c(resultados_vel_lineal["MAE"], resultados_vel_rf["MAE"], resultados_vel_xgb["MAE"],
          resultados_exc_lineal["MAE"], resultados_exc_rf["MAE"], resultados_exc_xgb["MAE"],
          resultados_risk_lineal["MAE"], resultados_risk_rf["MAE"], resultados_risk_xgb["MAE"])
)


knitr::kable(tabla_resultados, caption = "Comparación de rendimiento de modelos de regresión en validación")

Comparación de rendimiento de modelos de regresión en validación

Modelo	RMSE	Rsquared	MAE
Velocidad promedio - Lineal	10.1897855	0.2002195	8.1816445
Velocidad promedio - Random Forest	8.0279455	0.5050522	5.8080060
Velocidad promedio - XGBoost	9.2143836	0.3572529	7.3046444
Índice de exceso - Lineal	0.2228191	0.2023049	0.1582864
Índice de exceso - Random Forest	0.1737999	0.5151767	0.1065803
Índice de exceso - XGBoost	0.1998678	0.3688992	0.1394732
Índice de riesgo - Lineal	13.7617194	0.0453389	7.4272081
Índice de riesgo - Random Forest	13.5536390	0.0754314	7.1992456
Índice de riesgo - XGBoost	13.6129634	0.0656938	7.1726211

# PASO CLAVE: Crear un set de prueba limpio y completo
# Para asegurar que los datos de predicción y los de referencia tengan la misma longitud,
# elimino cualquier fila de df_test que tenga NAs en las variables que el modelo necesita.
# Esto garantiza que predict() no omita filas silenciosamente.

# Extraemos las variables que el modelo glm REALMENTE necesita
vars_necesarias <- all.vars(formula(modelo_logit))

# Creamos el set de prueba limpio para este modelo
df_test_clean <- df_test %>%
  dplyr::select(all_of(vars_necesarias)) %>%
  tidyr::drop_na()

# Ahora, df_test_clean es nuestro set de prueba validado.
#---------------------------------------------------------

# --- 1. Generar predicciones para cada modelo USANDO EL DATASET LIMPIO ---

# Modelo Logístico (glm)
pred_clase_logit <- predict(modelo_logit, newdata = df_test_clean, type = "response")
pred_clase_logit <- factor(ifelse(pred_clase_logit > 0.5, "Sí", "No"), levels = c("Sí", "No"))
pred_prob_logit  <- predict(modelo_logit, newdata = df_test_clean, type = "response")

# Modelos de Caret (Random Forest y XGBoost)
pred_clase_rf <- predict(rf_log, newdata = df_test_clean)
pred_prob_rf  <- predict(rf_log, newdata = df_test_clean, type = "prob")

pred_clase_xgb <- predict(xgb_logit, newdata = df_test_clean)
pred_prob_xgb  <- predict(xgb_logit, newdata = df_test_clean, type = "prob")

# --- 2. Matriz de Confusión: Ahora las longitudes coincidirán ---

# La variable de referencia ahora viene del dataset limpio: df_test_clean$siniestros_dummy
cat("### Matriz de Confusión - Modelo Logístico\n")

## ### Matriz de Confusión - Modelo Logístico

confusionMatrix(data = pred_clase_logit, reference = df_test_clean$siniestros_dummy, positive = "Sí")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   Sí   No
##         Sí  521  726
##         No 1387  834
##                                           
##                Accuracy : 0.3907          
##                  95% CI : (0.3744, 0.4072)
##     No Information Rate : 0.5502          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : -0.1852         
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 0.2731          
##             Specificity : 0.5346          
##          Pos Pred Value : 0.4178          
##          Neg Pred Value : 0.3755          
##              Prevalence : 0.5502          
##          Detection Rate : 0.1502          
##    Detection Prevalence : 0.3596          
##       Balanced Accuracy : 0.4038          
##                                           
##        'Positive' Class : Sí              
## 

cat("\n### Matriz de Confusión - Random Forest\n")

## 
## ### Matriz de Confusión - Random Forest

confusionMatrix(data = pred_clase_rf, reference = df_test_clean$siniestros_dummy, positive = "Sí")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   Sí   No
##         Sí 1399  751
##         No  509  809
##                                           
##                Accuracy : 0.6367          
##                  95% CI : (0.6204, 0.6527)
##     No Information Rate : 0.5502          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.2554          
##                                           
##  Mcnemar's Test P-Value : 1.126e-11       
##                                           
##             Sensitivity : 0.7332          
##             Specificity : 0.5186          
##          Pos Pred Value : 0.6507          
##          Neg Pred Value : 0.6138          
##              Prevalence : 0.5502          
##          Detection Rate : 0.4034          
##    Detection Prevalence : 0.6200          
##       Balanced Accuracy : 0.6259          
##                                           
##        'Positive' Class : Sí              
## 

cat("\n### Matriz de Confusión - XGBoost\n")

## 
## ### Matriz de Confusión - XGBoost

confusionMatrix(data = pred_clase_xgb, reference = df_test_clean$siniestros_dummy, positive = "Sí")

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction   Sí   No
##         Sí 1387  776
##         No  521  784
##                                           
##                Accuracy : 0.626           
##                  95% CI : (0.6097, 0.6421)
##     No Information Rate : 0.5502          
##     P-Value [Acc > NIR] : < 2.2e-16       
##                                           
##                   Kappa : 0.233           
##                                           
##  Mcnemar's Test P-Value : 1.753e-12       
##                                           
##             Sensitivity : 0.7269          
##             Specificity : 0.5026          
##          Pos Pred Value : 0.6412          
##          Neg Pred Value : 0.6008          
##              Prevalence : 0.5502          
##          Detection Rate : 0.3999          
##    Detection Prevalence : 0.6237          
##       Balanced Accuracy : 0.6148          
##                                           
##        'Positive' Class : Sí              
## 

# --- 3. Curva ROC y AUC ---

# La variable de referencia también viene del dataset limpio
roc_logit <- roc(df_test_clean$siniestros_dummy, pred_prob_logit)

## Setting levels: control = Sí, case = No

## Setting direction: controls < cases

roc_rf    <- roc(df_test_clean$siniestros_dummy, pred_prob_rf$Sí)

## Setting levels: control = Sí, case = No

## Setting direction: controls > cases

roc_xgb   <- roc(df_test_clean$siniestros_dummy, pred_prob_xgb$Sí)

## Setting levels: control = Sí, case = No
## Setting direction: controls > cases

# Imprimir los valores de AUC
cat(paste0("\nAUC Modelo Logístico: ", round(auc(roc_logit), 4)))

## 
## AUC Modelo Logístico: 0.6328

cat(paste0("\nAUC Random Forest: ", round(auc(roc_rf), 4)))

## 
## AUC Random Forest: 0.6754

cat(paste0("\nAUC XGBoost: ", round(auc(roc_xgb), 4)))

## 
## AUC XGBoost: 0.6696

# Graficar las curvas ROC juntas para comparar visualmente
ggroc(list(Logit = roc_logit, RF = roc_rf, XGBoost = roc_xgb)) +
  ggtitle("Curva ROC Comparativa de Modelos de Clasificación") +
  theme_minimal()

# --- 2. Recopilar métricas de CLASIFICACIÓN ---
# Extraemos los resultados de las matrices de confusión y los cálculos de AUC

# Modelo Logístico
cm_logit <- confusionMatrix(data = pred_clase_logit, reference = df_test_clean$siniestros_dummy, positive = "Sí")
auc_logit <- auc(roc_logit)

# Random Forest
cm_rf <- confusionMatrix(data = pred_clase_rf, reference = df_test_clean$siniestros_dummy, positive = "Sí")
auc_rf <- auc(roc_rf)

# XGBoost
cm_xgb <- confusionMatrix(data = pred_clase_xgb, reference = df_test_clean$siniestros_dummy, positive = "Sí")
auc_xgb <- auc(roc_xgb)

# Crear el data.frame para los modelos de clasificación
tabla_resultados_clf <- data.frame(
  Modelo = c("Siniestros - Logístico", "Siniestros - Random Forest", "Siniestros - XGBoost"),
  Accuracy = c(cm_logit$overall["Accuracy"], cm_rf$overall["Accuracy"], cm_xgb$overall["Accuracy"]),
  Sensitivity = c(cm_logit$byClass["Sensitivity"], cm_rf$byClass["Sensitivity"], cm_xgb$byClass["Sensitivity"]),
  Specificity = c(cm_logit$byClass["Specificity"], cm_rf$byClass["Specificity"], cm_xgb$byClass["Specificity"]),
  AUC = c(auc_logit, auc_rf, auc_xgb)
)


# --- 3. Unir ambas tablas en una tabla final ---
# dplyr::bind_rows es perfecto para esto, ya que rellena con NA las columnas que no coinciden.
tabla_final <- dplyr::bind_rows(tabla_resultados, tabla_resultados_clf)

# Reordenar las columnas para una mejor presentación
tabla_final <- tabla_final %>%
  dplyr::select(Modelo, RMSE, Rsquared, MAE, Accuracy, Sensitivity, Specificity, AUC)

# --- 4. Presentar la tabla con formato profesional ---
library(kableExtra)

## Warning: package 'kableExtra' was built under R version 4.4.3

## 
## Attaching package: 'kableExtra'

## The following object is masked from 'package:dplyr':
## 
##     group_rows

knitr::kable(
  tabla_final, 
  caption = "Comparación de Rendimiento de Todos los Modelos en Validación",
  digits = 4 # Redondear a 4 decimales
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  # Colorear las celdas NA para que sean menos distractoras
  column_spec(2:4, background = ifelse(is.na(tabla_final$RMSE), "#F0F0F0", "white")) %>%
  column_spec(5:8, background = ifelse(is.na(tabla_final$Accuracy), "#F0F0F0", "white"))

Comparación de Rendimiento de Todos los Modelos en Validación

Modelo	RMSE	Rsquared	MAE	Accuracy	Sensitivity	Specificity	AUC
Velocidad promedio - Lineal	10.1898	0.2002	8.1816	NA	NA	NA	NA
Velocidad promedio - Random Forest	8.0279	0.5051	5.8080	NA	NA	NA	NA
Velocidad promedio - XGBoost	9.2144	0.3573	7.3046	NA	NA	NA	NA
Índice de exceso - Lineal	0.2228	0.2023	0.1583	NA	NA	NA	NA
Índice de exceso - Random Forest	0.1738	0.5152	0.1066	NA	NA	NA	NA
Índice de exceso - XGBoost	0.1999	0.3689	0.1395	NA	NA	NA	NA
Índice de riesgo - Lineal	13.7617	0.0453	7.4272	NA	NA	NA	NA
Índice de riesgo - Random Forest	13.5536	0.0754	7.1992	NA	NA	NA	NA
Índice de riesgo - XGBoost	13.6130	0.0657	7.1726	NA	NA	NA	NA
Siniestros - Logístico	NA	NA	NA	0.3907	0.2731	0.5346	0.6328
Siniestros - Random Forest	NA	NA	NA	0.6367	0.7332	0.5186	0.6754
Siniestros - XGBoost	NA	NA	NA	0.6260	0.7269	0.5026	0.6696

Visualizaciones para el modelo Random Forest del logit con la variable Dummy de siniestralidad

Cada gráfico a continuación muestra el efecto que tiene una sola variable predictora sobre la probabilidad de que ocurra un siniestro, mientras se mantienen todas las demás variables constantes en su valor promedio. Esto nos permite aislar y visualizar el impacto marginal de cada factor clave. Por ejemplo, podemos ver cómo cambia el riesgo a medida que aumenta el número de árboles o la distancia a la cámara más cercana.

# Visualización del efecto parcial de las variables más importantes en el modelo Random Forest
# dist_min_camara
pdp_camara <- partial(rf_log, pred.var = "dist_min_camara", prob = TRUE, which.class = "Sí")
# vel_promed
pdp_vel <- partial(rf_log, pred.var = "vel_promed", prob = TRUE, which.class = "Sí")
# indice_exceso
pdp_exceso <- partial(rf_log, pred.var = "indice_exceso", prob = TRUE, which.class = "Sí")
# mean_v_ref
pdp_vref <- partial(rf_log, pred.var = "mean_v_ref", prob = TRUE, which.class = "Sí")
# mean_denpo
pdp_denpo <- partial(rf_log, pred.var = "mean_denpo", prob = TRUE, which.class = "Sí")
# ancho
pdp_ancho <- partial(rf_log, pred.var = "ancho", prob = TRUE, which.class = "Sí")
# rutas_sitp
pdp_rutas <- partial(rf_log, pred.var = "rutas_sitp", prob = TRUE, which.class = "Sí")
# n_arboles
pdp_arboles <- partial(rf_log, pred.var = "n_arboles", prob = TRUE, which.class = "Sí")
# n_semaforo
pdp_semaforo <- partial(rf_log, pred.var = "n_semaforo", prob = TRUE, which.class = "Sí")
# carriles
pdp_carriles <- partial(rf_log, pred.var = "carriles", prob = TRUE, which.class = "Sí")

# --- Gráfico 1: Distancia a la cámara ---
ggplot(pdp_camara, aes(x = dist_min_camara, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial de la distancia a la cámara más cercana",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Distancia a la cámara más cercana (metros)",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal(base_size = 14)

## 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.

# --- Gráfico 2: Velocidad Promedio ---
ggplot(pdp_vel, aes(x = vel_promed, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial de la velocidad promedio",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Velocidad promedio (km/h)",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal(base_size = 14)

# --- Gráfico 3: Índice de Exceso de Velocidad ---
ggplot(pdp_exceso, aes(x = indice_exceso, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del índice de exceso de velocidad",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Índice de exceso de velocidad (Proporción en 24h)",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(labels = scales::percent) + # También en porcentaje
  theme_minimal(base_size = 14)

# --- Gráfico 4: Valor Catastral ---
ggplot(pdp_vref, aes(x = mean_v_ref, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del valor catastral promedio de la manzana",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Valor catastral promedio (Miles de milllones COP)",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  # Escala eje x en millones
  scale_x_continuous(labels = scales::dollar_format(prefix = "$", suffix = "M", scale = 1e-6),
                     breaks = scales::pretty_breaks(n = 10)) +
  theme_minimal(base_size = 14)

# --- Gráfico 5: Densidad Poblacional ---
ggplot(pdp_denpo, aes(x = mean_denpo, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial de la densidad poblacional de la manzana",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Densidad poblacional promedio (habs/ha)",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal(base_size = 14)

# --- Gráfico 6: Ancho de la calzada ---
ggplot(pdp_ancho, aes(x = ancho, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del ancho de la calzada",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Ancho (metros)",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal(base_size = 14)

# --- Gráfico 7: Rutas SITP ---
ggplot(pdp_rutas, aes(x = rutas_sitp, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del número de rutas SITP por el tramo vial",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Número de rutas SITP en el tramo",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal(base_size = 14)

# --- Gráfico 8: Número de Árboles ---
ggplot(pdp_arboles, aes(x = n_arboles, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del número de árboles en el tramo vial",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Número de árboles en el tramo",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal(base_size = 14)

# --- Gráfico 9: Número de Semáforos ---
ggplot(pdp_semaforo, aes(x = n_semaforo, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del número de semáforos en el tramo vial",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Número de semáforos en el tramo",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) +
  theme_minimal(base_size = 14)

# --- Gráfico 10: Número de Carriles ---
ggplot(pdp_carriles, aes(x = carriles, y = yhat)) +
  geom_line(color = "#0072B2", size = 1.2) +
  labs(
    title = "Efecto parcial del número de carriles de la calzada",
    subtitle = "Influencia sobre la probabilidad de que ocurra un siniestro",
    x = "Número de carriles en la calzada",
    y = "Probabilidad estimada"
  ) +
  scale_y_continuous(labels = scales::percent) +
  # Escala eje x solo números enteros
  scale_x_continuous(breaks = scales::pretty_breaks(n = 6)) +
  theme_minimal(base_size = 14)

Visualizaciones para el modelo de velocidad promedio (rf_vel)

Los siguientes gráficos muestran cómo cambia la velocidad promedio predicha al variar un predictor, manteniendo el resto constantes.

# Se calculan los PDP para el modelo de velocidad
pdp_vel_camara <- partial(rf_vel, pred.var = "dist_min_camara")
pdp_vel_denem <- partial(rf_vel, pred.var = "mean_denem")
pdp_vel_vref <- partial(rf_vel, pred.var = "mean_v_ref")
pdp_vel_denpo <- partial(rf_vel, pred.var = "mean_denpo")

# Se grafican usando la función personalizada (si no la tienes, créala)
# O usa este código directo con ggplot2
ggplot(pdp_vel_camara, aes(x = dist_min_camara, y = yhat)) + 
  geom_line(color = "#0072B2", size = 1.2) + 
  labs(title = "Efecto de la distancia a la cámara más cercana", 
       subtitle = "Influencia sobre la velocidad promedio",
       x = "Distancia a la cámara más cercana (m)", 
       y = "Velocidad Promedio (km/h)") + 
  theme_minimal()

ggplot(pdp_vel_denem, aes(x = mean_denem, y = yhat)) + 
  geom_line(color = "#0072B2", size = 1.2) + 
  labs(title = "Efecto de la Densidad de Empleos", 
       subtitle = "Influencia sobre la velocidad promedio",
       x = "Densidad de empleos", 
       y = "Velocidad Promedio (km/h)") + 
  theme_minimal()

ggplot(pdp_vel_vref, aes(x = mean_v_ref, y = yhat)) + 
  geom_line(color = "#0072B2", size = 1.2) + 
  labs(title = "Efecto del valor catastral promedio de la manzana", 
       subtitle = "Influencia sobre la velocidad promedio",
       x = "Valor Catastral Promedio (Miles de milllones COP)", 
       y = "Velocidad Promedio (km/h)") +
  scale_x_continuous(labels = scales::dollar_format(scale = 1e-6, suffix = "M")) +
  theme_minimal()

ggplot(pdp_vel_denpo, aes(x = mean_denpo, y = yhat)) + 
  geom_line(color = "#0072B2", size = 1.2) + 
  labs(title = "Efecto de la densidad poblacional de la manzana", 
       subtitle = "Influencia sobre la velocidad promedio",
       x = "Densidad Poblacional (habs/ha)", 
       y = "Velocidad Promedio (km/h)") + 
  theme_minimal()

Visualizaciones para el Modelo de Índice de Exceso (rf_exc)

De manera similar, se analiza el impacto de las variables más relevantes sobre la proporción de tiempo en que los conductores exceden el límite de velocidad.

# Se calculan los PDP para el modelo de índice de exceso
pdp_exc_camara <- partial(rf_exc, pred.var = "dist_min_camara")
pdp_exc_vref <- partial(rf_exc, pred.var = "mean_v_ref")
pdp_exc_denpo <- partial(rf_exc, pred.var = "mean_denpo")
pdp_exc_denem <- partial(rf_exc, pred.var = "mean_denem")

# Se grafican
ggplot(pdp_exc_camara, aes(x = dist_min_camara, y = yhat)) + 
  geom_line(color = "#D55E00", size = 1.2) + 
  labs(title = "Efecto de la distancia a la cámara más cercana",
       subtitle = "Influencia sobre el índice de exceso de velocidad",
       x = "Distancia a la cámara más cercana (m)", 
       y = "Índice de Exceso (24h)") + 
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

ggplot(pdp_exc_vref, aes(x = mean_v_ref, y = yhat)) + 
  geom_line(color = "#D55E00", size = 1.2) + 
  labs(title = "Efecto del valor catastral promedio de la manzana",
       subtitle = "Influencia sobre el índice de exceso de velocidad",
       x = "Valor Catastral Promedio (Miles de milllones COP)", 
       y = "Índice de Exceso") +
  scale_x_continuous(labels = scales::dollar_format(scale = 1e-6, suffix = "M")) +
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

ggplot(pdp_exc_denpo, aes(x = mean_denpo, y = yhat)) + 
  geom_line(color = "#D55E00", size = 1.2) + 
  labs(title = "Efecto de la densidad poblacional de la manzana", 
       subtitle = "Influencia sobre el índice de exceso de velocidad",
       x = "Densidad Poblacional (habs/ha)", 
       y = "Índice de Exceso") + 
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

ggplot(pdp_exc_denem, aes(x = mean_denem, y = yhat)) + 
  geom_line(color = "#D55E00", size = 1.2) + 
  labs(title = "Efecto de la Densidad de Empleos", 
       subtitle = "Influencia sobre el índice de exceso de velocidad",
       x = "Densidad de Empleos (promedio)", 
       y = "Índice de Exceso") + 
  scale_y_continuous(labels = scales::percent) +
  theme_minimal()

Resultados finales

# Predecir la probabilidad de siniestro para todos los tramos viales de Bogotá
predicciones_finales <- predict(rf_log, newdata = vias, type = "prob")

# Crear un nuevo dataframe con los resultados
resultados_finales <- data.frame(
  pk_id_calz = vias$pk_id_calz,
  nombre_corredor = vias$tipos_via, 
  probabilidad_siniestro = predicciones_finales$Sí,
  siniestros_reportados = vias$n_her + vias$n_muer, # Para comparar con la realidad
  tipomallae = vias$tipomallae 
) %>%
  arrange(desc(probabilidad_siniestro)) # Ordenar de mayor a menor riesgo

# Ver los 20 tramos viales con mayor probabilidad de siniestro
head(resultados_finales, 20)

##    pk_id_calz                            nombre_corredor probabilidad_siniestro
## 1      517598      AVENIDA AGOBERTO MEJIA CIFUENTES (  )              0.9171496
## 2      412253               AVENIDA PRIMERO DE MAYO (  )              0.9167669
## 3      511884      AVENIDA AGOBERTO MEJIA CIFUENTES (  )              0.9166820
## 4    24122041               AVENIDA PRIMERO DE MAYO (  )              0.9161936
## 5      515526                AVENIDA CIUDAD DE CALI (  )              0.9149267
## 6      412009               AVENIDA PRIMERO DE MAYO (  )              0.9124526
## 7    91027062               AVENIDA DE LAS AMERICAS (  )              0.9123372
## 8      530780                AVENIDA CIUDAD DE CALI (  )              0.9119728
## 9    91027053               AVENIDA DE LAS AMERICAS (  )              0.9111455
## 10   91027054               AVENIDA DE LAS AMERICAS (  )              0.9088832
## 11     506787           AVENIDA CIUDAD DE LIMA ( CL 19 )              0.9085982
## 12     515718                AVENIDA CIUDAD DE CALI (  )              0.9068099
## 13     518663                AVENIDA CIUDAD DE CALI (  )              0.9052888
## 14     472139 AVENIDA DEL CONGRESO EUCARISTICO ( AK 68 )              0.9051555
## 15   91027063               AVENIDA DE LAS AMERICAS (  )              0.9045948
## 16     506785           AVENIDA CIUDAD DE LIMA ( CL 19 )              0.9041367
## 17     506797           AVENIDA CIUDAD DE LIMA ( CL 19 )              0.9039677
## 18     515500                AVENIDA CIUDAD DE CALI (  )              0.9030491
## 19   91012171    AVENIDA ALBERTO LLERAS CAMARGO ( AK 7 )              0.9029581
## 20     515519                AVENIDA CIUDAD DE CALI (  )              0.9026675
##    siniestros_reportados tipomallae
## 1                     10   Arterial
## 2                      4   Arterial
## 3                      7   Arterial
## 4                     14   Arterial
## 5                      1   Arterial
## 6                      8   Arterial
## 7                      2    Troncal
## 8                     16   Arterial
## 9                     12    Troncal
## 10                     3    Troncal
## 11                     2   Arterial
## 12                     0   Arterial
## 13                     8   Arterial
## 14                     5   Arterial
## 15                    17    Troncal
## 16                     7   Arterial
## 17                    13   Arterial
## 18                     3   Arterial
## 19                     6    Troncal
## 20                     7   Arterial

Cierre del clúster paralelo

stopCluster(cl)