Análisis Estadístico de la Siniestralidad Vial en los Departamentos de Colombia (2015-2019)
Laura Sofía Vásquez Lozano Y Ana Ximena Rosas Medina
2025-05-12
Librerías necesarias
library(dplyr) # Para manipulación de datos
library(ggplot2) # Para visualización de datos
library(readxl) # Para leer archivos Excel
library(plyr) #usar funciones como count
library(knitr) #Para una mejor impresion
library(psych)
library(gridExtra)
library(corrplot) # Para visualizar la matriz de correlaciones
library(FactoMineR)
library(factoextra)
library(GPArotation) # para rotaciones
library(Metrics)SECTORES CRITICOS DE SINIESTRALIDAD VIAL (2015-2019)
Carga del archivo:
datos <- read_excel("BaseDatos_PIF.xlsx")
# Resumen detallado de los datos
describe(datos)## vars n mean sd median trimmed mad min max range
## __id* 1 316 158.50 91.37 158.50 158.50 117.13 1.00 316.00 315.00
## gipvalue 2 316 0.02 0.03 0.01 0.02 0.01 0.00 0.10 0.10
## divipola* 3 316 72.93 37.40 79.00 74.17 48.93 1.00 133.00 132.00
## id_mt* 4 316 6.87 8.21 1.00 5.48 0.00 1.00 25.00 24.00
## longitud 5 316 -75.17 1.44 -75.43 -75.18 1.62 -78.79 -71.72 7.06
## gizscore 6 316 3.19 1.64 2.70 2.90 1.08 1.66 14.54 12.88
## entidad* 7 316 2.50 1.39 3.00 2.50 1.48 1.00 4.00 3.00
## departamento* 8 316 12.24 7.05 10.50 12.42 8.15 1.00 22.00 21.00
## municipio* 9 316 64.29 39.81 59.00 63.36 53.37 1.00 132.00 131.00
## fallecidos 10 316 11.33 7.07 10.00 10.66 4.45 1.00 52.00 51.00
## tramo* 11 293 39.52 25.03 35.00 38.79 31.13 1.00 85.00 84.00
## nombre* 12 141 12.17 6.33 11.00 12.21 7.41 1.00 24.00 23.00
## latitud 13 316 5.67 2.75 4.75 5.54 2.37 0.82 11.52 10.70
## pr* 14 316 58.23 33.28 59.50 59.11 43.74 1.00 105.00 104.00
## skew kurtosis se
## __id* 0.00 -1.21 5.14
## gipvalue 1.13 -0.02 0.00
## divipola* -0.24 -1.17 2.10
## id_mt* 1.08 -0.35 0.46
## longitud 0.05 -0.47 0.08
## gizscore 2.49 9.93 0.09
## entidad* -0.06 -1.87 0.08
## departamento* -0.01 -1.20 0.40
## municipio* 0.15 -1.28 2.24
## fallecidos 1.46 4.46 0.40
## tramo* 0.21 -1.18 1.46
## nombre* -0.01 -1.14 0.53
## latitud 0.45 -0.77 0.15
## pr* -0.09 -1.26 1.871. Desbalanceo e Imputación de Datos
En esta primera etapa, realizamos un análisis preliminar para identificar dos aspectos clave del conjunto de datos: el desbalanceo entre clases y la presencia de valores faltantes. El desbalanceo se refiere a una distribución desigual en las categorías de la variable dependiente, lo cual puede afectar el rendimiento de los modelos predictivos.
#Ver cuántos NA hay en cada columna
colSums(is.na(datos))## __id gipvalue divipola id_mt longitud gizscore
## 0 0 0 0 0 0
## entidad departamento municipio fallecidos tramo nombre
## 0 0 0 0 23 175
## latitud pr
## 0 0Vemos que tramo tiene 23 valores faltantes y nombre tiene 175 valores faltantes, lo cual es bastante relevante. En cambio, las demás columnas no tienen valores faltantes. Ahora, pasemos al siguiente paso de imputar los valores faltantes. Para eso, tenemos varias opciones:
°tramo: Si esta variable es categórica, podríamos imputarla con la moda (el valor más frecuente).
°nombre: Este campo parece contener identificadores o categorías, por lo que también se podría imputar con la moda.
# Función para obtener la moda
getmode <- function(x) {
ux <- unique(x[!is.na(x)])
ux[which.max(tabulate(match(x, ux)))]
}
# Imputar NA en 'tramo' con la moda
moda_tramo <- getmode(datos$tramo)
datos$tramo[is.na(datos$tramo)] <- moda_tramo
# Imputar NA en 'nombre' con la moda
moda_nombre <- getmode(datos$nombre)
datos$nombre[is.na(datos$nombre)] <- moda_nombre
# Función para obtener la moda
getmode <- function(x) {
ux <- unique(x[!is.na(x)])
ux[which.max(tabulate(match(x, ux)))]
}
# Imputación para columna 'pr' (filtrando primero los "sd")
# Paso 1: Reemplazar valores "sd" por NA para que no afecten la moda
datos$pr[datos$pr == "sd"] <- NA
# Paso 2: Calcular la moda sin incluir NA ni "sd"
moda_pr <- getmode(datos$pr)
# Paso 3: Imputar los NA (antes "sd" o vacíos) con la moda correcta
datos$pr[is.na(datos$pr)] <- moda_pr
# Verificamos que ya no hay "sd" y todos los NA han sido imputados
table(datos$pr, useNA="always")##
## 0 1 10 102 103 104 106 107 109 11 111 114 117 118 12 126
## 8 6 5 2 1 1 1 2 1 1 1 1 2 1 7 1
## 128 129 13 14 15 16 17 19 2 21 22 23 24 25 26 27
## 1 1 6 3 6 4 1 4 5 1 1 3 2 4 3 3
## 28 29 3 30 31 32 33 34 35 36 37 38 4 40 41 42
## 2 3 6 2 1 1 1 4 2 2 4 4 2 2 2 2
## 43 45 46 47 48 49 5 50 51 52 53 54 55 56 57 58
## 2 2 1 2 4 3 4 3 2 4 2 1 3 1 2 6
## 59 6 60 61 62 63 64 65 67 68 69 7 70 71 72 73
## 2 9 3 3 5 1 1 1 4 2 3 1 3 3 4 1
## 74 75 76 77 78 79 8 80 81 82 84 85 86 87 88 89
## 2 2 1 1 4 2 51 1 1 3 2 3 2 3 1 1
## 9 91 92 93 94 96 97 98 <NA>
## 6 2 2 2 1 2 2 2 0Y para ver cómo quedaron los valores imputados, vamos a hacer un pequeño resumen de las columnas afectadas:
# Ver cuántas veces aparece la moda ahora
table(datos$tramo)##
## Aguachica - San Roque
## 6
## Arboletes-Montería
## 2
## Barbosa - Bucaramanga
## 5
## Barranca - Aguachica
## 9
## Barranquilla - Cartagena
## 2
## Barranquilla - Malambo - Palmar de Varela
## 4
## Bogotá - Barbosa
## 1
## Bogotá - Mosquera
## 2
## Bogotá - Tunja
## 4
## Bogotá - Villavicencio
## 4
## Bogotá - Villeta
## 3
## Bosconia - Cienaga
## 7
## Buenaventura-Citronela
## 4
## Buga- La Victoria
## 13
## Calarcá - Ibagué
## 4
## Cali-Dagua
## 6
## Cali - Jamundí - La Primavera
## 5
## Cali - Palmira
## 2
## Campoalegre-Neiva
## 3
## Candelaria-Buga
## 6
## Carmen de Bolívar-Bosconia
## 1
## Carreto-Cartagena
## 1
## Carreto - Barranquilla
## 2
## Cartagena-Barranquilla (Cordialidad)
## 1
## Caucasia-Planeta Rica
## 1
## Cereté - San Pelayo
## 1
## Chigorodó-Turbo
## 8
## Cienaga - Barranquilla
## 4
## Duitama-Sogamoso
## 2
## El Cerrito - Obando - Cali
## 1
## El Hatillo - Girardota - Copacabana - Bello - Med*
## 38
## El Zulia-Cúcuta
## 1
## Florencia - Guayabal - Altamira
## 4
## Florencia - San Vicente del Caguán
## 1
## Girardot - Bogotá
## 14
## Granada - Villavicencio
## 6
## Guaduas-Puerto Salgar
## 1
## Ibagué- Girardot
## 2
## Ibagué - Venadillo - Lérida
## 1
## Ipiales-Pedregal
## 2
## Ipiales - Aldana - Guachucal
## 2
## La Gloria- El Tarra-Campo Seis
## 2
## La Paila - Calarcá
## 1
## La Paz - Maicao
## 1
## La Victoria-La Virginia
## 1
## Loboguerrero - Buga
## 2
## Lorica - Chinú
## 1
## Medellín - Santuario
## 9
## Montería-La Ye
## 4
## Mosquera - Facatativa - Los Andes
## 3
## Mosquera - Funza
## 4
## Neiva-Girardot
## 5
## Palmira - El Palo
## 5
## Pasto-Popayán
## 4
## Pasto - La Unión - Mercaderes
## 2
## Pasto - Sandona - Yacuanquer
## 2
## Piedecuesta-Girón
## 1
## Planeta Rica - Montería
## 2
## Planeta Rica - Sincelejo
## 2
## Popayán-Santander de Quilichao
## 12
## Puente Rumichaca-Ipiales
## 3
## Puerta del Hierro-Magangué
## 1
## Puerto Araújo-Cimitarra
## 1
## Puerto Salgar - Barranca
## 2
## San Alberto - Bucaramanga
## 1
## San Juan de Arama-Granada
## 1
## San Pelayo-Lorica
## 2
## San Roque-La Paz
## 2
## San Roque - Bosconia
## 2
## Santa Marta-Maicao-Paraguachón
## 7
## Santafé de Antioquía - Chigorodó
## 2
## Santander de Quilichao-Candelaria
## 4
## Sincelejo - Toluviejo
## 2
## Sincerin - Arjona - Turbaco - Cartagena
## 3
## Tame - San Lope
## 1
## Tierralta-Puerto Libertador-Montelíbano
## 1
## Tumaco-Pedregal
## 7
## Tunja - Duitama
## 2
## Variante Pereira
## 2
## Villavicencio-Puente Arimena
## 5
## Villavicencio - Yopal
## 1
## Villeta-Honda
## 4
## Viterbo-Ansermanuevo
## 2
## Yopal - Tame
## 2
## Yumbo - Yotoco - Buga
## 7table(datos$nombre)##
## Armenia -Pereira - Manizales
## 2
## Autopistas al Mar 2
## 10
## Bogotá – La Vega – Villeta (BVLPV)
## 3
## Bogotá (Fontibón) – Faca – Los Alpes
## 5
## Briceño - Tunja - Sogamoso
## 8
## Cartagena – Barranquilla – “Circunvalar de la Prosperidad”
## 2
## Córdoba – Sucre
## 8
## Desarrollo Vial del Oriente de Medellín – DEVIMED
## 9
## Girardot – Ibagué – Cajamarca
## 3
## IP – Antioquia – Bolívar
## 8
## IP – Ampliación a Tercer Carril Doble Calzada Bogotá - Girardot
## 189
## IP – Cambao Manizales
## 1
## IP – Chirajara - Fundadores
## 4
## IP – Neiva – Espinal - Girardot
## 6
## IP Malla Vial del Meta
## 9
## Popayán – Santander de Quilichao
## 12
## Puerta de Hierro – Palmar de Varela y Carreto – Cruz del Viso
## 2
## Rumichaca – Pasto
## 5
## Ruta Caribe
## 9
## Ruta del Sol 1
## 1
## Ruta del Sol 3
## 10
## Santa Marta - Paraguachón
## 7
## Santana - Mocoa - Neiva
## 2
## Villavicencio - Yopal
## 1Verificar desbalance entre “fatales” y “no fatales”
Primero necesitamos una variable categórica que indique si un siniestro fue fatal (por ejemplo, si hay al menos 1 fallecido).
datos$fatal <- ifelse(datos$fallecidos > 0, "Fatal", "No Fatal")
table(datos$fatal)##
## Fatal
## 316prop.table(table(datos$fatal))##
## Fatal
## 1No podemos hacer sobremuestreo ni submuestreo porque no hay dos clases. Solo se tiene eventos con fallecidos. Entonces, esta parte la saltamos porque no aplica en este caso.
2. Valores atípicos (outliers) en las variables numéricas:
Vamos a detectar y tratar valores atípicos (outliers) en las variables numéricas. Los outliers son valores extremos que pueden afectar el análisis estadístico y la visualización de datos, especialmente en medidas como la media, la desviación estándar y en modelos como regresiones.
¿Cómo lo haremos?
Usaremos boxplots (diagramas de caja) para visualizar outliers. También usaremos el criterio del rango intercuartílico (IQR) para detectarlos numéricamente.
#Crear boxplots para variables numéricas
# Boxplot para la variable fallecidos
ggplot(datos, aes(y = fallecidos)) +
geom_boxplot(fill = "#FF9999") +
labs(title = "Boxplot de Fallecidos", y = "Número de Fallecidos")
# Boxplot para gizscore
ggplot(datos, aes(y = gizscore)) +
geom_boxplot(fill = "#99CCFF") +
labs(title = "Boxplot de GIZScore", y = "GIZScore")
# Boxplot para longitud
ggplot(datos, aes(y = longitud)) +
geom_boxplot(fill = "#CCFFCC") +
labs(title = "Boxplot de Longitud", y = "Longitud")
# Boxplot para latitud
ggplot(datos, aes(y = latitud)) +
geom_boxplot(fill = "#FFFFCC") +
labs(title = "Boxplot de Latitud", y = "Latitud")
Análisis de los Boxplots:
En el boxplot de fallecidos, la mayoría de los valores están concentrados en torno a la mediana (aproximadamente 10), con una distribución bastante concentrada. Sin embargo, se observan varios outliers (puntos fuera de los bigotes). Esto indica que hay algunos casos con un número de fallecidos mucho mayor que la mayoría, lo que podría ser relevante para un análisis más detallado.
Para el Gizscore, se observa una distribución algo más extendida. Al igual que en el caso anterior, se presentan varios outliers en los valores más altos, indicando una dispersión considerable en los puntajes de Gizscore.
En los boxplots de longitud y latitud, no hay valores atípicos. Las cajas de bigotes están bastante concentradas, lo que indica que los datos son más homogéneos en esas variables.
Aplicar el IQR (Rango Intercuartílico) para detectar los outliers de manera más formal y ver si encontramos algún valor adicional:
Usaremos el rango intercuartílico (IQR) para detectar los outliers de forma automática. El IQR es la diferencia entre el cuartil 3 (Q3) y cuartil 1 (Q1) de los datos.
Si un valor está por debajo de Q1 - 1.5 x IQR o por encima de Q3 + 1.5 x IQR, se considera un outlier.
# Paso 3: Detectar outliers usando IQR
# Función para encontrar los outliers
find_outliers <- function(variable) {
Q1 <- quantile(variable, 0.25) # Primer cuartil
Q3 <- quantile(variable, 0.75) # Tercer cuartil
IQR <- Q3 - Q1 # Rango intercuartílico
# Detectar outliers
outliers <- which(variable < (Q1 - 1.5 * IQR) | variable > (Q3 + 1.5 * IQR))
return(outliers) # Devuelve los índices de los outliers
}
# Detectamos outliers en las variables de interés
outliers_fallecidos <- find_outliers(datos$fallecidos)
outliers_gizscore <- find_outliers(datos$gizscore)
outliers_longitud <- find_outliers(datos$longitud)
outliers_latitud <- find_outliers(datos$latitud)
# Mostrar los índices de los outliers
outliers_fallecidos## [1] 28 96 133 176 184 191 200 203 217 219 240 274 295 297 313outliers_gizscore## [1] 28 34 133 176 184 200 203 217 219 240 295 313outliers_longitud## integer(0)outliers_latitud## integer(0)Análisis e interpretación:
Se identificaron 15 observaciones atípicas en la variable fallecidos y 12 en gizscore, lo cual coincide con lo observado visualmente en los boxplots. Estas observaciones corresponden a siniestros con una cantidad de fallecidos o una puntuación inusualmente alta, que se alejan del patrón general del resto del dataset.
En cambio, las variables longitud y latitud no presentan outliers, lo que sugiere que las ubicaciones geográficas de los siniestros están bastante concentradas dentro de un rango típico y no hay coordenadas anómalas que generen sospechas.
Decisión sobre Outliers
Los outliers identificados en fallecidos y
gizscore no se eliminan, ya que:
1. Representan la suma de fallecidos a lo largo de varios años en
un mismo tramo (valores altos plausibles).
2. Reflejan de manera válida la concentración espacial de
siniestralidad.
3. Conservarlos permite resaltar los tramos con mayor riesgo y no
sesgar el análisis descriptivo.
Nota: Dejar estos valores extremos ayuda a identificar los tramos más críticos y orientar posibles intervenciones de seguridad vial.
3.Transformaciones de Datos
En muchas ocasiones, las variables numéricas presentan distribuciones muy sesgadas o con colas largas, lo que afecta tanto las medidas de tendencia central (como la media) como las técnicas posteriores (correlaciones, PCA, regresiones).
Ya limpiamos los missing y detectamos (y decidimos conservar) los outliers. Para que nuestras estadísticas y gráficos no estén sesgados por distribuciones muy asimétricas o valores extremos, transformamos ciertas variables antes de calcular medias, hacer correlaciones o PCA. En nuestro caso, aplicaremos dos tipos de transformaciones sobre las variables que presentan colas largas (fallecidos y gizscore) y luego estandarizaremos todas las variables numéricas (fallecidos, gizscore, longitud, latitud) para dejarlas con media 0 y desviación estándar 1. Esto nos preparará para los análisis descriptivos y multivariados posteriores:
# 1) Logaritmo y raíz cuadrada para estabilizar varianzas
datos$log_fallecidos <- log(datos$fallecidos) # convierte la escala a logaritmo
datos$sqrt_fallecidos <- sqrt(datos$fallecidos) # convierte la escala a raíz cuadrada
datos$log_gizscore <- log(datos$gizscore) # idem para gizscore
datos$sqrt_gizscore <- sqrt(datos$gizscore)
# 2) Estandarización (z-score) de variables numéricas
vars_num <- c("fallecidos","gizscore","longitud","latitud")
datos_z <- scale(datos[, vars_num]) # media 0, sd 1
# Guardamos las columnas estandarizadas con sufijo "_z"
colnames(datos_z) <- paste0(vars_num, "_z")
datos <- cbind(datos, datos_z)Transformación de Variables Continuas
Con el objetivo de mejorar la simetría de las distribuciones y
mitigar la influencia de valores extremos, se aplicaron transformaciones
logarítmicas y de raíz cuadrada sobre las variables continuas
fallecidos y gizscore.
A continuación, se presentan los histogramas comparativos de las versiones originales y transformadas:
# Creamos histogramas antes y después para 'fallecidos'
hist_f1 <- ggplot(datos, aes(x = fallecidos)) +
geom_histogram(fill = "#FF9999", bins = 30) +
ggtitle("Fallecidos - Original")
hist_f2 <- ggplot(datos, aes(x = log_fallecidos)) +
geom_histogram(fill = "#FF6666", bins = 30) +
ggtitle("Fallecidos - Logarítmica")
hist_f3 <- ggplot(datos, aes(x = sqrt_fallecidos)) +
geom_histogram(fill = "#FF3333", bins = 30) +
ggtitle("Fallecidos - Raíz cuadrada")
# Histogramas para 'gizscore'
hist_g1 <- ggplot(datos, aes(x = gizscore)) +
geom_histogram(fill = "#99CCFF", bins = 30) +
ggtitle("GIZ Score - Original")
hist_g2 <- ggplot(datos, aes(x = log_gizscore)) +
geom_histogram(fill = "#6699FF", bins = 30) +
ggtitle("GIZ Score - Logarítmica")
hist_g3 <- ggplot(datos, aes(x = sqrt_gizscore)) +
geom_histogram(fill = "#3366FF", bins = 30) +
ggtitle("GIZ Score - Raíz cuadrada")
# Mostrar todos los histogramas en conjunto
grid.arrange(hist_f1, hist_f2, hist_f3,
hist_g1, hist_g2, hist_g3,
ncol = 3)
datos$log_fallecidos_z <- scale(datos$log_fallecidos)
datos$log_gizscore_z <- scale(datos$log_gizscore)
datos$longitud_z <- scale(datos$longitud)
datos$latitud_z <- scale(datos$latitud)Fallecidos
Distribución Original: Se observa una fuerte asimetría positiva, con muchos valores cercanos a cero y pocos valores extremadamente altos. Esto indica una distribución sesgada, lo que puede afectar los supuestos de normalidad requeridos para ciertos análisis estadísticos.
Transformación Logarítmica: Redujo significativamente la asimetría, concentrando más los datos hacia el centro. Aunque aún se nota cierto sesgo, la forma es más simétrica que la original.
Transformación de Raíz Cuadrada: También mejora la simetría, pero en menor medida que la transformación logarítmica. La distribución continúa mostrando colas más pesadas.
Conclusión: Para la variable
fallecidos, la transformación logarítmica parece ser la más
adecuada, ya que reduce de forma más efectiva el sesgo.
GIZ Score
Distribución Original: La variable
gizscoretambién presenta un sesgo positivo fuerte, con la mayoría de los datos concentrados en valores bajos y algunos casos extremos.Transformación Logarítmica: Disminuye la asimetría y produce una distribución más balanceada, aunque se siguen observando colas.
Transformación de Raíz Cuadrada: Suaviza el sesgo, pero no tanto como la logarítmica. La concentración de frecuencias sigue alta en valores bajos.
Conclusión: La transformación logarítmica ofrece
una mejor normalización de la variable gizscore que la raíz
cuadrada.
Decisión Final
Dado que el objetivo es acercar las distribuciones a una forma
más simétrica y manejable para análisis como el Análisis de Componentes
Principales (PCA) o modelos que asumen normalidad, se optará por
conservar las versiones logarítmicas de ambas variables
(log_fallecidos y log_gizscore)
para los análisis posteriores.*
Las versiones originales y con raíz cuadrada se conservarán únicamente como referencia, pero no se utilizarán directamente en los modelos.
4. Estadística Descriptiva
En esta sección se realizara un análisis exploratorio de las variables cuantitativas del conjunto de datos, con el objetivo de comprender su comportamiento, distribución y posibles transformaciones necesarias para futuros análisis multivariados.
Variables cualitativas:
cat_vars <- c("entidad", "departamento", "municipio", "tramo", "nombre")
# Creamos una lista vacía para guardar todas las tablas
resumenes_cualitativos <- list()
# Generamos las tablas de frecuencia
for (var in cat_vars) {
freq_abs <- table(datos[[var]])
freq_rel <- prop.table(freq_abs)
freq_abs_acum <- cumsum(freq_abs)
freq_rel_acum <- cumsum(freq_rel)
resumen <- data.frame(
Categoria = names(freq_abs),
Frec_Absoluta = as.numeric(freq_abs),
Frec_Relativa = round(as.numeric(freq_rel), 4),
Frec_Abs_Acum = as.numeric(freq_abs_acum),
Frec_Rel_Acum = round(as.numeric(freq_rel_acum), 4)
)
# Guardamos cada tabla en una lista con nombre
resumenes_cualitativos[[var]] <- resumen
}
# Imprimir cada tabla con título bonito
for (var in names(resumenes_cualitativos)) {
cat("\n\n### Frecuencia de la variable:", var, "\n\n")
print(kable(head(resumenes_cualitativos[[var]], 10),
caption = paste("Resumen de", var)))
}##
##
## ### Frecuencia de la variable: entidad
##
##
##
## Table: Resumen de entidad
##
## |Categoria | Frec_Absoluta| Frec_Relativa| Frec_Abs_Acum| Frec_Rel_Acum|
## |:------------|-------------:|-------------:|-------------:|-------------:|
## |ANI | 141| 0.4462| 141| 0.4462|
## |CIUDAD | 1| 0.0032| 142| 0.4494|
## |DEPTO | 48| 0.1519| 190| 0.6013|
## |INVIAS-OTROS | 126| 0.3987| 316| 1.0000|
##
##
## ### Frecuencia de la variable: departamento
##
##
##
## Table: Resumen de departamento
##
## |Categoria | Frec_Absoluta| Frec_Relativa| Frec_Abs_Acum| Frec_Rel_Acum|
## |:------------|-------------:|-------------:|-------------:|-------------:|
## |ANTIOQUIA | 35| 0.1108| 35| 0.1108|
## |ARAUCA | 5| 0.0158| 40| 0.1266|
## |ATLÁNTICO | 8| 0.0253| 48| 0.1519|
## |BOLÍVAR | 7| 0.0222| 55| 0.1741|
## |BOYACÁ | 4| 0.0127| 59| 0.1867|
## |CAQUETÁ | 2| 0.0063| 61| 0.1930|
## |CAUCA | 19| 0.0601| 80| 0.2532|
## |CESAR | 22| 0.0696| 102| 0.3228|
## |CÓRDOBA | 15| 0.0475| 117| 0.3703|
## |CUNDINAMARCA | 41| 0.1297| 158| 0.5000|
##
##
## ### Frecuencia de la variable: municipio
##
##
##
## Table: Resumen de municipio
##
## |Categoria | Frec_Absoluta| Frec_Relativa| Frec_Abs_Acum| Frec_Rel_Acum|
## |:---------------|-------------:|-------------:|-------------:|-------------:|
## |ACACÍAS | 1| 0.0032| 1| 0.0032|
## |AGUACHICA | 1| 0.0032| 2| 0.0063|
## |ALGARROBO | 1| 0.0032| 3| 0.0095|
## |ALVARADO | 1| 0.0032| 4| 0.0127|
## |ANDALUCÍA | 1| 0.0032| 5| 0.0158|
## |APARTADÓ | 3| 0.0095| 8| 0.0253|
## |ARJONA | 1| 0.0032| 9| 0.0285|
## |BARBOSA | 2| 0.0063| 11| 0.0348|
## |BARRANCABERMEJA | 1| 0.0032| 12| 0.0380|
## |BELLO | 14| 0.0443| 26| 0.0823|
##
##
## ### Frecuencia de la variable: tramo
##
##
##
## Table: Resumen de tramo
##
## |Categoria | Frec_Absoluta| Frec_Relativa| Frec_Abs_Acum| Frec_Rel_Acum|
## |:-----------------------------------------|-------------:|-------------:|-------------:|-------------:|
## |Aguachica - San Roque | 6| 0.0190| 6| 0.0190|
## |Arboletes-Montería | 2| 0.0063| 8| 0.0253|
## |Barbosa - Bucaramanga | 5| 0.0158| 13| 0.0411|
## |Barranca - Aguachica | 9| 0.0285| 22| 0.0696|
## |Barranquilla - Cartagena | 2| 0.0063| 24| 0.0759|
## |Barranquilla - Malambo - Palmar de Varela | 4| 0.0127| 28| 0.0886|
## |Bogotá - Barbosa | 1| 0.0032| 29| 0.0918|
## |Bogotá - Mosquera | 2| 0.0063| 31| 0.0981|
## |Bogotá - Tunja | 4| 0.0127| 35| 0.1108|
## |Bogotá - Villavicencio | 4| 0.0127| 39| 0.1234|
##
##
## ### Frecuencia de la variable: nombre
##
##
##
## Table: Resumen de nombre
##
## |Categoria | Frec_Absoluta| Frec_Relativa| Frec_Abs_Acum| Frec_Rel_Acum|
## |:----------------------------------------------------------|-------------:|-------------:|-------------:|-------------:|
## |Armenia -Pereira - Manizales | 2| 0.0063| 2| 0.0063|
## |Autopistas al Mar 2 | 10| 0.0316| 12| 0.0380|
## |Bogotá – La Vega – Villeta (BVLPV) | 3| 0.0095| 15| 0.0475|
## |Bogotá (Fontibón) – Faca – Los Alpes | 5| 0.0158| 20| 0.0633|
## |Briceño - Tunja - Sogamoso | 8| 0.0253| 28| 0.0886|
## |Cartagena – Barranquilla – “Circunvalar de la Prosperidad” | 2| 0.0063| 30| 0.0949|
## |Córdoba – Sucre | 8| 0.0253| 38| 0.1203|
## |Desarrollo Vial del Oriente de Medellín – DEVIMED | 9| 0.0285| 47| 0.1487|
## |Girardot – Ibagué – Cajamarca | 3| 0.0095| 50| 0.1582|
## |IP – Antioquia – Bolívar | 8| 0.0253| 58| 0.1835|cat_summary <- data.frame(
Variable = character(),
Moda = character(),
stringsAsFactors = FALSE
)
for (var in cat_vars) {
moda <- names(which.max(table(datos[[var]])))
cat_summary <- rbind(cat_summary, data.frame(
Variable = var,
Moda = moda
))
}
kable(cat_summary, caption = "Moda de variables cualitativas")| Variable | Moda |
|---|---|
| entidad | ANI |
| departamento | VALLE DEL CAUCA |
| municipio | BELLO |
| tramo | El Hatillo - Girardota - Copacabana - Bello - Med* |
| nombre | IP – Ampliación a Tercer Carril Doble Calzada Bogotá - Girardot |
—————————————————————————————————————————————————————————————————–
Para representar visualmente las modas de las variables cualitativas, se creó un gráfico de barras donde se muestra la categoría más frecuente de cada variable junto con su frecuencia absoluta. Esto permite ver rápidamente qué entidades, departamentos, municipios y tramos tienen mayor concentración de siniestros, siendo un insumo clave para identificar patrones y priorizar zonas de atención:
# Cambiar colores a tonos pastel
pastel_colors <- c("#F7A8B8", "#F3E0D2", "#B5D9A2", "#D2A6F5", "#B9E3F2")
# Crear tabla con las modas y sus frecuencias absolutas
modas_df <- data.frame(
Variable = c("Entidad", "Departamento", "Municipio", "Tramo", "Nombre"),
Categoria = c("ANI", "VALLE DEL CAUCA", "BELLO",
"El Hatillo-Girardota-Copacabana-Bello-Med*",
"IP–AmpliaciónaTercerCarrilDobleCalzadaBogotá-Girardot"),
Frecuencia = c(141, 60, 14, 9, 189) # Frecuencia ajustada para tramo y nombre
)
# Crear gráfico de barras con colores pastel y la moda más frecuente
ggplot(modas_df, aes(x = Variable, y = Frecuencia, fill = Variable)) +
geom_col(show.legend = FALSE) +
geom_text(aes(label = Categoria), vjust = -0.5, size = 3, angle = 15, color = "black") +
scale_fill_manual(values = pastel_colors) +
labs(
title = "Moda por Variable Cualitativa",
subtitle = "Categorías más frecuentes en cada variable",
x = "Variable",
y = "Frecuencia Absoluta"
) +
theme_minimal() +
theme(plot.title = element_text(size = 14, face = "bold"))
En la parte de estadística descriptiva para las variables cualitativas, hicimos un resumen completo usando tablas de frecuencia para entender mejor cómo se distribuyen los datos. Se incluyeron las frecuencias absolutas y relativas, además de las acumuladas (tanto absolutas como relativas), lo que nos permite ver no solo cuántas veces aparece cada categoría, sino también cómo se van acumulando en el total. Para no hacerlo tan largo, mostramos solo las 10 categorías más comunes en cada variable (como entidad, departamento, municipio, tramo y nombre). Esto es muy útil porque nos deja ver rápidamente cuáles son las más frecuentes, y nos ayuda a identificar patrones o concentraciones de datos. También calculamos la moda de cada variable, es decir, el valor que más se repite. Por ejemplo, la entidad más común es ANI, el departamento más frecuente es Valle del Cauca, el municipio con más registros es Bello, y en el caso de tramo y nombre, las categorías más repetidas fueron El Hatillo - Girardota - Copacabana - Bello - Med y IP – Ampliación a Tercer Carril Doble Calzada Bogotá - Girardot, respectivamente. Esto nos da una idea clara de dónde y en qué proyectos se concentra la mayor parte de los datos del conjunto, lo cual es clave para futuros análisis o decisiones.
Variables cuantitativas:
# Variables cuantitativas (usa las que realmente te interesan y tienen sentido)
quant_vars <- c("fallecidos", "gizscore", "latitud", "longitud",
"log_fallecidos", "sqrt_fallecidos", "log_gizscore",
"sqrt_gizscore", "fallecidos_z", "gizscore_z",
"longitud_z", "latitud_z")
# Creamos un dataframe vacío para guardar los resultados
summary_quant <- data.frame(
Variable = character(),
Media = numeric(),
Mediana = numeric(),
Moda = numeric(),
Varianza = numeric(),
Desviacion = numeric(),
Rango = numeric(),
P25 = numeric(),
P50 = numeric(),
P75 = numeric(),
stringsAsFactors = FALSE
)
# Función para calcular la moda
get_moda <- function(x) {
t <- table(x)
as.numeric(names(t)[which.max(t)])
}
# Llenamos el dataframe
for (var in quant_vars) {
x <- datos[[var]]
media <- mean(x, na.rm = TRUE)
mediana <- median(x, na.rm = TRUE)
moda <- get_moda(x)
varianza <- var(x, na.rm = TRUE)
desviacion <- sd(x, na.rm = TRUE)
rango <- max(x, na.rm = TRUE) - min(x, na.rm = TRUE)
p25 <- quantile(x, 0.25, na.rm = TRUE)
p50 <- quantile(x, 0.50, na.rm = TRUE)
p75 <- quantile(x, 0.75, na.rm = TRUE)
summary_quant <- rbind(summary_quant, data.frame(
Variable = var,
Media = round(media, 2),
Mediana = round(mediana, 2),
Moda = round(moda, 2),
Varianza = round(varianza, 2),
Desviacion = round(desviacion, 2),
Rango = round(rango, 2),
P25 = round(p25, 2),
P50 = round(p50, 2),
P75 = round(p75, 2)
))
}
# Mostrar todo bien bonito
kable(summary_quant, caption = "Resumen estadístico de variables cuantitativas")| Variable | Media | Mediana | Moda | Varianza | Desviacion | Rango | P25 | P50 | P75 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 25% | fallecidos | 11.33 | 10.00 | 9.00 | 50.02 | 7.07 | 51.00 | 7.00 | 10.00 | 14.00 |
| 25%1 | gizscore | 3.19 | 2.70 | 1.81 | 2.67 | 1.64 | 12.88 | 2.11 | 2.70 | 3.88 |
| 25%2 | latitud | 5.67 | 4.75 | 3.73 | 7.55 | 2.75 | 10.70 | 3.68 | 4.75 | 7.80 |
| 25%3 | longitud | -75.17 | -75.43 | -76.64 | 2.07 | 1.44 | 7.06 | -76.32 | -75.43 | -73.95 |
| 25%4 | log_fallecidos | 2.19 | 2.30 | 2.20 | 0.62 | 0.79 | 3.95 | 1.95 | 2.30 | 2.64 |
| 25%5 | sqrt_fallecidos | 3.20 | 3.16 | 3.00 | 1.11 | 1.05 | 6.21 | 2.65 | 3.16 | 3.74 |
| 25%6 | log_gizscore | 1.07 | 0.99 | 0.59 | 0.17 | 0.41 | 2.17 | 0.75 | 0.99 | 1.36 |
| 25%7 | sqrt_gizscore | 1.74 | 1.64 | 1.35 | 0.16 | 0.39 | 2.52 | 1.45 | 1.64 | 1.97 |
| 25%8 | fallecidos_z | 0.00 | -0.19 | -0.33 | 1.00 | 1.00 | 7.21 | -0.61 | -0.19 | 0.38 |
| 25%9 | gizscore_z | 0.00 | -0.30 | -0.84 | 1.00 | 1.00 | 7.88 | -0.66 | -0.30 | 0.42 |
| 25%10 | longitud_z | 0.00 | -0.18 | -1.02 | 1.00 | 1.00 | 4.91 | -0.80 | -0.18 | 0.85 |
| 25%11 | latitud_z | 0.00 | -0.34 | -0.71 | 1.00 | 1.00 | 3.89 | -0.73 | -0.34 | 0.77 |
————————————————————————————————————————————————————-
se presentan histogramas para algunas de las variables cuantitativas más importantes, como fallecidos y gizscore, lo que permite una mejor comprensión visual de su distribución y frecuencia:
# Histograma de fallecidos
ggplot(datos, aes(x = fallecidos)) +
geom_histogram(binwidth = 1, fill = "#FF9999", color = "black", alpha = 0.7) +
labs(title = "Distribución de Fallecidos", x = "Fallecidos", y = "Frecuencia") +
theme_minimal()
# Histograma de Gizscore
ggplot(datos, aes(x = gizscore)) +
geom_histogram(binwidth = 0.5, fill = "#99CCFF", color = "black", alpha = 0.7) +
labs(title = "Distribución de Gizscore", x = "Gizscore", y = "Frecuencia") +
theme_minimal()
En esta parte hicimos un resumen completo de las variables cuantitativas del conjunto de datos, usando medidas como la media, mediana, moda, varianza, desviación estándar, rango y percentiles (P25, P50 y P75). Estas medidas nos ayudaron a entender cómo se comportan los datos en general, si están concentrados, dispersos o si hay valores extremos. Por ejemplo, vimos que algunas variables como fallecidos y gizscore tenían diferencias entre la media, la mediana y la moda, lo que muestra que los datos están un poco sesgados y no son tan simétricos. Por eso también hicimos transformaciones como el logaritmo y la raíz cuadrada para que los datos fueran más estables y comparables. Además, normalizamos varias variables (las que terminan en _z) para que todas estuvieran en una misma escala, lo cual es útil para análisis más avanzados. Que algunas variables tengan valores negativos no es un error, simplemente significa que esos datos están por debajo del promedio. En general, este análisis nos dio una buena base para entender mejor los datos y prepararlos para los siguientes pasos del estudio.
5. Análisis Bidimensional y Correlaciones
Ahora vamos a cruzar variables entre sí para ver si existe alguna relación. Esto nos va a ayudar a identificar, por ejemplo, si en ciertos departamentos hay más fallecidos, o si ciertas entidades están relacionadas con mayor siniestralidad.
También es útil para ver si hay correlación numérica directa entre variables, como por ejemplo:
gizscore vs. fallecidos
latitud vs. fallecidos
longitud vs. gizscore
Matriz de Correlaciones
Vamos a iniciar el análisis bidimensional observando cómo se relacionan entre sí algunas variables numéricas clave, ya transformadas y estandarizadas. Para esto usamos la matriz de correlación de Pearson, que nos permite identificar relaciones lineales.
Se incluyen: log_fallecidos_z,
log_gizscore_z, longitud_z y
latitud_z.
El gráfico resultante (corrplot) nos muestra qué pares de variables tienen relaciones positivas o negativas más fuertes. Esto nos ayuda a decidir en qué relaciones profundizar con gráficos de dispersión:
# Seleccionamos las variables z-transformadas
vars_corr <- datos %>%
select(log_fallecidos_z, log_gizscore_z, longitud_z, latitud_z)
# Calculamos la matriz de correlaciones
matriz_cor <- cor(vars_corr, use = "complete.obs", method = "pearson")
# Visualizamos con un mapa de calor
corrplot(matriz_cor, method = "color", type = "upper",
tl.col = "black", tl.srt = 45,
addCoef.col = "black", number.cex = 0.7,
col = colorRampPalette(c("#B3E5FC", "#0288D1"))(200),
mar = c(0, 0, 1, 0))
Análisis de la matriz de correlaciones
En la matriz de correlación observamos cómo se relacionan entre sí las variables numéricas estandarizadas (z-score). La relación más fuerte se da entre log_fallecidos_z y log_gizscore_z, con una correlación de 0.51, lo que indica una asociación positiva moderada: a mayor puntuación en GIZ (zonas más críticas), tiende a haber más fallecidos. Esto es coherente, ya que el GIZ mide la peligrosidad de los tramos viales.
Las demás relaciones son bastante débiles o casi nulas. Por ejemplo, log_fallecidos_z y longitud_z tienen una correlación de solo 0.01, lo que sugiere que la ubicación longitudinal no influye en la cantidad de fallecidos. La única relación un poco más destacada aparte de la anterior es entre longitud_z y latitud_z (0.53), pero esto tiene más que ver con la geografía natural del país (ubicaciones tienden a estar alineadas espacialmente), y no con los accidentes en sí.
Este análisis nos muestra que la variable más útil para relacionar con fallecidos es log_gizscore_z. Será clave representarla luego en los diagramas de dispersión para profundizar más.
————————————————————————————————————————————————————-
Tablas de contingencia entre variables cualitativas
Es una tabla que cruza dos variables cualitativas para ver cómo se distribuyen juntas. Por ejemplo, cuántos tramos en el departamento Valle del Cauca están a cargo de ANI, o si cierta entidad aparece más en ciertos municipios.
# Tabla de frecuencias absolutas entre entidad y departamento
tabla_contingencia <- table(datos$entidad, datos$departamento)
tabla_contingencia##
## ANTIOQUIA ARAUCA ATLÁNTICO BOLÍVAR BOYACÁ CAQUETÁ CAUCA CESAR
## ANI 19 0 7 6 4 0 12 6
## CIUDAD 0 0 0 0 0 0 0 1
## DEPTO 16 1 0 0 0 0 0 0
## INVIAS-OTROS 0 4 1 1 0 2 7 15
##
## CÓRDOBA CUNDINAMARCA HUILA LA GUAJIRA MAGDALENA META NARIÑO
## ANI 13 30 2 5 6 11 5
## CIUDAD 0 0 0 0 0 0 0
## DEPTO 0 5 0 0 4 0 0
## INVIAS-OTROS 2 6 4 1 0 3 14
##
## NORTE DE SANTANDER QUINDÍO RISARALDA SANTANDER SUCRE TOLIMA
## ANI 0 2 0 0 3 10
## CIUDAD 0 0 0 0 0 0
## DEPTO 1 0 0 8 0 0
## INVIAS-OTROS 3 4 2 9 1 0
##
## VALLE DEL CAUCA
## ANI 0
## CIUDAD 0
## DEPTO 13
## INVIAS-OTROS 47# Tabla de proporciones (por fila)
prop_por_fila <- prop.table(tabla_contingencia, margin = 1)
round(prop_por_fila, 2) # Redondeamos a 2 decimales##
## ANTIOQUIA ARAUCA ATLÁNTICO BOLÍVAR BOYACÁ CAQUETÁ CAUCA CESAR
## ANI 0.13 0.00 0.05 0.04 0.03 0.00 0.09 0.04
## CIUDAD 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
## DEPTO 0.33 0.02 0.00 0.00 0.00 0.00 0.00 0.00
## INVIAS-OTROS 0.00 0.03 0.01 0.01 0.00 0.02 0.06 0.12
##
## CÓRDOBA CUNDINAMARCA HUILA LA GUAJIRA MAGDALENA META NARIÑO
## ANI 0.09 0.21 0.01 0.04 0.04 0.08 0.04
## CIUDAD 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## DEPTO 0.00 0.10 0.00 0.00 0.08 0.00 0.00
## INVIAS-OTROS 0.02 0.05 0.03 0.01 0.00 0.02 0.11
##
## NORTE DE SANTANDER QUINDÍO RISARALDA SANTANDER SUCRE TOLIMA
## ANI 0.00 0.01 0.00 0.00 0.02 0.07
## CIUDAD 0.00 0.00 0.00 0.00 0.00 0.00
## DEPTO 0.02 0.00 0.00 0.17 0.00 0.00
## INVIAS-OTROS 0.02 0.03 0.02 0.07 0.01 0.00
##
## VALLE DEL CAUCA
## ANI 0.00
## CIUDAD 0.00
## DEPTO 0.27
## INVIAS-OTROS 0.37# Tabla de proporciones (por columna)
prop_por_columna <- prop.table(tabla_contingencia, margin = 2)
round(prop_por_columna, 2)##
## ANTIOQUIA ARAUCA ATLÁNTICO BOLÍVAR BOYACÁ CAQUETÁ CAUCA CESAR
## ANI 0.54 0.00 0.88 0.86 1.00 0.00 0.63 0.27
## CIUDAD 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05
## DEPTO 0.46 0.20 0.00 0.00 0.00 0.00 0.00 0.00
## INVIAS-OTROS 0.00 0.80 0.12 0.14 0.00 1.00 0.37 0.68
##
## CÓRDOBA CUNDINAMARCA HUILA LA GUAJIRA MAGDALENA META NARIÑO
## ANI 0.87 0.73 0.33 0.83 0.60 0.79 0.26
## CIUDAD 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## DEPTO 0.00 0.12 0.00 0.00 0.40 0.00 0.00
## INVIAS-OTROS 0.13 0.15 0.67 0.17 0.00 0.21 0.74
##
## NORTE DE SANTANDER QUINDÍO RISARALDA SANTANDER SUCRE TOLIMA
## ANI 0.00 0.33 0.00 0.00 0.75 1.00
## CIUDAD 0.00 0.00 0.00 0.00 0.00 0.00
## DEPTO 0.25 0.00 0.00 0.47 0.00 0.00
## INVIAS-OTROS 0.75 0.67 1.00 0.53 0.25 0.00
##
## VALLE DEL CAUCA
## ANI 0.00
## CIUDAD 0.00
## DEPTO 0.22
## INVIAS-OTROS 0.78Análisis:
ANI (Agencia Nacional de Infraestructura) es la entidad con mayor presencia general, especialmente en departamentos como Cundinamarca (30 tramos), Antioquia, Córdoba y Tolima.
INVIAS-OTROS, por su parte, tiene una fuerte presencia en Valle del Cauca (47 tramos), Nariño (14), Cauca, y otros más dispersos.
DEPTO (probablemente entidades departamentales) opera principalmente en Antioquia (16), Cundinamarca, y Valle del Cauca.
°Cuando observamos la distribución relativa:
Por fila, vemos cómo cada entidad se reparte entre los departamentos. Por ejemplo, el 27% de los tramos operados por DEPTOS están en Valle del Cauca, mientras que el 37% de los tramos de INVIAS-OTROS están también allí.
Por columna, se nota qué entidad domina cada departamento. Ejemplo: INVIAS-OTROS gestiona el 78% de los tramos en el Valle del Cauca, mientras que ANI gestiona el 87% en Córdoba.
Nota: El test chi-cuadrado te dirá si hay dependencia entre ambas variables (es decir, si una influye sobre la otra o si se distribuyen de manera independiente):
# Test de independencia entre entidad y departamento
chisq.test(tabla_contingencia)##
## Pearson's Chi-squared test
##
## data: tabla_contingencia
## X-squared = 235.5, df = 63, p-value < 2.2e-16Interpretación:
° El p-value es mucho menor a 0.05, lo que indica que rechazamos la hipótesis nula de independencia. Es decir, sí existe una relación significativa entre la entidad responsable y el departamento donde ocurre el siniestro.
° En palabras simples: la distribución de las entidades no es aleatoria entre departamentos, sino que hay un patrón claro de asociación (como ya vimos en la tabla).
° La advertencia (“Chi-squared approximation may be incorrect”) solo indica que algunas frecuencias en la tabla son bajas, lo que puede hacer que el test no sea 100% preciso. Pero no es grave si los tamaños de muestra son medianamente grandes, como aquí.
6. Gráficos Multivariados
Diagramas de dispersión con línea de tendencia lineal
En este paso buscamos entender visualmente si existe una relación entre la cantidad de fallecidos y otras variables numéricas, usando gráficos de dispersión acompañados por líneas de regresión lineal.
¿Cómo lo haremos? Se utilizaran dos scatterplots con las siguientes combinaciones:
1. GIZScore (log y z-score) vs. Fallecidos (log y z-score) 2. Latitud (z-score) vs. Fallecidos (log y z-score)
Ambas variables fueron previamente transformadas y estandarizadas para permitir una mejor comparación visual y facilitar el análisis estadístico:
ggplot(datos, aes(x = log_gizscore_z, y = log_fallecidos_z)) +
geom_point(alpha = 0.6, size = 2, color = "#69b3a2") +
geom_smooth(method = "lm", se = TRUE, color = "#404080") +
labs(
title = "Relación entre GIZScore y Fallecidos (log z score)",
x = "GIZScore (log y estandarizado)",
y = "Fallecidos (log y estandarizado)"
) +
theme_minimal()
ggplot(datos, aes(x = latitud_z, y = log_fallecidos_z)) +
geom_point(alpha = 0.6, size = 2, color = "#F7A8B8") +
geom_smooth(method = "lm", se = TRUE, color = "#803040") +
labs(
title = "Relación entre Latitud y Fallecidos (z‑score)",
x = "Latitud (estandarizada)",
y = "Fallecidos (log y estandarizado)"
) +
theme_minimal()
————————————————————————————————————————————————————-
Scatterplot multivariado (tipo burbuja)
Este gráfico nos permite ver dónde ocurren los siniestros, cuántos fallecidos hubo y qué entidad estaba a cargo:
Esto te permite ver:
Dónde ocurren los siniestros (ubicación).
Cuántos fallecidos hubo en cada punto (tamaño).
Qué entidad estaba a cargo (color).
# Filtrar las dos entidades más frecuentes (opcional)
top2 <- names(sort(table(datos$entidad), TRUE))[1:2]
datos2 <- datos %>% filter(entidad %in% top2)
ggplot(datos2, aes(x = longitud, y = latitud, size = fallecidos, color = entidad)) +
geom_point(alpha = 0.5) +
scale_color_brewer(palette = "Set2") +
scale_size_continuous(range = c(2, 10), breaks = c(10, 30, 50)) +
labs(
title = "Siniestros viales por Entidad (Top 2)",
subtitle = "Tamaño: fallecidos | Color: entidad",
x = "Longitud", y = "Latitud"
) +
theme_minimal() +
theme(legend.position = "right")
————————————————————————————————————————————————————-
Boxplots comparativos (cuantitativa vs. cualitativa)
Vamos a comparar la distribución de fallecidos (o su versión transformada) entre las distintas categorías de una variable cualitativa. Por ejemplo:
Fallecidos por Entidad
Fallecidos por Departamento
Los boxplots nos mostrarán médianas, rangos y posibles outliers, ayudándonos a ver rápidamente si alguna entidad o departamento destaca por tener más o menos siniestralidad.
ggplot(datos, aes(x = entidad, y = fallecidos)) +
geom_boxplot(fill = "#A3C4DC", outlier.color = "#FF6F61") +
labs(
title = "Distribución de Fallecidos por Entidad",
x = "Entidad",
y = "Número de Fallecidos"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1)
)
Nota: Durante la exploración de datos se identificó una categoría “CIUDAD” en la variable
entidadcon un único registro, que no corresponde a ninguna de las entidades oficiales. Para mantener la consistencia y evitar distorsiones en el análisis comparativo, se procedera a eliminar esta observación:
# Eliminar la fila donde entidad == "CIUDAD"
datos <- datos %>%
filter(entidad != "CIUDAD")ggplot(datos, aes(x = departamento, y = fallecidos)) +
geom_boxplot(fill = "#C7CEEA", outlier.color = "#FF6F61") +
labs(
title = "Distribución de Fallecidos por Departamento",
x = "Departamento",
y = "Número de Fallecidos"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1, size = 6)
)
Interpretación
Boxplot por Entidad:
Observamos que la ANI tiende a tener un rango de fallecidos ligeramente más alto (mediana y cuartiles superior) que las demás entidades. Esto sugiere que los tramos a cargo de la ANI presentan, en promedio, más víctimas mortales.Boxplot por Departamento:
El Valle del Cauca muestra una mediana de fallecidos por tramo significativamente más elevada que otros departamentos, con un outlier en la parte alta. Esto indica puntos críticos de siniestralidad que requerirían atención prioritaria.
Mapa de calor de siniestralidad vial
ggplot(datos, aes(x = longitud, y = latitud)) +
stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = 0.6, contour = TRUE) +
geom_point(aes(x = longitud, y = latitud), color = "black", alpha = 0.3, size = 1) +
scale_fill_gradient(low = "#FFE5E5", high = "#B80000") +
labs(
title = "Mapa de calor de siniestralidad vial",
subtitle = "Zonas con mayor concentración de siniestros fatales",
x = "Longitud",
y = "Latitud",
fill = "Densidad"
) +
coord_fixed(ratio = 1.3) +
theme_minimal(base_size = 12)
Este gráfico nos muestra las zonas con mayor concentración de siniestros fatales en Colombia, usando coordenadas de longitud y latitud. Se empleó un mapa de calor espacial donde los tonos rojos más intensos indican mayor densidad de incidentes. Además, se sobrepusieron los puntos reales de los siniestros para visualizar con más claridad la distribución exacta.
Podemos observar claramente varios focos rojos, especialmente en zonas cercanas al centro y suroccidente del país, lo que indica que allí se presentan con mayor frecuencia los accidentes con fallecidos. Esta visualización ayuda a identificar zonas críticas que podrían priorizarse para intervenciones viales o mayor vigilancia.
El gráfico combina la suavidad del mapa de calor con la precisión de los puntos individuales, lo que permite un análisis más completo de los patrones espaciales de siniestralidad.
7. Reducción de Dimensiones (PCA/EFA)
Cuando tenemos muchas variables numéricas y algunas están correlacionadas entre sí, el análisis puede volverse redundante o difícil de interpretar. Aquí es donde entra la Reducción de Dimensiones. Con PCA y EFA podemos:
Detectar patrones más generales y compactos.
Simplificar sin perder información importante.
Visualizar agrupaciones o relaciones ocultas entre los datos.
Preparar el terreno para modelos más avanzados (opcional).
usaremos principalmente las variables numéricas transformadas y estandarizadas (las z-score), como:
variables_pca <- datos %>% select(log_fallecidos_z, log_gizscore_z, longitud_z, latitud_z)
Estas están listas porque ya:
Se transformaron (logarítmica para normalizar)
Se estandarizaron (media 0, desviación 1)
Están limpias y sin outliers extremos eliminados
Análisis de Componentes Principales (PCA)
# Seleccionar variables numéricas estandarizadas (log y z-score)
variables_pca <- datos %>%
select(log_fallecidos_z, log_gizscore_z, longitud_z, latitud_z)
# Ejecutar PCA
pca_resultado <- PCA(variables_pca, scale.unit = TRUE, graph = FALSE)
# Tabla con la varianza explicada por cada componente
pca_resultado$eig## eigenvalue percentage of variance cumulative percentage of variance
## comp 1 1.5932609 39.83152 39.83152
## comp 2 1.4736704 36.84176 76.67328
## comp 3 0.4821900 12.05475 88.72803
## comp 4 0.4508788 11.27197 100.00000Gráfico de varianza explicada (gráfico del codo)
fviz_eig(pca_resultado, addlabels = TRUE, ylim = c(0, 80)) +
labs(title = "Varianza explicada por componente (PCA)",
x = "Componentes principales", y = "Porcentaje de varianza")
Análisis:
El análisis de componentes principales (PCA) mostró que los dos primeros componentes explican el 76.75% de la varianza total del conjunto de datos (39.9% y 36.9% respectivamente). Esto indica que con solo dos componentes es posible representar gran parte de la información contenida en las variables originales (log_fallecidos_z, log_gizscore_z, longitud_z y latitud_z), reduciendo la dimensionalidad sin perder mucho contenido informativo.
Los componentes 3 y 4 aportan una varianza mucho menor (12.1% y 11.2%), por lo que no serían necesarios para análisis exploratorios o visuales básicos. En conclusión, se opta por conservar únicamente los dos primeros componentes principales, lo cual simplifica el análisis y permite una interpretación más clara.
————————————————————————————————————————————————————-
Mostrar la tabla de cargas de los componentes (loadings)
# Seleccionamos solo las variables numéricas estandarizadas
vars_pca <- datos[, c("log_fallecidos_z", "log_gizscore_z", "longitud_z", "latitud_z")]
# Ejecutamos PCA
pca_result <- princomp(vars_pca, cor = TRUE)
# Mostrar cargas factoriales de los dos primeros componentes
loadings <- pca_result$loadings
round(unclass(loadings)[, 1:2], 3)## Comp.1 Comp.2
## log_fallecidos_z 0.520 0.468
## log_gizscore_z 0.437 0.555
## longitud_z 0.455 -0.553
## latitud_z 0.576 -0.408¿Qué indica esto?
Componente 1 (Dim1 – 39.9% de varianza):
Las variables con mayor carga son latitud_z (0.582) y log_fallecidos_z (0.513), seguidas por longitud_z (0.464) y log_gizscore_z (0.426).
Interpretación: Este componente parece representar una combinación entre ubicación geográfica (latitud/longitud) y nivel de siniestralidad.
Componente 2 (Dim2 – 36.9% de varianza):
Las mayores cargas vienen de log_gizscore_z (0.564) y log_fallecidos_z (0.476), con longitud_z (-0.545) también fuerte, pero en sentido opuesto.
Interpretación: Este eje podría reflejar un contraste entre nivel de desarrollo (GIZScore) y localización longitudinal.
Ambas dimensiones tienen aportes similares y explican en conjunto más del 76% de la varianza total, lo que es excelente.
————————————————————————————————————————————————————-
Análisis final del PCA y visualización multivariada
fviz_pca_biplot(pca_result,
repel = TRUE,
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
col.ind = "gray40",
label = "var", # Solo muestra etiquetas de variables
title = "Biplot de Componentes Principales")
La imagen muestra que:
El gráfico de componentes principales (PCA) muestra cómo las
cuatro variables seleccionadas (log_fallecidos_z,
log_gizscore_z, latitud_z,
longitud_z) se agrupan y se relacionan entre sí en un nuevo
espacio reducido de dos dimensiones.
La Dimensión 1 (39.9%) parece capturar una
combinación de ubicación geográfica, ya que tanto latitud_z
como longitud_z apuntan hacia la misma dirección, indicando
una contribución importante. Esto sugiere que esta dimensión está
asociada a la posición espacial de los siniestros
Por otro lado, la Dimensión 2 (36.9%) está dominada por
log_fallecidos_z y log_gizscore_z, que también
apuntan juntas pero en una dirección diferente. Esta dimensión parece
asociarse más con la gravedad del siniestro o el nivel de impacto
(muertes y vulnerabilidad territorial).
En conjunto, estas dos dimensiones explican un 76.8% de la varianza total, lo cual es bastante alto y nos permite simplificar los datos sin perder demasiada información. El biplot nos ayuda a visualizar cómo se distribuyen los puntos (observaciones) según estas nuevas dimensiones y cómo cada variable influye en esa distribución.
Este tipo de gráfico nos permite detectar patrones espaciales y agrupaciones relacionadas con la severidad de los siniestros, que puede ser muy útil para futuras estrategias de prevención.
Análisis Factorial Exploratorio (AFE)
# ── selecciona solo las columnas numéricas que vas a factorizar ─────────────
vars_factor <- datos[, c("log_fallecidos_z",
"log_gizscore_z",
"longitud_z",
"latitud_z")]
# verifica que no haya NAs
sum(is.na(vars_factor)) ## [1] 0————————————————————————————————————————————————————-
Prueba de Esfericidad de Bartlett E Índice KMO (Kaiser-Meyer-Olkin):
# ── tests previos ────────────────────────────────────────────────────────────
cortest.bartlett(vars_factor) # Ho: matriz = identidad → p < .05 para seguir## $chisq
## [1] 209.6876
##
## $p.value
## [1] 1.641526e-42
##
## $df
## [1] 6KMO(vars_factor) # MSA global y por variable## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = vars_factor)
## Overall MSA = 0.5
## MSA for each item =
## log_fallecidos_z log_gizscore_z longitud_z latitud_z
## 0.51 0.50 0.49 0.50¿Qué significa?
Esta prueba verifica si las variables están lo suficientemente relacionadas como para aplicar AFE. El valor de p es muchísimo menor que 0.05, lo que indica que sí hay correlaciones significativas entre las variables.
El KMO mide si los datos son adecuados para el AFE. Un valor cercano a 1 es ideal. En este caso obtuvimos 0.5, que es el mínimo aceptable.
————————————————————————————————————————————————————-
Análisis Paralelo (gráfico de scree plot)
# ── parallel analysis ───────────────────────────────────────────────────────
fa.parallel(vars_factor, fa = "fa", n.iter = 100, show.legend = FALSE)
## Parallel analysis suggests that the number of factors = 2 and the number of components = NA¿Qué se ve?
Las dos primeras componentes (factores) tienen eigenvalores por encima de la línea roja → entonces se justifica mantener 2 factores
————————————————————————————————————————————————————-
Cargas, Comunalidades y Unicidades
# ── modelo factorial ────────────────────────────────────────────────────────
efa_res <- fa(vars_factor,
nfactors = 2, # cambia si tu paralelo sugiere 1
rotate = "varimax",# rotación ortogonal (carga simple)
fm = "ml") # método de estimación: máxima verosimilitud
print(efa_res, cut = 0.30) # muestra solo cargas ≥ .30## Factor Analysis using method = ml
## Call: fa(r = vars_factor, nfactors = 2, rotate = "varimax", fm = "ml")
## Standardized loadings (pattern matrix) based upon correlation matrix
## ML1 ML2 h2 u2 com
## log_fallecidos_z 0.72 0.52 0.48 1.0
## log_gizscore_z 0.72 0.51 0.49 1.0
## longitud_z 0.73 0.55 0.45 1.0
## latitud_z 0.73 0.55 0.45 1.1
##
## ML1 ML2
## SS loadings 1.09 1.05
## Proportion Var 0.27 0.26
## Cumulative Var 0.27 0.53
## Proportion Explained 0.51 0.49
## Cumulative Proportion 0.51 1.00
##
## Mean item complexity = 1
## Test of the hypothesis that 2 factors are sufficient.
##
## df null model = 6 with the objective function = 0.67 with Chi Square = 209.69
## df of the model are -1 and the objective function was 0
##
## The root mean square of the residuals (RMSR) is 0
## The df corrected root mean square of the residuals is NA
##
## The harmonic n.obs is 315 with the empirical chi square 0 with prob < NA
## The total n.obs was 315 with Likelihood Chi Square = 0 with prob < NA
##
## Tucker Lewis Index of factoring reliability = 1.03
## Fit based upon off diagonal values = 1
## Measures of factor score adequacy
## ML1 ML2
## Correlation of (regression) scores with factors 0.84 0.83
## Multiple R square of scores with factors 0.71 0.68
## Minimum correlation of possible factor scores 0.41 0.37¿Qué significa?
- Cargas (≥0.7) → indican una fuerte relación entre variable y factor.
- h2 (comunalidades) → qué tanto del comportamiento de esa variable está explicado por los factores. Como están todas por encima de 0.5, están bien explicadas.
- u2 (unicidad) → es el “resto”, lo que no se explica por los factores. Si es bajo, mejor.
————————————————————————————————————————————————————-
Análisis Paralelo (gráfico de scree plot)
# ── diagrama de factores ────────────────────────────────────────────────────
fa.diagram(efa_res,
main = "Mapa de Cargas Factoriales",
digits = 2, # redondeo
simple = FALSE) # muestra cargas cruzadas si las hay
ML1 agrupa la información geográfica → latitud_z y longitud_z
ML2 agrupa la siniestralidad → log_fallecidos_z y log_gizscore_z
Esto tiene mucho sentido lógico, porque separa los factores de ubicación de los relacionados con la gravedad de los siniestros.
————————————————————————————————————————————————————-
Interpretación final del EFA:
Para entender mejor la relación entre las variables y reducir la complejidad del análisis, aplicamos un Análisis Factorial Exploratorio (EFA) sobre las variables estandarizadas: log_fallecidos_z, log_gizscore_z, latitud_z y longitud_z.
Primero, comprobamos que era válido usar EFA: la prueba de esfericidad de Bartlett resultó significativa (p < 0.05), lo que indica que sí hay relaciones entre las variables. El índice KMO fue exactamente 0.5, el valor mínimo aceptado, así que aunque el modelo no es muy fuerte, sí es aplicable.
A partir del análisis paralelo (gráfico de líneas), se sugirió conservar 2 factores principales. Esta decisión fue coherente con los resultados, ya que las variables se agruparon de manera lógica:
- Factor 1 (ML1) agrupa latitud_z y longitud_z, lo que representa la localización geográfica de los siniestros.
- Factor 2 (ML2) agrupa log_fallecidos_z y log_gizscore_z, asociadas a la gravedad de los siniestros.
Las cargas factoriales de todas las variables fueron altas (alrededor de 0.7), lo que significa que los factores explican bien cada variable. Además, las comunalidades fueron todas mayores a 0.5, lo que indica que los factores explican más del 50% de la variabilidad de cada una.
En resumen, el EFA nos permitió identificar dos dimensiones clave en los datos: una relacionada con la ubicación de los siniestros y otra con su severidad. Esta reducción facilita el análisis posterior y permite una comprensión más clara de los patrones generales.
8. Modelo Predictivo:
Un modelo predictivo es una herramienta de análisis estadístico o de machine learning (aprendizaje automático) que utiliza datos históricos para predecir un resultado futuro. En otras palabras, predice lo que podría pasar bajo ciertas condiciones, basándose en patrones detectados en los datos previos.
¿Para qué serviría?
Este modelo serviría para obtener información predictiva sobre las situaciones de siniestralidad vial, lo que puede ser muy útil para la planificación de políticas públicas, prevención de accidentes, asignación de recursos en zonas de alto riesgo, etc.
# Crear el modelo de regresión lineal
modelo <- lm(fallecidos ~ log_gizscore_z + longitud_z + latitud_z + entidad + departamento, data = datos)
# Ver el resumen del modelo
summary(modelo)##
## Call:
## lm(formula = fallecidos ~ log_gizscore_z + longitud_z + latitud_z +
## entidad + departamento, data = datos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -19.4625 -2.4265 0.3614 2.2785 18.0528
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.9485 1.2992 5.348 1.81e-07 ***
## log_gizscore_z 5.1656 0.2911 17.743 < 2e-16 ***
## longitud_z -1.8381 1.4564 -1.262 0.20794
## latitud_z 4.3252 1.9363 2.234 0.02627 *
## entidadDEPTO -0.3881 0.9874 -0.393 0.69457
## entidadINVIAS-OTROS 1.4686 0.8369 1.755 0.08036 .
## departamentoARAUCA 5.1717 4.8081 1.076 0.28299
## departamentoATLÁNTICO -1.4869 3.8952 -0.382 0.70294
## departamentoBOLÍVAR -1.7877 3.3681 -0.531 0.59599
## departamentoBOYACÁ 6.2209 3.7620 1.654 0.09929 .
## departamentoCAQUETÁ 9.2110 4.9288 1.869 0.06266 .
## departamentoCAUCA 7.0701 3.4069 2.075 0.03885 *
## departamentoCESAR 2.5866 3.5599 0.727 0.46807
## departamentoCÓRDOBA 0.1602 2.1279 0.075 0.94003
## departamentoCUNDINAMARCA 6.2819 2.1121 2.974 0.00319 **
## departamentoHUILA 10.7194 3.6447 2.941 0.00354 **
## departamentoLA GUAJIRA 0.4588 5.5713 0.082 0.93443
## departamentoMAGDALENA -0.1578 4.2000 -0.038 0.97005
## departamentoMETA 9.5518 2.9753 3.210 0.00148 **
## departamentoNARIÑO 6.7851 5.0831 1.335 0.18299
## departamentoNORTE DE SANTANDER 4.9658 4.6114 1.077 0.28245
## departamentoQUINDÍO 0.5965 2.6158 0.228 0.81978
## departamentoRISARALDA 2.9088 3.7661 0.772 0.44054
## departamentoSANTANDER 2.7607 3.1067 0.889 0.37495
## departamentoSUCRE 0.6638 3.2453 0.205 0.83808
## departamentoTOLIMA 6.8763 2.5553 2.691 0.00754 **
## departamentoVALLE DEL CAUCA 4.1310 2.5231 1.637 0.10267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.8 on 288 degrees of freedom
## Multiple R-squared: 0.5789, Adjusted R-squared: 0.5409
## F-statistic: 15.23 on 26 and 288 DF, p-value: < 2.2e-16# Hacer predicciones
predicciones <- predict(modelo, newdata = datos)
# Comparar las predicciones con los valores reales
plot(datos$fallecidos, predicciones, main = "Comparación entre valores reales y predicciones", xlab = "Fallecidos reales", ylab = "Predicciones")
Diagnóstico de residuos:
Sirve para verificar si el modelo cumple con los supuestos clásicos (normalidad, homocedasticidad, independencia).
par(mfrow=c(2,2))
plot(modelo)
- Evaluación del ajuste:
Calcular el RMSE (error cuadrático medio) para ver qué tan bien predice el modelo
rmse(datos$fallecidos, predict(modelo))## [1] 4.58933Análisis del Modelo de Regresión Lineal Múltiple
Para cerrar el análisis, decidimos construir un modelo predictivo que nos permitiera estimar cuántas personas podrían fallecer en un siniestro vial, dependiendo de ciertas características. Para eso usamos una regresión lineal múltiple, que básicamente es una forma de predecir una variable (en este caso, fallecidos) a partir de varias otras.
¿Qué variables usamos?
log_gizscore_z: una medida transformada del riesgo vial.
longitud_z y latitud_z: que indican en qué lugar ocurrió el siniestro (coordenadas).
entidad: quién está a cargo de ese tramo de vía (ANI, INVÍAS, etc.).
departamento: en qué departamento ocurrió.
Estas variables se usan para ver qué tanto influyen en la cantidad de fallecidos y si nos ayudan a predecir con cierta precisión.
Resultados importantes del modelo
Intercepto (6.94): si todas las variables estuvieran en cero, el modelo empezaría prediciendo unos 7 fallecidos. Pero claro, esto solo sirve como punto de partida.
log_gizscore_z: fue la variable más importante. Tiene un efecto positivo y muy significativo. ¿Qué significa? Que a mayor riesgo (más alto el GIZScore), mayor cantidad de fallecidos. Es coherente y esperable.
latitud_z también fue significativa, lo que sugiere que la ubicación sí influye en la siniestralidad.
longitud_z no fue tan significativa, pero igual se incluye porque puede aportar contexto.
Entidades y departamentos: algunas categorías sí mostraron relación (como CAUCA, HUILA, META o TOLIMA), mientras que otras no aportan tanto individualmente.
¿Qué tan bueno fue el modelo?
R² ajustado = 0.539: eso quiere decir que el modelo explica el 53.9% de lo que pasa con los fallecidos. No es perfecto, pero es una muy buena base para este tipo de datos reales.
- Error estándar residual = 4.8: en promedio, el modelo se equivoca por unos 4-5 fallecidos. Está dentro de lo aceptable.
- MSE (error cuadrático medio) = 4.58: cuanto más bajo este valor, mejor. Este valor muestra que el modelo es bastante estable.
¿Y los gráficos qué muestran?
Gráfico 1 (real vs predicho): nos muestra qué tanto se parecen las predicciones a los datos reales. En general, se ve una tendencia positiva, lo que es bueno. Pero también vemos que el modelo a veces se equivoca un poco más en los casos extremos (cuando hay muy pocos o muchos fallecidos).
Gráficos de diagnóstico: sirven para verificar si el modelo cumple los supuestos básicos (como que los errores sean normales y estén bien repartidos). En este caso, no hay problemas graves, así que el modelo se puede considerar válido.
Conclusión final del modelo
Incluir este modelo fue una decisión inteligente y acertada. No solo nos ayudó a entender mejor cómo se relacionan las variables, sino que también nos permitió predecir la cantidad de fallecidos con una base estadística sólida. Aunque no es perfecto (ningún modelo lo es), sí aporta valor al análisis y refuerza lo que ya habíamos visto en la exploración anterior.
Sirve como una herramienta para:
Identificar factores clave que aumentan la siniestralidad.
Proponer intervenciones futuras según entidad o zona geográfica.
M- ostrar que sí es posible predecir, al menos en parte, cuántas vidas podrían estar en riesgo.
————————————————————————————————————————————————————-
Conclusiones Finales del Proyecto:
Exploración inicial y limpieza de datos
Analizamos las variables disponibles, transformamos algunas (como el logaritmo de fallecidos y GIZScore) para mejorar la interpretación, y estandarizamos variables numéricas. Esto nos permitió trabajar en una escala comparable y evitar sesgos.
Estadística descriptiva
Vimos que la mayoría de siniestros tienen entre 5 y 15 fallecidos, aunque hay casos extremos con hasta 50. Las entidades responsables más comunes fueron ANI e INVIAS-OTROS, y los departamentos con más casos fueron Cundinamarca y Valle del Cauca.
Análisis de correlaciones
Encontramos una relación moderada positiva entre el riesgo vial (GIZScore) y la cantidad de fallecidos. Otras variables como longitud y latitud mostraron relaciones más débiles pero útiles para el análisis espacial.
Gráficos multivariados
Usamos mapas, boxplots y mapas de calor para visualizar cómo varían los fallecidos según entidad, ubicación o nivel de riesgo. El heatmap mostró zonas del país con más concentración de siniestros fatales, especialmente entre el centro y el norte.
Tablas de contingencia y test de independencia
Al cruzar entidad y departamento, confirmamos que hay dependencia significativa entre ambos, lo que indica que ciertas entidades están más asociadas a ciertos territorios y esto podría influir en la cantidad de fallecidos.
Reducción de dimensiones (PCA y EFA) Con PCA, descubrimos que las variables se agrupan en dos dimensiones principales:
Una que representa la ubicación geográfica (longitud y latitud).
Otra que refleja la intensidad del siniestro (log de fallecidos y riesgo). Luego, EFA confirmó esta agrupación, mostrando una clara división entre lo espacial y lo mortal.
Modelo predictivo (Regresión lineal múltiple)
Creamos un modelo para predecir el número de fallecidos usando variables como el riesgo vial, ubicación y administración vial.
El modelo explica más del 50% de la variabilidad.
La variable más influyente fue el riesgo GIZScore.
El modelo no es perfecto, pero funciona bastante bien para hacer estimaciones razonables.
Conclusión general
Este proyecto nos permitió entender de forma más profunda dónde, cómo y bajo qué condiciones ocurren los siniestros viales más graves en Colombia. Combinamos análisis descriptivo, espacial, multivariado y predictivo para ofrecer una visión completa del problema.
La clave está en que los factores geográficos, el nivel de riesgo y la gestión vial sí influyen en la cantidad de fallecidos, y con modelos como este, se pueden generar alertas, tomar decisiones y, lo más importante, salvar vidas.