Modelo RF FIA
Tomás Villescas
2025-09-25
Modelo FIA: Random Forest
Este documento presenta el desarrollo y la validación de un conjunto de modelos predictivos de Machine Learning, enmarcados en el proyecto FIA Speed Camera Enforcement. El objetivo principal de este análisis es identificar y cuantificar los factores de infraestructura, operacionales y socioeconómicos que influyen en la seguridad y operación vial en la ciudad de Bogotá. Para ello, se busca predecir cuatro variables clave en los tramos viales de la ciudad: la velocidad de operación (medida con Google API), la proporción del exceso de velocidad en la velocidad promedio de una vía (proporción de 24h), un índice de riesgo ponderado por siniestralidad y, fundamentalmente, la probabilidad de que ocurra un siniestro vial con heridos o fallecidos.
La base de datos utilizada para este fin es una tabla georreferenciada que consolida información detallada de 11,561 tramos viales y sus manzanas adyacentes. Esta fuente de datos, proveniente de diversas entidades y procesada para este proyecto, integra variables de diseño geométrico (ancho, número de carriles), de entorno (densidad poblacional y de empleos, valor catastral), operacionales (velocidades por hora, presencia de rutas SITP) y de seguridad (histórico de siniestros). El enfoque metodológico se centra en el algoritmo Random Forest, reconocido por su alta precisión predictiva y su capacidad para modelar relaciones complejas y no lineales.
Librerías y datos
La primera fase del análisis consiste en la configuración del entorno
de trabajo en R. Esto incluye la carga de todas las librerías necesarias
para el procesamiento de datos, la modelación y la visualización. Se
establece una semilla de reproducibilidad (set.seed) para
garantizar que los resultados que dependen de procesos aleatorios, como
la partición de datos y el entrenamiento de los modelos, sean
consistentes entre ejecuciones.
Adicionalmente, se configura un clúster de procesamiento en paralelo
para acelerar significativamente los cálculos computacionalmente
intensivos, como el ajuste de hiperparámetros y el entrenamiento de los
modelos de Random Forest. Para este informe en particular, se cargan los
resultados de los modelos ya entrenados desde un archivo
(modelos_entrenados.RData) para permitir una generación
rápida del documento sin necesidad de volver a ejecutar los cálculos
pesados.
library(tidyverse)
library(readxl)
library(janitor)
library(caret)
library(pROC)
library(ranger)
library(doParallel)
library(pdp)
library(kableExtra)
set.seed(123)
# --- Configuración Robusta para Paralelización y Knitr ---
# 1. Crear el clúster
cl <- makePSOCKcluster(max(1, parallel::detectCores() - 1))
registerDoParallel(cl)
# 2. (SOLUCIÓN) Cargar las librerías necesarias en cada proceso del clúster
# Los workers no heredan las librerías automáticamente al tejer
clusterEvalQ(cl, {
library(caret)
library(ranger)
})## [[1]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[2]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[3]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[6]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[7]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[8]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[9]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[10]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[11]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[12]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[13]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[14]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"
##
## [[15]]
## [1] "ranger" "caret" "lattice" "ggplot2" "stats" "graphics"
## [7] "grDevices" "utils" "datasets" "methods" "base"# 3. Asegurar que el clúster se detenga al final, incluso si hay un error
on.exit({
try(stopCluster(cl), silent = TRUE)
try(registerDoSEQ(), silent = TRUE) # Vuelve a modo no-paralelo
}, add = TRUE)
# --- Cargar modelos entrenados previamente ---
load("modelos_entrenados.RData")El éxito de cualquier modelo de Machine Learning depende fundamentalmente de la calidad de los datos de entrada. Por ello, se realiza un proceso exhaustivo de ingeniería y limpieza de variables. Este paso es crucial para transformar los datos crudos en un formato estructurado y optimizado para el algoritmo de Random Forest.
Las tareas realizadas en este bloque incluyen:
Limpieza de Nombres: Estandarización de los nombres de las columnas a un formato consistente.
Conversión de Tipos: Asegurar que las variables numéricas y categóricas tengan el tipo de dato correcto.
Creación de Variables Derivadas: Generación de nuevas variables con mayor poder predictivo a partir de las existentes. Esto incluye la creación de un indice_exceso estandarizado, una variable binaria (
siniestros_dummy) que indica la presencia o ausencia de siniestros y la variabledist_min_camaraque calcula la distancia a la cámara de fotodetección más cercana.Manejo de Datos Faltantes: Identificación y tratamiento de valores
SIN DATO.Formato de Factores: Conversión de variables categóricas clave a factores y establecimiento de niveles de referencia para una correcta interpretación por parte del modelo.
# --- Cargar Datos ---
vias <- read_xlsx("Tabla_corregida_poligonos.xlsx") %>%
clean_names %>%
mutate(across(where(is.character), ~na_if(.x, "SIN DATO")))
# --- Ingeniería y Limpieza de Variables ---
vias <- vias %>%
mutate(
# --- Variables Derivadas ---
n_muer = as.integer(n_muer),
n_her = as.integer(n_her),
estrato_fi = as.integer(estrato_fi),
# Índice de exceso: proporción de tiempo sobre 50 km/h (limitado entre 0 y 1)
# Se asume que la variable 'exceso' ya es la proporción. Si no, ajustar.
indice_exceso = pmin(pmax(exceso, 0), 1),
# Dummy de siniestros (numérico)
siniestros_dummy_num = as.integer(replace_na(n_muer, 0) + replace_na(n_her, 0) > 0),
# --- CORRECCIÓN: Variable de distancia a la cámara más cercana ---
# Se usa pmin para encontrar la mínima distancia entre todas las cámaras
dist_min_camara = suppressWarnings(pmin(camara_1, camara_2, camara_3, camara_4, camara_5, na.rm = TRUE)),
# --- Limpieza de Datos ---
estaashto = if_else(estaashto == "SIN DATO", "BUENO", estaashto),
# --- Factores y Niveles de Referencia ---
# Se convierten las variables categóricas a factores para el modelo
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"),
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: factor con niveles "Sí"/"No" ---
siniestros_dummy = factor(if_else(siniestros_dummy_num == 1, "Sí", "No"), levels = c("Sí", "No"))
)
vias <- vias %>%
mutate(
indice_exceso = case_when(
"exceso" %in% names(vias) ~ pmin(pmax(as.numeric(exceso), 0), 1),
"hora_mas_5" %in% names(vias) ~ pmin(pmax(as.numeric(hora_mas_5)/24, 0), 1),
TRUE ~ NA_real_
)
)Definición de fórmulas y datos para el modelo
En este bloque, realizamos la partición de los datos en conjuntos de
entrenamiento (df_train) y prueba (df_test).
También realizamos una imputación de medianas para manejar los valores
faltantes de manera robusta. Estos conjuntos de datos serán la base para
todo el análisis posterior.
# --- Núcleo común de predictores ---
x_comunes <- as.formula("~ tipomallae + ancho + carriles + ciclorruta +
mean_v_ref + mean_denpo + mean_denem + n_arboles +
n_semaforo + rutas_sitp + estaashto +
dist_min_camara + longitud_k + ancho_carril")
# --- Partición de datos (70% entrenamiento, 30% prueba) ---
set.seed(123)
idx_split <- createDataPartition(vias$siniestros_dummy, p = .7, list = FALSE)
df_train <- vias[idx_split, ]
df_test <- vias[-idx_split, ]
# --- Imputación de Medianas ---
# Se calculan las medianas SOLO en el set de entrenamiento
num_cols <- names(select(df_train, where(is.numeric)))
medianas_train <- sapply(df_train[num_cols], \(x) median(x, na.rm = TRUE))
df_train[num_cols] <- lapply(df_train[num_cols], \(x) ifelse(is.na(x), medianas_train[names(x)], x))
df_test[num_cols] <- lapply(df_test[num_cols], \(x) ifelse(is.na(x), medianas_train[names(x)], x))
# --- Función para preparar dataframes limpios por modelo ---
prep_df <- function(data, form){
vars <- all.vars(form)
faltan <- setdiff(vars, names(data))
if (length(faltan)) stop("Faltan variables en la BD: ", paste(faltan, collapse = ", "))
data %>% dplyr::select(all_of(vars)) %>% tidyr::drop_na()
}
# --- Controles de entrenamiento ---
# Para modelos de Regresión
ctrl_reg <- trainControl(method = "cv", number = 5, verboseIter = FALSE, allowParallel = TRUE)
# Para modelos de Clasificación
ctrl_clf <- trainControl(method = "cv", number = 5,
classProbs = TRUE, summaryFunction = twoClassSummary,
verboseIter = FALSE, allowParallel = TRUE)
# --- Grilla de parámetros más exigente para Random Forest ---
tuneGrid_rf <- expand.grid(
.mtry = c(5, 10, 15),
.splitrule = "variance",
.min.node.size = c(15, 25)
)
# Grilla RF (clasificación) — más consistente que tuneLength
tuneGrid_rf_clf <- expand.grid(
.mtry = c(5, 10, 15),
.splitrule = c("gini","extratrees"),
.min.node.size = c(1, 5, 10)
)Análisis de sensibilidad
Antes de construir los modelos finales, realizamos un análisis de sensibilidad sistemático para tomar decisiones informadas sobre tres componentes clave:
- Métrica de Velocidad: Determinamos cuál de las métricas disponibles
(
percen_70,vel_promed,median_vel) tiene mayor poder predictivo sobre la seguridad vial. - Ponderación del Riesgo: Evaluamos diferentes formas de calcular el
indice_riesgopara encontrar la definición que mejor se alinea con las características viales. - Selección de Predictores: Utilizamos una técnica avanzada de Machine Learning, Recursive Feature Elimination (RFE), para identificar el subconjunto de variables más potente para cada uno de nuestros cuatro modelos objetivo.
Este proceso garantiza que los modelos finales se construyan sobre la base más sólida y empíricamente validada posible.
Sensibilidad para la métrica de velocidad
Evaluamos cómo cambia el rendimiento de los modelos de Siniestros y Riesgo al utilizar diferentes métricas de velocidad. El objetivo es identificar la variable de velocidad que maximiza el AUC (para clasificación) y el R-cuadrado (para regresión).
# --- ANÁLISIS DE SENSIBILIDAD 1: MÉTRICA DE VELOCIDAD ---
# Definimos las variables a probar y preparamos una lista para los resultados
velocity_vars <- c("percen_70", "vel_promed", "median_vel")
results_list_vel <- list()
# Bucle para probar cada variable de velocidad
for (vel_var in velocity_vars) {
cat("\n--- Probando métrica de velocidad:", vel_var, "---\n")
# Fórmulas dinámicas para riesgo y siniestros
form_risk_temp <- reformulate(c(all.vars(x_comunes)[-1], vel_var, "indice_exceso"), response = "indice_riesgo")
form_logit_temp <- reformulate(c(all.vars(x_comunes)[-1], vel_var, "indice_exceso"), response = "siniestros_dummy")
# Esto garantiza que no se pasen NAs a la función train()
train_risk_temp_clean <- prep_df(df_train, form_risk_temp)
train_log_temp_clean <- prep_df(df_train, form_logit_temp)
# Entrenamos modelos más rápidos para este test
rf_risk_temp <- train(form_risk_temp, data = train_risk_temp_clean, method = "ranger", trControl = ctrl_reg, tuneGrid = tuneGrid_rf, metric = "RMSE", num.trees = 500)
rf_log_temp <- train(form_logit_temp, data = train_log_temp_clean, method = "ranger", trControl = ctrl_clf, tuneGrid = tuneGrid_rf_clf, metric = "ROC", num.trees = 500)
# Evaluamos usando los datos de prueba limpios para cada fórmula
df_test_risk_clean <- prep_df(df_test, form_risk_temp)
df_test_log_clean <- prep_df(df_test, form_logit_temp)
res_risk <- postResample(predict(rf_risk_temp, df_test_risk_clean), df_test_risk_clean$indice_riesgo)
pred_prob <- predict(rf_log_temp, df_test_log_clean, type = "prob")
roc_obj <- roc(response = df_test_log_clean$siniestros_dummy, predictor = pred_prob$Sí, levels = c("No", "Sí"))
# Guardamos los resultados
results_list_vel[[vel_var]] <- data.frame(
Metrica_Velocidad = vel_var,
Rsquared_Risk = res_risk["Rsquared"],
AUC_Siniestros = as.numeric(auc(roc_obj))
)
}Se presentan los resultados.
# Combinamos y mostramos los resultados en una tabla limpia
sensitivity_results_vel <- bind_rows(results_list_vel)
kable(sensitivity_results_vel, digits = 4, caption = "Rendimiento de Modelos por Métrica de Velocidad") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Metrica_Velocidad | Rsquared_Risk | AUC_Siniestros | |
|---|---|---|---|
| Rsquared…1 | percen_70 | 0.0356 | 0.6089 |
| Rsquared…2 | vel_promed | 0.0334 | 0.6089 |
| Rsquared…3 | median_vel | 0.0361 | 0.6088 |
# Decisión automática basada en los resultados para usarla después
BEST_VEL_VAR <- sensitivity_results_vel %>%
slice_max(AUC_Siniestros, n = 1) %>%
pull(Metrica_Velocidad)
cat("\nDecisión: La mejor métrica de velocidad seleccionada es '", BEST_VEL_VAR, "'.\n", sep = "")##
## Decisión: La mejor métrica de velocidad seleccionada es 'vel_promed'.Sensibilidad para el índice de riesgo
Utilizando la mejor métrica de velocidad encontrada, probamos
diferentes ponderaciones para la fórmula del indice_riesgo
(muertos * peso + heridos), incluyendo una relación 1:1 que
simplemente suma el total de víctimas. Las alternativas utilizadas se
basan en la valoración económica de siniestralidad hecha por Bocarejo y
Díaz en 2010 (30.25:1), la valoración utilizada por el puntaje de los
operadores de TransMilenio (9:1) y la utilizada por la Secretaría
Distrital de Movilidad en su Potencial de Vidas Salvadas (6:1).
# --- ANÁLISIS DE SENSIBILIDAD 2: PONDERACIÓN DEL ÍNDICE DE RIESGO ---
# Definimos los pesos a probar, incluyendo la relación 1:1
risk_weights <- c(30.25, 9, 6, 1)
results_risk_list <- list()
for (weight in risk_weights) {
cat("\n--- Probando ponderación de riesgo:", weight, ":1 ---\n")
# Recalculamos el índice de riesgo con el nuevo peso
df_train_temp <- df_train %>% mutate(indice_riesgo_temp = n_muer * weight + n_her)
df_test_temp <- df_test %>% mutate(indice_riesgo_temp = n_muer * weight + n_her)
# Creamos la fórmula usando la mejor métrica de velocidad
form_risk_temp <- reformulate(c(all.vars(x_comunes)[-1], BEST_VEL_VAR, "indice_exceso"), response = "indice_riesgo_temp")
# --- CORRECCIÓN: Crear dataframe de entrenamiento limpio ANTES de entrenar ---
train_risk_temp_clean <- prep_df(df_train_temp, form_risk_temp)
# Entrenamos y evaluamos
rf_risk_temp <- train(form_risk_temp, data = train_risk_temp_clean, method = "ranger", trControl = ctrl_reg, tuneGrid = tuneGrid_rf, metric = "RMSE", num.trees = 500)
df_test_risk_clean <- prep_df(df_test_temp, form_risk_temp)
res_risk <- postResample(predict(rf_risk_temp, df_test_risk_clean), df_test_risk_clean$indice_riesgo_temp)
# Guardamos resultados
results_risk_list[[as.character(weight)]] <- data.frame(
Ponderacion = paste0(weight, ":1"),
RMSE = res_risk["RMSE"],
Rsquared = res_risk["Rsquared"]
)
}Se presentan los resultados.
# Mostramos los resultados
sensitivity_results_risk <- bind_rows(results_risk_list)
kable(sensitivity_results_risk, digits = 4, caption = "Rendimiento del Modelo de Riesgo por Ponderación") %>%
kable_styling(bootstrap_options = "striped", full_width = FALSE)| Ponderacion | RMSE | Rsquared | |
|---|---|---|---|
| RMSE…1 | 30.25:1 | 13.8893 | 0.0347 |
| RMSE…2 | 9:1 | 6.6219 | 0.0521 |
| RMSE…3 | 6:1 | 5.7787 | 0.0552 |
| RMSE…4 | 1:1 | 4.6901 | 0.0485 |
# Decisión automática de la mejor ponderación
BEST_RISK_WEIGHT <- sensitivity_results_risk %>%
slice_max(Rsquared, n = 1) %>%
pull(Ponderacion) %>%
sub(":1", "", .) %>%
as.numeric()
cat("\nDecisión: La mejor ponderación para el Índice de Riesgo es '", BEST_RISK_WEIGHT, ":1'.\n", sep = "")##
## Decisión: La mejor ponderación para el Índice de Riesgo es '6:1'.RFE para selección de variables
Finalmente, identificamos el subconjunto óptimo de predictores para cada uno de los cuatro modelos. RFE entrena modelos de forma iterativa, eliminando las variables menos importantes en cada paso para encontrar el balance perfecto entre simplicidad y poder predictivo.
# --- ANÁLISIS DE SENSIBILIDAD 3: SELECCIÓN DE PREDICTORES CON RFE ---
# Paso 1: Definir el universo de predictores, incluyendo la nueva 'ancho_carril'
all_candidate_predictors <- c(
"tipomallae", "estaashto", "aashto_bue", "estrato_fi", "mean_v_ref",
"mean_denpo", "mean_denem", "mean_estra", "tipo", "carriles", "sitp",
"pres_anden", "n_arboles", "n_semaforo", "rutas_sitp", "ciclorruta",
"flujo_1", "camara_100", "camara_200", "camara_300", "camara_500",
"longitud_k", "vel_promed", "percen_70", "median_vel",
"indice_exceso", "dist_min_camara", "ancho_carril"
)
# Paso 2: Definir el control de RFE (reutilizable)
rfe_control <- rfeControl(functions = rfFuncs, method = "cv", number = 5, verbose = FALSE) # verbose=FALSE para un output más limpio
subset_sizes <- c(5, 8, 10, 12, 15, 20, 25)
# Función auxiliar para ejecutar un ciclo de RFE
run_rfe <- function(data, predictors, response_var, sizes, control) {
cat("\n--- Iniciando RFE para:", response_var, "---\n")
train_data <- data %>% select(all_of(predictors), all_of(response_var)) %>% drop_na()
x_train <- train_data %>% select(-all_of(response_var))
y_train <- train_data[[response_var]]
set.seed(123)
rfe_profile <- rfe(x = x_train, y = y_train, sizes = sizes, rfeControl = control)
return(rfe_profile)
}
# --- Ejecución de RFE para cada modelo ---
# RFE para Velocidad
predictors_for_vel <- setdiff(all_candidate_predictors, c("pecen_70", "median_vel", "indice_exceso", "indice_riesgo"))
rfe_profile_vel <- run_rfe(df_train, predictors_for_vel, "vel_promed", subset_sizes, rfe_control)
best_predictors_vel <- predictors(rfe_profile_vel)
# RFE para Índice de Exceso
predictors_for_exc <- setdiff(all_candidate_predictors, c("vel_promed", "percen_70", "median_vel", "indice_riesgo"))
rfe_profile_exc <- run_rfe(df_train, predictors_for_exc, "indice_exceso", subset_sizes, rfe_control)
best_predictors_exc <- predictors(rfe_profile_exc)
# RFE para Índice de Riesgo (usando el mejor peso)
df_train_risk_final <- df_train %>% mutate(indice_riesgo_final = n_muer * BEST_RISK_WEIGHT + n_her)
rfe_profile_risk <- run_rfe(df_train_risk_final, all_candidate_predictors, "indice_riesgo_final", subset_sizes, rfe_control)
best_predictors_risk <- predictors(rfe_profile_risk)
# RFE para Siniestros
rfe_profile_log <- run_rfe(df_train, all_candidate_predictors, "siniestros_dummy", subset_sizes, rfe_control)
best_predictors_log <- predictors(rfe_profile_log)Se presentan los resultados.
# --- Presentación de resultados de RFE ---
print(rfe_profile_vel)##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (5 fold)
##
## Resampling performance over subset size:
##
## Variables RMSE Rsquared MAE RMSESD RsquaredSD MAESD Selected
## 5 9.181 0.4660 6.590 0.3012 0.02238 0.1686
## 8 8.952 0.4910 6.262 0.3003 0.01988 0.1688
## 10 8.963 0.4898 6.189 0.3480 0.02547 0.1547
## 12 8.927 0.4945 6.165 0.3010 0.01908 0.1630
## 15 8.878 0.5013 6.163 0.2918 0.01977 0.1492
## 20 8.865 0.5041 6.182 0.2777 0.01761 0.1375
## 24 8.858 0.5044 6.173 0.2857 0.01864 0.1555 *
##
## The top 5 variables (out of 24):
## longitud_k, rutas_sitp, dist_min_camara, mean_denem, mean_denpoplot(rfe_profile_vel, type = c("g", "o"), main = "Perfil RFE - Velocidad promedio")
cat("Predictores óptimos para la velocidad promedio:", paste(best_predictors_vel, collapse = ", "), "\n")## Predictores óptimos para la velocidad promedio: longitud_k, rutas_sitp, dist_min_camara, mean_denem, mean_denpo, mean_v_ref, mean_estra, ancho_carril, tipo, carriles, flujo_1, estaashto, n_arboles, tipomallae, estrato_fi, ciclorruta, n_semaforo, camara_500, sitp, camara_300, aashto_bue, camara_200, camara_100, pres_andenprint(rfe_profile_exc)##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (5 fold)
##
## Resampling performance over subset size:
##
## Variables RMSE Rsquared MAE RMSESD RsquaredSD MAESD Selected
## 5 0.1793 0.4856 0.1131 0.005860 0.01393 0.003123
## 8 0.1764 0.5014 0.1080 0.006859 0.02602 0.003487
## 10 0.1709 0.5311 0.1039 0.007301 0.02366 0.003721
## 12 0.1690 0.5418 0.1023 0.009823 0.03712 0.005242
## 15 0.1689 0.5435 0.1046 0.008422 0.03033 0.004057
## 20 0.1686 0.5457 0.1051 0.008875 0.03292 0.004245
## 24 0.1682 0.5475 0.1047 0.009227 0.03452 0.004449 *
##
## The top 5 variables (out of 24):
## longitud_k, dist_min_camara, rutas_sitp, mean_denem, mean_denpoplot(rfe_profile_exc, type = c("g", "o"), main = "Perfil RFE - Índice de exceso")
cat("Predictores óptimos para el índice de exceso:", paste(best_predictors_exc, collapse = ", "), "\n")## Predictores óptimos para el índice de exceso: longitud_k, dist_min_camara, rutas_sitp, mean_denem, mean_denpo, mean_estra, mean_v_ref, flujo_1, ancho_carril, carriles, tipo, tipomallae, n_arboles, estaashto, estrato_fi, ciclorruta, camara_500, n_semaforo, camara_300, sitp, camara_200, camara_100, aashto_bue, pres_andenprint(rfe_profile_risk)##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (5 fold)
##
## Resampling performance over subset size:
##
## Variables RMSE Rsquared MAE RMSESD RsquaredSD MAESD Selected
## 5 6.302 0.02893 3.977 0.2409 0.010611 0.08303
## 8 6.330 0.03258 4.023 0.2414 0.013054 0.07146
## 10 6.338 0.03214 4.044 0.2185 0.010641 0.07383
## 12 6.328 0.03151 4.048 0.2224 0.004981 0.06077
## 15 6.271 0.03610 4.036 0.2620 0.004783 0.05507
## 20 6.236 0.03744 4.025 0.2597 0.003835 0.07097
## 25 6.227 0.03822 4.027 0.2509 0.004106 0.07207
## 28 6.223 0.03707 4.030 0.2561 0.002260 0.06779 *
##
## The top 5 variables (out of 28):
## rutas_sitp, mean_v_ref, longitud_k, mean_estra, mean_denemplot(rfe_profile_risk, type = c("g", "o"), main = "Perfil RFE - Índice de riesgo (6:1)")
cat("Predictores óptimos para el índice de riesgo:", paste(best_predictors_risk, collapse = ", "), "\n")## Predictores óptimos para el índice de riesgo: rutas_sitp, mean_v_ref, longitud_k, mean_estra, mean_denem, mean_denpo, dist_min_camara, vel_promed, median_vel, tipomallae, percen_70, tipo, ancho_carril, ciclorruta, indice_exceso, estrato_fi, flujo_1, n_semaforo, n_arboles, carriles, camara_500, camara_200, camara_100, camara_300, aashto_bue, estaashto, sitp, pres_andencat("\n--- Resultados Finales de RFE ---\n")##
## --- Resultados Finales de RFE ---print(rfe_profile_log)##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (5 fold)
##
## Resampling performance over subset size:
##
## Variables Accuracy Kappa AccuracySD KappaSD Selected
## 5 0.5540 0.09253 0.010328 0.02372
## 8 0.5571 0.09854 0.019483 0.03870
## 10 0.5596 0.10273 0.009391 0.01824
## 12 0.5632 0.10735 0.011133 0.02294
## 15 0.5600 0.10029 0.007372 0.01481
## 20 0.5676 0.11558 0.011162 0.02315
## 25 0.5694 0.11871 0.012883 0.02593 *
## 28 0.5659 0.11057 0.010475 0.02134
##
## The top 5 variables (out of 25):
## rutas_sitp, mean_v_ref, carriles, longitud_k, mean_estraplot(rfe_profile_log, type = c("g", "o"), main = "Perfil RFE - Siniestros")
cat("Predictores óptimos para la probabilidad de siniestros:", paste(best_predictors_log, collapse = ", "), "\n")## Predictores óptimos para la probabilidad de siniestros: rutas_sitp, mean_v_ref, carriles, longitud_k, mean_estra, tipo, mean_denpo, dist_min_camara, median_vel, ancho_carril, percen_70, vel_promed, mean_denem, tipomallae, indice_exceso, flujo_1, n_arboles, estrato_fi, estaashto, aashto_bue, pres_anden, ciclorruta, camara_500, n_semaforo, camara_300Modelo Random Forest
Selección final de variables y construcción de fórmulas
Armados con los resultados del análisis de sensibilidad, ahora construimos la versión final y optimizada de nuestros cuatro modelos de Random Forest. Utilizamos los subconjuntos de predictores identificados por RFE y la ponderación de riesgo óptima para asegurar que cada modelo esté afinado para su tarea específica.
# Modelo de Velocidad: Seleccionamos los 20 mejores.
best_predictors_vel_final <- head(rfe_profile_vel$variables$var, 20)
cat("\nDecisión Velocidad: Seleccionando los 20 mejores predictores para parsimonia.\n")##
## Decisión Velocidad: Seleccionando los 20 mejores predictores para parsimonia.print(best_predictors_vel_final)## [1] "longitud_k" "rutas_sitp" "dist_min_camara" "mean_denem"
## [5] "mean_v_ref" "mean_estra" "mean_denpo" "ancho_carril"
## [9] "tipo" "carriles" "flujo_1" "estaashto"
## [13] "n_arboles" "tipomallae" "estrato_fi" "ciclorruta"
## [17] "n_semaforo" "camara_500" "sitp" "camara_300"# Modelo de Índice de Exceso: Seleccionamos los 20 mejores.
best_predictors_exc_final <- head(rfe_profile_exc$variables$var, 20)
cat("\nDecisión Exceso: Seleccionando los 20 mejores predictores para parsimonia.\n")##
## Decisión Exceso: Seleccionando los 20 mejores predictores para parsimonia.print(best_predictors_exc_final)## [1] "longitud_k" "dist_min_camara" "rutas_sitp" "mean_denem"
## [5] "mean_estra" "mean_v_ref" "mean_denpo" "flujo_1"
## [9] "tipo" "ancho_carril" "carriles" "tipomallae"
## [13] "n_arboles" "estaashto" "estrato_fi" "ciclorruta"
## [17] "camara_500" "camara_300" "n_semaforo" "camara_100"# Modelo de Índice de Riesgo: Usamos la recomendación completa de RFE (28).
best_predictors_risk_final <- predictors(rfe_profile_risk)
cat("\nDecisión Riesgo: Usando la selección completa de RFE de", length(best_predictors_risk_final), "predictores.\n")##
## Decisión Riesgo: Usando la selección completa de RFE de 28 predictores.print(best_predictors_risk_final)## [1] "rutas_sitp" "mean_v_ref" "longitud_k" "mean_estra"
## [5] "mean_denem" "mean_denpo" "dist_min_camara" "vel_promed"
## [9] "median_vel" "tipomallae" "percen_70" "tipo"
## [13] "ancho_carril" "ciclorruta" "indice_exceso" "estrato_fi"
## [17] "flujo_1" "n_semaforo" "n_arboles" "carriles"
## [21] "camara_500" "camara_200" "camara_100" "camara_300"
## [25] "aashto_bue" "estaashto" "sitp" "pres_anden"# Modelo de Siniestros: Seleccionamos los 20 mejores.
best_predictors_log_final <- head(rfe_profile_log$variables$var, 20)
cat("Decisión Siniestros: Seleccionando los 20 mejores predictores para parsimonia.\n")## Decisión Siniestros: Seleccionando los 20 mejores predictores para parsimonia.print(best_predictors_log_final)## [1] "rutas_sitp" "mean_v_ref" "longitud_k" "carriles"
## [5] "mean_estra" "ancho_carril" "median_vel" "percen_70"
## [9] "dist_min_camara" "vel_promed" "tipomallae" "indice_exceso"
## [13] "mean_denem" "mean_denpo" "tipo" "ciclorruta"
## [17] "flujo_1" "estrato_fi" "n_arboles" "pres_anden"# --- Construcción de las fórmulas finales ---
# --- Paso 1: Actualizar el Índice de Riesgo con la Ponderación Óptima ---
df_train <- df_train %>% mutate(indice_riesgo = n_muer * BEST_RISK_WEIGHT + n_her)
df_test <- df_test %>% mutate(indice_riesgo = n_muer * BEST_RISK_WEIGHT + n_her)
# --- Paso 2: Crear las Fórmulas Finales basadas en las selecciones parsimoniosas ---
form_vel_final <- reformulate(best_predictors_vel_final, response = "vel_promed")
form_exc_final <- reformulate(best_predictors_exc_final, response = "indice_exceso")
form_risk_final <- reformulate(best_predictors_risk_final, response = "indice_riesgo")
form_log_final <- reformulate(best_predictors_log_final, response = "siniestros_dummy")
# --- PASO 3 (CORRECCIÓN): Crear DataFrames de Entrenamiento Limpios para cada Modelo ---
# Este paso garantiza que no haya NAs en los datos que se pasan a train()
train_vel_final <- prep_df(df_train, form_vel_final)
train_exc_final <- prep_df(df_train, form_exc_final)
train_risk_final <- prep_df(df_train, form_risk_final)
train_log_final <- prep_df(df_train, form_log_final)Entrenamiento de los modelos
Ahora, entrenamos los cuatro modelos de Random Forest utilizando las fórmulas y conjuntos de datos optimizados. Implementamos una validación cruzada robusta y una búsqueda de hiperparámetros para asegurar que cada modelo esté afinado para su variable objetivo.
# --- Entrenamiento de los 4 modelos ---
set.seed(123)
# 1. Modelo de Velocidad (percen_70)
rf_vel_final <- train(form_vel_final, data = train_vel_final, method = "ranger", trControl = ctrl_reg,
preProcess = c("center","scale"), tuneGrid = tuneGrid_rf,
metric = "RMSE", importance = "permutation", num.trees = 1200)
rf_exc_final <- train(form_exc_final, data = train_exc_final, method = "ranger", trControl = ctrl_reg,
preProcess = c("center","scale"), tuneGrid = tuneGrid_rf,
metric = "RMSE", importance = "permutation", num.trees = 1200)
rf_risk_final <- train(form_risk_final, data = train_risk_final, method = "ranger", trControl = ctrl_reg,
preProcess = c("center","scale"), tuneGrid = tuneGrid_rf,
metric = "RMSE", importance = "permutation", num.trees = 1200)
rf_log_final <- train(form_log_final, data = train_log_final, method = "ranger", trControl = ctrl_clf,
preProcess = c("center","scale"), tuneGrid = tuneGrid_rf_clf,
metric = "ROC", importance = "permutation", num.trees = 1200)Resultados de los modelos:
# --- Imprimir resultados y gráficos de importancia ---
print("Resultados del Modelo de Velocidad:")## [1] "Resultados del Modelo de Velocidad:"print(rf_vel_final)## Random Forest
##
## 8000 samples
## 20 predictor
##
## Pre-processing: centered (33), scaled (33)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6399, 6400, 6400, 6401, 6400
## Resampling results across tuning parameters:
##
## mtry min.node.size RMSE Rsquared MAE
## 5 15 9.189760 0.4770460 6.644357
## 5 25 9.346304 0.4628437 6.834598
## 10 15 9.001416 0.4910951 6.394927
## 10 25 9.119354 0.4805338 6.563635
## 15 15 8.927965 0.4973509 6.299782
## 15 25 9.057711 0.4846918 6.481528
##
## 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 = 15, splitrule = variance
## and min.node.size = 15.plot(varImp(rf_vel_final), top = 20, main = "Importancia de Variables - Velocidad (P70)")
print("Resultados del Modelo de Índice de Exceso:")## [1] "Resultados del Modelo de Índice de Exceso:"print(rf_exc_final)## Random Forest
##
## 8000 samples
## 20 predictor
##
## Pre-processing: centered (33), scaled (33)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6400, 6400, 6400, 6400, 6400
## Resampling results across tuning parameters:
##
## mtry min.node.size RMSE Rsquared MAE
## 5 15 0.1753380 0.5173224 0.1149928
## 5 25 0.1788924 0.5006179 0.1189244
## 10 15 0.1726207 0.5256950 0.1102338
## 10 25 0.1753657 0.5127386 0.1136159
## 15 15 0.1717573 0.5288762 0.1085663
## 15 25 0.1743372 0.5160441 0.1117886
##
## 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 = 15, splitrule = variance
## and min.node.size = 15.plot(varImp(rf_exc_final), top = 20, main = "Importancia de Variables - Índice de Exceso")
print("Resultados del Modelo de Índice de Riesgo:")## [1] "Resultados del Modelo de Índice de Riesgo:"print(rf_risk_final)## Random Forest
##
## 8000 samples
## 28 predictor
##
## Pre-processing: centered (41), scaled (41)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6401, 6400, 6400, 6399, 6400
## Resampling results across tuning parameters:
##
## mtry min.node.size RMSE Rsquared MAE
## 5 15 6.104993 0.04358503 3.939708
## 5 25 6.085078 0.04584511 3.919782
## 10 15 6.144101 0.04138142 3.975705
## 10 25 6.113536 0.04368341 3.951104
## 15 15 6.163432 0.04105818 3.993410
## 15 25 6.133040 0.04265168 3.969235
##
## 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 = 5, splitrule = variance
## and min.node.size = 25.plot(varImp(rf_risk_final), top = 20, main = "Importancia de Variables - Índice de Riesgo")
print("Resultados del Modelo de Clasificación de Siniestros:")## [1] "Resultados del Modelo de Clasificación de Siniestros:"print(rf_log_final)## Random Forest
##
## 8000 samples
## 20 predictor
## 2 classes: 'Sí', 'No'
##
## Pre-processing: centered (28), scaled (28)
## Resampling: Cross-Validated (5 fold)
## Summary of sample sizes: 6401, 6399, 6401, 6400, 6399
## Resampling results across tuning parameters:
##
## mtry splitrule min.node.size ROC Sens Spec
## 5 gini 1 0.5749186 0.6616770 0.4595019
## 5 gini 5 0.5787524 0.6700706 0.4550505
## 5 gini 10 0.5833884 0.6830056 0.4441967
## 5 extratrees 1 0.5741978 0.6580455 0.4500424
## 5 extratrees 5 0.5809524 0.6759727 0.4422500
## 5 extratrees 10 0.5894586 0.7027497 0.4316778
## 10 gini 1 0.5733358 0.6532797 0.4634005
## 10 gini 5 0.5768731 0.6643983 0.4558862
## 10 gini 10 0.5810129 0.6764234 0.4508777
## 10 extratrees 1 0.5726531 0.6501021 0.4600602
## 10 extratrees 5 0.5781597 0.6673475 0.4461454
## 10 extratrees 10 0.5859679 0.6859557 0.4378020
## 15 gini 1 0.5733909 0.6489690 0.4675718
## 15 gini 5 0.5763323 0.6605411 0.4542145
## 15 gini 10 0.5794733 0.6741571 0.4517106
## 15 extratrees 1 0.5708081 0.6471524 0.4642319
## 15 extratrees 5 0.5781148 0.6637165 0.4558870
## 15 extratrees 10 0.5848640 0.6823253 0.4455902
##
## ROC was used to select the optimal model using the largest value.
## The final values used for the model were mtry = 5, splitrule = extratrees
## and min.node.size = 10.plot(varImp(rf_log_final), top = 20, main = "Importancia de Variables - Siniestros")
Evaluación de los modelos
Con los modelos optimizados ya entrenados, el paso final y más
crucial es evaluar su rendimiento en el conjunto de datos de prueba
(df_test). Este conjunto contiene información que los
modelos no han visto durante su entrenamiento, proporcionándonos una
medida honesta y realista de su capacidad para generalizar y hacer
predicciones en nuevos escenarios viales.
Para esta evaluación, primero creamos “sets de prueba limpios” para cada modelo, asegurando que solo utilizamos las filas que no tienen valores faltantes para los predictores específicos de ese modelo. Esto garantiza una comparación justa y sin errores.
# --- Paso 1: Crear sets de prueba limpios para cada modelo final ---
# Esto asegura que las predicciones y las observaciones tengan la misma longitud.
vars_vel_final <- all.vars(form_vel_final)
vars_exc_final <- all.vars(form_exc_final)
vars_risk_final <- all.vars(form_risk_final)
vars_log_final <- all.vars(form_log_final)
# Creamos un dataframe de prueba limpio y específico para cada modelo
df_test_vel_clean <- df_test %>% select(all_of(vars_vel_final)) %>% drop_na()
df_test_exc_clean <- df_test %>% select(all_of(vars_exc_final)) %>% drop_na()
df_test_risk_clean <- df_test %>% select(all_of(vars_risk_final)) %>% drop_na()
df_test_log_clean <- df_test %>% select(all_of(vars_log_final)) %>% drop_na()
# --- Paso 2: Generar Predicciones y Evaluar Modelos de Regresión ---
# Generamos las predicciones usando los modelos finales y los sets de prueba limpios
pred_vel <- predict(rf_vel_final, df_test_vel_clean)
pred_exc <- predict(rf_exc_final, df_test_exc_clean)
pred_risk <- predict(rf_risk_final, df_test_risk_clean)
# Calculamos las métricas de rendimiento (RMSE, Rsquared, MAE)
# CORRECCIÓN: La variable objetivo para el modelo de velocidad es 'percen_70'
resultados_vel <- postResample(pred = pred_vel, obs = df_test_vel_clean$percen_70)
resultados_exc <- postResample(pred = pred_exc, obs = df_test_exc_clean$indice_exceso)
resultados_risk <- postResample(pred = pred_risk, obs = df_test_risk_clean$indice_riesgo)
# --- Paso 3: Evaluar el Modelo de Clasificación ---
# Predicciones de clase ('Sí'/'No') y de probabilidad
pred_clase_rf <- predict(rf_log_final, df_test_log_clean)
pred_prob_rf <- predict(rf_log_final, df_test_log_clean, type = "prob")
# Matriz de Confusión para métricas como Accuracy, Sensitivity, etc.
cm_rf <- confusionMatrix(data = pred_clase_rf, reference = df_test_log_clean$siniestros_dummy, positive = "Sí")
cm_rf## Confusion Matrix and Statistics
##
## Reference
## Prediction Sí No
## Sí 1327 856
## No 561 684
##
## Accuracy : 0.5866
## 95% CI : (0.5699, 0.6032)
## No Information Rate : 0.5508
## P-Value [Acc > NIR] : 1.234e-05
##
## Kappa : 0.1497
##
## Mcnemar's Test P-Value : 5.710e-15
##
## Sensitivity : 0.7029
## Specificity : 0.4442
## Pos Pred Value : 0.6079
## Neg Pred Value : 0.5494
## Prevalence : 0.5508
## Detection Rate : 0.3871
## Detection Prevalence : 0.6368
## Balanced Accuracy : 0.5735
##
## 'Positive' Class : Sí
## # Curva ROC y cálculo del AUC (Área Bajo la Curva)
roc_rf <- roc(response = df_test_log_clean$siniestros_dummy, predictor = pred_prob_rf$Sí, levels = c("No", "Sí"))## Setting direction: controls < casesauc_rf <- auc(roc_rf)
# --- Paso 4: Consolidar y Presentar Resultados en una Tabla Final ---
# Creamos dataframes separados para los resultados de regresión y clasificación
tabla_resultados_reg <- data.frame(
Modelo = c("Velocidad (P70) - RF", "Índice de exceso - RF", "Índice de riesgo - RF"),
RMSE = c(resultados_vel["RMSE"], resultados_exc["RMSE"], resultados_risk["RMSE"]),
Rsquared = c(resultados_vel["Rsquared"], resultados_exc["Rsquared"], resultados_risk["Rsquared"]),
MAE = c(resultados_vel["MAE"], resultados_exc["MAE"], resultados_risk["MAE"])
)
tabla_resultados_clf <- data.frame(
Modelo = "Siniestros - RF",
Accuracy = cm_rf$overall["Accuracy"],
Sensitivity = cm_rf$byClass["Sensitivity"],
Specificity = cm_rf$byClass["Specificity"],
AUC = as.numeric(auc_rf)
)
# Unimos ambas tablas en una tabla final para una visión completa
tabla_final <- bind_rows(tabla_resultados_reg, tabla_resultados_clf) %>%
select(Modelo, RMSE, Rsquared, MAE, Accuracy, Sensitivity, Specificity, AUC)
# Presentamos la tabla con un formato profesional
kable(tabla_final, digits = 4, caption = "Rendimiento Final de Modelos Optimizados en Validación") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
column_spec(2:4, background = ifelse(is.na(tabla_final$RMSE), "#F8F8F8", "white")) %>%
column_spec(5:8, background = ifelse(is.na(tabla_final$Accuracy), "#F8F8F8", "white"))| Modelo | RMSE | Rsquared | MAE | Accuracy | Sensitivity | Specificity | AUC | |
|---|---|---|---|---|---|---|---|---|
| …1 | Velocidad (P70) - RF | 8.9434 | 0.4833 | 6.3434 | NA | NA | NA | NA |
| …2 | Índice de exceso - RF | 0.1731 | 0.5075 | 0.1079 | NA | NA | NA | NA |
| …3 | Índice de riesgo - RF | 5.7208 | 0.0624 | 3.7955 | NA | NA | NA | NA |
| Accuracy | Siniestros - RF | NA | NA | NA | 0.5866 | 0.7029 | 0.4442 | 0.601 |
Análisis de Dependencia Parcial (PDP)
Con los modelos finales entrenados, ahora exploramos cómo las variables predictoras más importantes influyen en los resultados. Los Gráficos de Dependencia Parcial (PDP) nos permiten aislar el efecto de una sola variable sobre la predicción del modelo, manteniendo todas las demás variables constantes en su valor promedio. Esto nos ayuda a entender las relaciones que el modelo ha aprendido de los datos, respondiendo a preguntas como: ¿cómo cambia la probabilidad de siniestro a medida que aumenta el número de rutas SITP?
Para este análisis, crearemos visualizaciones para cada uno de nuestros cuatro modelos, centrándonos en las variables que resultaron más relevantes en el análisis de RFE.
# --- PREPARACIÓN PARA LAS VISUALIZACIONES ---
# Creamos una grilla de predicción para enfocar el análisis de la distancia a la cámara
# en los primeros 500 metros, que es donde se espera el mayor impacto.
pg_0_500 <- data.frame(dist_min_camara = seq(0, 500, by = 10))
# FUNCIÓN AUXILIAR DE PLOTEO
# Esta función crea un ggplot estandarizado para cualquier PDP,
# manejando tanto variables numéricas (líneas) como categóricas (barras).
create_pdp_plot <- function(pdp_data, title, x_lab, y_lab, color = "#0072B2") {
# Extraemos el nombre de la variable predictora de la primera columna
pred_var_name <- names(pdp_data)[1]
p <- ggplot(pdp_data, aes(x = .data[[pred_var_name]], y = yhat))
# Si la variable es un factor (categórica), usamos barras
if (is.factor(pdp_data[[pred_var_name]])) {
p <- p + geom_col(fill = color, width = 0.7)
} else { # Si es numérica, usamos una línea
p <- p + geom_line(color = color, linewidth = 1.2)
}
p <- p + labs(
title = title,
x = x_lab,
y = y_lab
) +
theme_minimal(base_size = 14)
return(p)
}Visualización de PDP para el modelo de velocidad (P70)
Estos gráficos muestran cómo cambia la velocidad (percentil 70) estimada al variar diferentes características de la vía y su entorno.
# --- Cálculo de PDPs para el Modelo de Velocidad ---
pdp_vel_denpo <- partial(rf_vel_final, pred.var = "mean_denpo")
pdp_vel_dist <- partial(rf_vel_final, pred.var = "dist_min_camara", pred.grid = pg_0_500)
pdp_vel_denem <- partial(rf_vel_final, pred.var = "mean_denem")
pdp_vel_anc_c <- partial(rf_vel_final, pred.var = "ancho_carril")
pdp_vel_estaa <- partial(rf_vel_final, pred.var = "estaashto")
pdp_vel_carril <- partial(rf_vel_final, pred.var = "carriles")
pdp_vel_sitp <- partial(rf_vel_final, pred.var = "rutas_sitp")
pdp_vel_semaf <- partial(rf_vel_final, pred.var = "n_semaforo")
pdp_vel_estra <- partial(rf_vel_final, pred.var = "mean_estra")
pdp_vel_vref <- partial(rf_vel_final, pred.var = "mean_v_ref")
pdp_vel_arbol <- partial(rf_vel_final, pred.var = "n_arboles")# --- Generación de Gráficos ---
plot_color_vel <- "darkgreen"
# Gráfico 1: Longitud del Corredor
create_pdp_plot(pdp_vel_denpo, "Efecto parcial de la densidad poblacional", "habs/ha", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 2: Distancia a Cámara
create_pdp_plot(pdp_vel_dist, "Efecto parcial de la distancia a cámara", "Distancia a cámara más cercana (m)", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 3: Densidad de Empleos
create_pdp_plot(pdp_vel_denem, "Efecto parcial de la densidad de empleos", "Densidad de empleos (empleos/ha)", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 4: Estado de la Vía (AASHTO)
create_pdp_plot(pdp_vel_estaa, "Efecto parcial del estado de la vía", "Estado de la vía (AASHTO)", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 5: Número de Carriles
create_pdp_plot(pdp_vel_carril, "Efecto parcial del número de carriles", "Número de carriles", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 6: Rutas SITP
create_pdp_plot(pdp_vel_sitp, "Efecto parcial del número de rutas SITP", "Número de rutas SITP en el tramo", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 7: Semáforos
create_pdp_plot(pdp_vel_semaf, "Efecto parcial del número de semáforos", "Número de semáforos en el tramo", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 8: Estratificación Socioeconómica
create_pdp_plot(pdp_vel_estra, "Efecto parcial de la estratificación socioeconómica", "Estrato socioeconómico", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
# Gráfico 9: Valor catastral promedio
create_pdp_plot(pdp_vel_vref, "Efecto parcial del valor catastral promedio", "Valor catastral promedio (Miles de millones - COP)", "Velocidad promedio estimada (km/h)", color = plot_color_vel) + scale_x_continuous(labels = scales::dollar_format(scale = 1e-6, suffix = "M"))
# Gráfico 10: Árboles
create_pdp_plot(pdp_vel_arbol, "Efecto parcial del número de árboles", "Número de árboles en el tramo", "Velocidad promedio estimada (km/h)", color = plot_color_vel)
Visualización de PDP para el modelo de índice de exceso
Estos gráficos muestran cómo cambia la proporción de tiempo estimada que se excede el límite de velocidad al variar los predictores más importantes.
# --- Cálculo de PDPs para el Modelo de Exceso ---
pdp_exc_denpo <- partial(rf_exc_final, pred.var = "mean_denpo")
pdp_exc_dist <- partial(rf_exc_final, pred.var = "dist_min_camara", pred.grid = pg_0_500)
pdp_exc_rutas <- partial(rf_exc_final, pred.var = "rutas_sitp")
pdp_exc_denem <- partial(rf_exc_final, pred.var = "mean_denem")
pdp_exc_anc_c <- partial(rf_exc_final, pred.var = "ancho_carril")
pdp_exc_vref <- partial(rf_exc_final, pred.var = "mean_v_ref")
pdp_exc_estra <- partial(rf_exc_final, pred.var = "mean_estra")
pdp_exc_sitp <- partial(rf_exc_final, pred.var = "rutas_sitp")
pdp_exc_carril<- partial(rf_exc_final, pred.var = "carriles")
pdp_exc_estaa <- partial(rf_exc_final, pred.var = "estaashto")
pdp_exc_arbol <- partial(rf_exc_final, pred.var = "n_arboles")
pdp_exc_semaf <- partial(rf_exc_final, pred.var = "n_semaforo")
pdp_exc_cam500<- partial(rf_exc_final, pred.var = "camara_500")# --- Generación de Gráficos ---
plot_color_exc <- "#D55E00"
# Gráfico 1: Densidad poblacional
create_pdp_plot(pdp_exc_denpo, "Efecto parcial de la densidad poblacional", "Densidad poblacional (habs/ha)", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 2: Distancia a Cámara
create_pdp_plot(pdp_exc_dist, "Efecto parcial de la distancia a la cámara", "Distancia a la cámara más cercana (m)", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 3: Rutas SITP
create_pdp_plot(pdp_exc_rutas, "Efecto parcial de rutas SITP", "Número de rutas SITP en el tramo", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 4: Ancho de Carril
create_pdp_plot(pdp_exc_anc_c, "Efecto parcial del ancho de carril", "Ancho promedio de Carril (m)", "Índice de Exceso Estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 5: Valor Catastral
create_pdp_plot(pdp_exc_vref, "Efecto parcial del valor catastral promedio", "Valor catastral promedio (Miles de millones - COP)", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent) + scale_x_continuous(labels = scales::dollar_format(scale = 1e-6, suffix = "M"))
# Gráfico 6: Estratificación Socioeconómica
create_pdp_plot(pdp_exc_estra, "Efecto parcial de la estratificación socioeconómica", "Estrato socioeconómico", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 7: Número de Carriles
create_pdp_plot(pdp_exc_carril, "Efecto parcial del número de carriles", "Número de carriles", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 8: Estado de la Vía (AASHTO)
create_pdp_plot(pdp_exc_estaa, "Efecto parcial del estado de la vía", "Estado de la vía (AASHTO)", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 9: Árboles
create_pdp_plot(pdp_exc_arbol, "Efecto parcial del número de árboles", "Número de árboles en el tramo", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 10: Semáforos
create_pdp_plot(pdp_exc_semaf, "Efecto parcial del número de semáforos", "Número de semáforos en el tramo", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
# Gráfico 11: Cámaras (500m)
create_pdp_plot(pdp_exc_cam500, "Efecto parcial de cámaras (500m)", "Existencia de una cámara en un radio de 500m", "Índice de exceso estimado", color = plot_color_exc) + scale_y_continuous(labels = scales::percent)
Visualización de PDP para el modelo de índice de riesgo
Estos gráficos muestran cómo cambia el Índice de Riesgo estimado (usando la ponderación 6:1) al variar los predictores más importantes.
# --- Cálculo de PDPs para el Modelo de Riesgo ---
pdp_risk_velp <- partial(rf_risk_final, pred.var = "vel_promed")
pdp_risk_medv <- partial(rf_risk_final, pred.var = "median_vel")
pdp_risk_p70 <- partial(rf_risk_final, pred.var = "percen_70")
pdp_risk_dist <- partial(rf_risk_final, pred.var = "dist_min_camara", pred.grid = pg_0_500)
pdp_risk_vref <- partial(rf_risk_final, pred.var = "mean_v_ref")
pdp_risk_denpo <- partial(rf_risk_final, pred.var = "mean_denpo")
pdp_risk_denem <- partial(rf_risk_final, pred.var = "mean_denem")
pdp_risk_estra <- partial(rf_risk_final, pred.var = "mean_estra")
pdp_risk_rutas <- partial(rf_risk_final, pred.var = "rutas_sitp")
pdp_risk_carril<- partial(rf_risk_final, pred.var = "carriles")
pdp_risk_anc_c <- partial(rf_risk_final, pred.var = "ancho_carril")
pdp_risk_estaa <- partial(rf_risk_final, pred.var = "estaashto")
pdp_risk_semaf <- partial(rf_risk_final, pred.var = "n_semaforo")
pdp_risk_arbol <- partial(rf_risk_final, pred.var = "n_arboles")
pdp_risk_exc <- partial(rf_risk_final, pred.var = "indice_exceso")
pdp_risk_cam200 <- partial(rf_risk_final, pred.var = "camara_200")# --- Generación de Gráficos ---
plot_color_risk <- "darkred"
# Gráfico 1: Velocidad Promedio
create_pdp_plot(pdp_risk_velp, "Efecto parcial de la velocidad promedio sobre el índice de riesgo", "Velocidad promedio (km/h)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 2: Velocidad Mediana
create_pdp_plot(pdp_risk_medv, "Efecto parcial de la velocidad mediana sobre el índice de riesgo", "Velocidad mediana (km/h)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 3: Percentil 70 de Velocidad
create_pdp_plot(pdp_risk_p70, "Efecto parcial de la velocidad (P70) sobre el riesgo", "Velocidad (P70) (km/h)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 4: Distancia a Cámara
create_pdp_plot(pdp_risk_dist, "Efecto parcial de la distancia a la cámara más cercana sobre el índice de riesgo", "Distancia a la cámara más cercana (m)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 5: Valor Catastral
create_pdp_plot(pdp_risk_vref, "Efecto parcial del valor catastral promedio sobre el índice de riesgo", "Valor catastral promedio (Miles de millones - COP)", "Índice de riesgo estimado", color = plot_color_risk) + scale_x_continuous(labels = scales::dollar_format(scale = 1e-6, suffix = "M"))
# Gráfico 6: Densidad Poblacional
create_pdp_plot(pdp_risk_denpo, "Efecto parcial de la densidad poblacional sobre el índice de riesgo", "Densidad poblacional (habs/ha)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 7: Densidad de Empleos
create_pdp_plot(pdp_risk_denem, "Efecto parcial de la densidad de empleos sobre el índice de riesgo", "Densidad de empleos (empleos/ha)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 8: Estratificación Socioeconómica
create_pdp_plot(pdp_risk_estra, "Efecto parcial de la estratificación socioeconómica sobre el índice de riesgo", "Estrato socioeconómico", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 9: Rutas SITP
create_pdp_plot(pdp_risk_rutas, "Efecto parcial de rutas SITP sobre el índice de riesgo", "Número de rutas SITP en el tramo", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 10: Número de Carriles
create_pdp_plot(pdp_risk_carril, "Efecto parcial del número de carriles sobre el índice de riesgo", "Número de carriles", "Índice de riesgo estimado", color = plot_color_risk) +
scale_x_continuous(breaks = scales::pretty_breaks(n = 6))
# Gráfico 11: Ancho de Carril
create_pdp_plot(pdp_risk_anc_c, "Efecto parcial del ancho de carril sobre el índice de riesgo", "Ancho promedio de Carril (m)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 12: Estado de la Vía (AASHTO)
create_pdp_plot(pdp_risk_estaa, "Efecto parcial del estado de la vía sobre el índice de riesgo", "Estado de la vía (AASHTO)", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 13: Semáforos
create_pdp_plot(pdp_risk_semaf, "Efecto parcial del número de semáforos sobre el índice de riesgo", "Número de semáforos en el tramo", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 14: Árboles
create_pdp_plot(pdp_risk_arbol, "Efecto parcial del número de árboles sobre el índice de riesgo", "Número de árboles en el tramo", "Índice de riesgo estimado", color = plot_color_risk)
# Gráfico 15: Índice de Exceso
create_pdp_plot(pdp_risk_exc, "Efecto parcial del índice de exceso sobre el índice de riesgo", "Índice de exceso", "Índice de riesgo estimado", color = plot_color_risk) +
scale_x_continuous(labels = scales::percent)
# Gráfico 16: Cámaras (200m)
create_pdp_plot(pdp_risk_cam200, "Efecto parcial de cámaras (200m) sobre el índice de riesgo", "Existencia de una cámara en un radio de 200m", "Índice de riesgo estimado", color = plot_color_risk)
Visualización de PDP para el modelo de siniestros
Estos gráficos muestran cómo cambia la probabilidad estimada de que ocurra un siniestro al variar cada predictor.
# --- Cálculo de PDPs para el Modelo de Siniestros ---
# Usamos prob=TRUE y which.class="Sí" para obtener la probabilidad de la clase positiva.
pdp_log_velp <- partial(rf_log_final, pred.var = "vel_promed", prob = TRUE, which.class = "Sí")
pdp_log_medv <- partial(rf_log_final, pred.var = "median_vel", prob = TRUE, which.class = "Sí")
pdp_log_p70 <- partial(rf_log_final, pred.var = "percen_70", prob = TRUE, which.class = "Sí")
pdp_log_dist <- partial(rf_log_final, pred.var = "dist_min_camara", pred.grid = pg_0_500, prob = TRUE, which.class = "Sí")
pdp_log_vref <- partial(rf_log_final, pred.var = "mean_v_ref", prob = TRUE, which.class = "Sí")
pdp_log_denpo <- partial(rf_log_final, pred.var = "mean_denpo", prob = TRUE, which.class = "Sí")
pdp_log_denem <- partial(rf_log_final, pred.var = "mean_denem", prob = TRUE, which.class = "Sí")
pdp_log_estra <- partial(rf_log_final, pred.var = "mean_estra", prob = TRUE, which.class = "Sí")
pdp_log_rutas <- partial(rf_log_final, pred.var = "rutas_sitp", prob = TRUE, which.class = "Sí")
pdp_log_carril <- partial(rf_log_final, pred.var = "carriles", prob = TRUE, which.class = "Sí")
pdp_log_anc_c <- partial(rf_log_final, pred.var = "ancho_carril", prob = TRUE, which.class = "Sí")
pdp_log_arbol <- partial(rf_log_final, pred.var = "n_arboles", prob = TRUE, which.class = "Sí")
pdp_log_exc <- partial(rf_log_final, pred.var = "indice_exceso", prob = TRUE, which.class = "Sí")
pdp_log_flujo <- partial(rf_log_final, pred.var = "flujo_1", prob = TRUE, which.class = "Sí")# --- Generación de Gráficos ---
plot_color_log <- "#0072B2" # Azul
# Gráfico 1: Velocidad Promedio
create_pdp_plot(pdp_log_velp, "Efecto parcial de la velocidad promedio sobre la probabilidad de siniestro", "Velocidad promedio (km/h)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 2: Rutas SITP
create_pdp_plot(pdp_log_rutas, "Efecto parcial de rutas SITP sobre la probabilidad de siniestro", "Número de rutas SITP en el tramo", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 3: Valor Catastral
create_pdp_plot(pdp_log_vref, "Efecto parcial del valor catastral promedio sobre la probabilidad de siniestro", "Valor catastral promedio (Miles de millones - COP)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent) + scale_x_continuous(labels = scales::dollar_format(scale = 1e-6, suffix = "M"))
# Gráfico 4: Estratificación Socioeconómica
create_pdp_plot(pdp_log_estra, "Efecto parcial de la estratificación socioeconómica sobre la probabilidad de siniestro", "Estrato socioeconómico", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 5: Densidad Poblacional
create_pdp_plot(pdp_log_denpo, "Efecto parcial de la densidad poblacional sobre la probabilidad de siniestro", "Densidad poblacional (habs/ha)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 6: Número de Carriles
create_pdp_plot(pdp_log_carril, "Efecto parcial del número de carriles sobre la probabilidad de siniestro", "Número de carriles", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent) + scale_x_continuous(breaks = scales::pretty_breaks(n = 6))
# Gráfico 7: Índice de Exceso
create_pdp_plot(pdp_log_exc, "Efecto parcial del índice de exceso sobre la probabilidad de siniestro", "Índice de exceso (Proporción 24h)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent) + scale_x_continuous(labels = scales::percent)
# Gráfico 8: Flujo Vehicular
create_pdp_plot(pdp_log_flujo, "Efecto parcial del flujo vehicular sobre la probabilidad de siniestro", "Nivel de Flujo Vehicular", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 9: Distancia a Cámara
create_pdp_plot(pdp_log_dist, "Efecto parcial de la distancia a la cámara más cercana sobre la probabilidad de siniestro", "Distancia a la cámara más cercana (m)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 10: Velocidad (P70)
create_pdp_plot(pdp_log_p70, "Efecto parcial de la velocidad (P70) sobre la probabilidad de siniestro", "Velocidad (P70) (km/h)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 11: Velocidad Mediana
create_pdp_plot(pdp_log_medv, "Efecto parcial de la velocidad mediana sobre la probabilidad de siniestro", "Velocidad mediana (km/h)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 12: Ancho de Carril
create_pdp_plot(pdp_log_anc_c, "Efecto parcial del ancho de carril sobre la probabilidad de siniestro", "Ancho promedio de Carril (m)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 13: Árboles
create_pdp_plot(pdp_log_arbol, "Efecto parcial del número de árboles sobre la probabilidad de siniestro", "Número de árboles en el tramo", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
# Gráfico 14: Densidad de Empleos
create_pdp_plot(pdp_log_denem, "Efecto parcial de la densidad de empleos sobre la probabilidad de siniestro", "Densidad de empleos (empleos/ha)", "Probabilidad de siniestro estimada", color = plot_color_log) + scale_y_continuous(labels = scales::percent)
## Promedios ponderados para escenario base
Vamos a calcular el promedio ponderado (sum(valor_variable *
longitudtr) / longitud_total_red) de las cuatro variables que estamos
prediciendo para luego comparar con los escenarios que se van a modelar.
La longitud total de la red es la suma de la longitud de todos los
tramos en vias.
# --- Cálculo de promedios ponderados para el escenario base ---
longitud_total_red <- sum(vias$longitudtr)
promedio_ponderado <- function(var) {
sum(vias[[var]] * vias$longitudtr) / longitud_total_red
}
base_vel_prom <- promedio_ponderado("vel_promed")
base_ind_exc <- promedio_ponderado("indice_exceso")
base_ind_risk <- promedio_ponderado("indice_riesgo")
base_prob_sin <- promedio_ponderado("siniestros_dummy_num")
# Mostramos los resultados
data.frame(
Variable = c("Velocidad promedio", "Índice de exceso", "Índice de riesgo", "Probabilidad de siniestro"),
Promedio_Ponderado = c(base_vel_prom, base_ind_exc, base_ind_risk, base_prob_sin)
) %>%
kable(digits = 4, caption = "Promedios Ponderados Escenario Base") %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)| Variable | Promedio_Ponderado |
|---|---|
| Velocidad promedio | 37.2046 |
| Índice de exceso | 0.1797 |
| Índice de riesgo | 7.7100 |
| Probabilidad de siniestro | 0.6095 |
Guardar el entorno de trabajo
# --- Detener el clúster paralelo ---
stopCluster(cl)
save.image("modelos_entrenados.RData")