Resumen

La alta accidentalidad en la ciudad de Medellín, se ha convertido en un problema de salud pública, por esta razón, las entidades gubernamentales y el cuerpo de tránsito, requieren de campañas y metodologías fuertes con el fin de atacar dicha problemática. El objetivo principal de este trabajo, es realizar un aporte a dichas entidades gubernamentales, con el fin de atacar este problema social y de salud pública, donde tuvimos como base, los datos proporcionados por el portal MEData: MeData la cual, es una estrategia de datos de la ciudad de Medellín, que busca la apropiación, apertura y uso de los datos como herramienta de gobierno, de allí, se utilizó como ingrediente, la base de datos de incidentes viales en la ciudad de Medellín entre los años 2014 y 2021.

  1. Introducción

claro esta, que la estadística es una herramienta de suma importancia para la toma de decisiones y análisis que pueden contribuir a la solución de problemas actuales, por esta razón, este campo ha tomado un importante auge en la sociedad, las empresas requieren más estadísticos para realizar este tipo de tareas, elaboración, manejo y análisis de bases de datos. El objetivo de este trabajo es utilizar la Estadística y el desarrollo computacional para analizar los datos de incidentes georreferenciados de 2014 a 2021 de la ciudad de Medellín para predecir accidentes del año 2020 y 2021 con el fin de descubrir factores que se relacionan con los accidentes y aportar de forma significativa a las entidades públicas a tomar decisiones en temas de seguridad vial.

Enlaces de soporte

  • En el enlace se encuentra alojado el aplicativo web AccApp que permite visualizar la accidentalidad y pronosticar Enlace: Aplicacion web AccApp

  • En el enlace se encuentra alojado el Github principal que contiene todos los archivos y scripts utilizados para la construcción e implementación de este trabajo Enlace: Github principal

  • En el enlace se encuentra el video promocional de la aplicación AccApp que permite visualizar la accidentalidad y pronosticar Enlace: Video Promocional AccApp

  1. Metodología

Para la realización de este informe, se dispuso de una base de datos de incidentes viales en la ciudad de Medellín, entre los años 2014 y 2021, que se utilizaron como insumo principal para su respectivo análisis, con el fin de aplicar modelos predictivos y contribuir a la toma de decisiones.

  1. Descubrimiento y entendimiento de datos

La Alcaldía de Medellín en el portal GeoMedellín suministra las bases de datos anuales para incidentes de tránsito registrados por la Secretaria de Movilidad. Los datos se encuentran en MEData: MeData. Para el entendimiento del conjunto de datos se define el accidente de transito como: “evento, generalmente involuntario, generado al menos por un vehiculo en movimiento, que causa daños a personas y bienes involucrados en él que igualmente afecta la normal circulación de los vehiculos que se movilizan por la via o vias comprendidas en el lugar o dentro de la zona de influencia del hecho”

Las variables que se tienen en cuenta para el estudio son las siguientes:

Variable Descripción Tipo Detalle
X Coordenada X (longitud) de la ubicación del accidente float -
Y Coordenada Y (latitud) de la ubicación del accidente) float -
OBJECTID Identificación del registro integer -
RADICADO Identificación única del accidente ante la Secretaría de Movilidad string -
FECHA Fecha de la ocurrencia del accidente string -
HORA Hora aproximada del accidente datetime -
DIA Número día dentro del mes del accidente integer -
PERIODO Año de la ocurrencia del accidente integer -
CLASE Tipo de accidente String Clases: atropello, caída del ocupante, choque, incendio, volcamiento, otro
DIRECCION Dirección detallada de la ubicación del accidente String -
DIRECCION_ENC Dirección encasillada de la ubicación del accidente string Formato único de direcciones de la Alcaldía de Medellín
CBML Identificación única del lote más cercano a la ubicación del accidente string Acrónimo de comuna, barrio, manzana, lote
TIPO_GEOCOD Tipo de ubicación catastral string Formato de la oficina de catastro de la Alcaldía de Medellín
GRAVEDAD Gravedad del accidente string Clases: herido, muerto, solo daños
BARRIO Barrio donde ocurrió el accidente string 340 niveles dentro de la variable
COMUNA Comuna donde ocurrió el accidente string 85 niveles dentro de la variable
DISENO Tipo de diseño en donde ocurrió el accidente string Clases: ciclo ruta, glorieta, intersección, lote o predio, paso a nivel, paso elevado, paso inferior, pontón, puente, tramo de vía, túnel, vía peatonal
DIA_NOMBRE Nombre del día de la semana del accidente string -
MES Número de mes dentro del año del accidente integer -
MES_NOMBRE Nombre del mes del accidente string -
X_MAGNAMED - - -
Y_MAGNAMED - - -
LONGITUD Coordenada longitud de la ubicación del accidente float -
LATITUD Coordenada latitud de la ubicación del accidente float -

  1. Gestión de la base de datos

El conjunto de datos fue tratado y previamente depurado, ya que había presencia de datos vacíos y otras anomalías que afectaban el análisis, para esto, utilizamos librerías de R como Tidyverse para tratar la base.

Nuevas variables

Se crean nuevas variables para analizar de forma óptima los datos, en un ambiente más específico.

  • FESTIVIDAD: La cual esta caracterizada por “SI” y “NO”

Cabe resaltar, que para la obtención de la base de datos final, se realizó un procedimiento de depuración y manejo de datos detalladamente, el cual esta especificado en el archivo ManejoDatos.Rmd ubicado en el repositorio de Github, de forma resumida, se realizaron los siguientes procesos:

  • Corrección de vacíos identificados en la base de datos.
  • Arreglo de variables COMUNA y BARRIO, que están en CBML.
  • Eliminación de columnas repetidas.
  • Quitar tilde a nombre de barrio.
  • Separar variable LOCATION.
  • Renombrar columnas.
library(data.table);library(purrr);library(dplyr)        
library(plotly);library(tidyr);library(stringr)      
library(lubridate);library(sf);library(knitr);library(kableExtra)
medellin <- read.csv("incidentes_viales.csv", dec=",", header=T,sep=";", encoding = "UTF-8")
# Registros por día
medellin_dia <- medellin %>% group_by(DIA_SEMANA) %>% dplyr::summarize(total_accidentes = n())
# Ordenar de Lunes a Domingo
medellin_dia$DIA_SEMANA <- ordered(medellin_dia$DIA_SEMANA, 
                                   levels = c("LUNES","MARTES","MIERCOLES","JUEVES","VIERNES",
                                              "SABADO","DOMINGO"))
grafico_1 <- ggplot(data = medellin_dia, aes(x = DIA_SEMANA, y = total_accidentes)) +
  geom_bar(stat = "identity", alpha = 0.65, fill = "grey20", color = "black") + 
  geom_text(aes(y = total_accidentes, label = total_accidentes),
            position = position_dodge(width = 0.8), size = 2.8, vjust = 0.5, hjust = -0.1, col = "gray10") +
  xlab("Días") + 
  ylab("Total Accidentes") + 
  ggtitle("Total Accidentes por día de la semana") +
  theme_minimal() +
  ylim(c(0,40000)) +
  coord_flip()

  1. Validación y evaluación 

library(dummies)
library(stringr)
library(readxl)
library(sf)
library(dplyr)
library(lubridate)
library(ggplot2)
library(GGally)
library(car)
library(MLmetrics)
library(wordcloud)
library(gplots)
library(R.utils)
library(tm)
library(DescTools)
library(raster)
library(mclust)
library(rgdal)
library(raster)
library(geosphere)
library(NbClust)
library(factoextra)
library(vegan)
library(qpcR)
library(leaflet)
#Lectura de la Base de datos con la cual se trabajará en este proyecto
base_final <- read.csv("accident_final.csv", encoding="UTF-8")

Según el enunciado, tendremos en cuenta los años 2014 hasta 2017 para realizar entrenamiento y para validación los años 2018 y 2019

Modelo lineal 

#Modelo lineal
base_final$CLASE <- as.factor(as.character(base_final$CLASE))
datos_vl <- subset(base_final, (AÑO == '2018'))
base_final01 <- subset(base_final, (AÑO != '2018'))
base_final02 <- subset(base_final01, (AÑO != '2019'))
base_final03 <- subset(base_final02, (AÑO != '2020'))
datos_ml1 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   DISENO) %>% dplyr::count(name ="NRO_ACCID") 


ml1 <- lm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA+DISENO, data = datos_ml1)
promedio <- mean(datos_ml1$NRO_ACCID)
TSS <- sum((datos_ml1$NRO_ACCID - promedio)^2)
RSS <- RSS(ml1)
r2 <- 1-RSS/TSS
RSS2 <- anova(ml1)[4, 2]
r2 <- 1-RSS/TSS

Aquí podemos observar el MSE y el R^2, con el fin de determinar la potencia de dicho modelo.

Datos de entrenamiento año 2018 (Predicción y evaluación) 

ml1_data <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   DISENO) %>% dplyr::count(name = "NRO_ACCID") 
ml1_tr <- ml1_data[,-c(5)]
predicted <- round(predict(ml1, newdata=ml1_tr))
actual <- ml1_data$NRO_ACCID
ml1_mse <- MSE(predicted, actual) # MSE
ml1_mae <- MAE(predicted, actual) # MAE
ml1_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", ml1_mse, ml1_mae, ml1_r2)
## [1] "MSE: 116.151546, MAE: 5.882622, R2: 0.896806"

Datos de validación año 2018(Predicción y evaluación) 

ml1_2018 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                  DISENO) %>% dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(ml1, newdata=ml1_2018))
actual <- ml1_2018$NRO_ACCID
lm1_mse <- MSE(predicted, actual) # MSE
lm1_mae <- MAE(predicted, actual) # MAE
lm1_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm1_mse, lm1_mae, lm1_r2)
## [1] "MSE: 150.632552, MAE: 6.739856, R2: 0.728641"

En primera instancia, podemos observar que el \(R^2\) de los datos de validación para el año 2018, predice en un 72.86% y los datos de validación y entrenamiento del año 2018 tienen una variación en cuanto al MSE del 34.4%

Datos de validación para año 2019 (Predicción y evaluación) 

datos_vl02 <- subset(base_final, (AÑO == '2019'))
ml1_2019 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                  DISENO) %>% dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(ml1, newdata=ml1_2019))
actual <- ml1_2019$NRO_ACCID
ml1_mse <- MSE(predicted, actual) # MSE
ml1_mae <- MAE(predicted, actual) # MAE
ml1_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", ml1_mse, ml1_mae, ml1_r2)
## [1] "MSE: 157.439013, MAE: 7.161181, R2: 0.765329"

Según el \(R^2\) obtenido, con este modelo, para el año 2019 se predice en un 76.53%

Disminución de variables

Teniendo en cuenta solo las variables FESTIVIDAD y DIA_SEMANA 

datos_ml2 <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
lm2 <- lm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA, data = datos_ml2)

Predicción y evaluación para los datos de entrenamiento 

lm2_tr <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm2, newdata=lm2_tr))
actual <- lm2_tr$NRO_ACCID
lm2_mse <- MSE(predicted, actual) # MSE
lm2_mae <- MAE(predicted, actual) # MAE
lm2_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm2_mse, lm2_mae, lm2_r2)
## [1] "MSE: 687.801563, MAE: 20.398438, R2: 0.071931"

Predicción y evaluación para los datos de validación año 2018 

lm2_2018 <- datos_vl %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm2, newdata=lm2_2018))
actual <- lm2_2018$NRO_ACCID
lm2_mse <- MSE(predicted, actual) # MSE
lm2_mae <- MAE(predicted, actual) # MAE
lm2_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm2_mse, lm2_mae, lm2_r2)
## [1] "MSE: 629.906849, MAE: 18.904110, R2: 0.028546"

Se evidencia claramente, que se predice un porcentaje muy bajo con dicho modelo, únicamente se explica el 2.85% de los datos

Predicción y evaluación para los datos de validación año 2019 

lm2_2019 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm2, newdata=lm2_2019))
actual <- lm2_2019$NRO_ACCID
lm2_mse <- MSE(predicted, actual) # MSE
lm2_mae <- MAE(predicted, actual) # MAE
lm2_r2 <- R2_Score(predicted, actual) # R2
sprintf("MSE: %f, MAE: %f, R2: %f", lm2_mse, lm2_mae, lm2_r2)
## [1] "MSE: 778.819178, MAE: 22.430137, R2: 0.100687"

Para el año 2019 se predice en un 10%, porcentaje muy bajo, por lo cual, los resultados no fueron los esperados

Modelo lineal generalizado (Poisson) 

datos_lm3 <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
lm3 <- glm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA, family = "poisson", data = datos_lm3)

Predicción y evaluación para los datos de entrenamiento 

lm3_tr <- base_final03 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
lm3_tr_1 <- lm3_tr[,-4]
predicted <- round(predict(lm3, newdata=lm3_tr_1, type="response"))
actual <- lm3_tr$NRO_ACCID
lm3_mse <- MSE(predicted, actual) # MSE
lm3_mae <- MAE(predicted, actual) # MAE
lm3_r2 <- R2_Score(predicted, actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm3_mse, lm3_mae, lm3_r2)
## [1] "MSE: 687.635937, MAE: 20.395312, R2 Score: 0.072154"

Predicción y evaluación para los datos de validación en año 2018 

lm3_2018 <- datos_vl %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm3, newdata=lm3_2018, type="response")) 
actual <- lm3_2018$NRO_ACCID
lm3_mse <- MSE(predicted, actual) # MSE
lm3_mae <- MAE(predicted, actual) # MAE
lm3_r2 <- R2_Score(predicted, actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm3_mse, lm3_mae, lm3_r2)
## [1] "MSE: 629.742466, MAE: 18.898630, R2 Score: 0.028800"

En la etapa de entrenamiento, se obtiene un \(R^2\) de 7% y para la etapa de validación de 2.88%, lo cual indica que el modelo tiene una potencia muy baja, por lo que no se recomienda para el uso de predicciones.

Predicción y evaluación para los datos de validación en año 2019 

lm3_2019 <- datos_vl02 %>% group_by(FECHA,FESTIVIDAD, DIA_SEMANA) %>% 
  dplyr::count(name = "NRO_ACCID")
predicted <- round(predict(lm3, newdata=lm3_2019, type="response")) 
actual <- lm3_2019$NRO_ACCID
lm3_mse <- MSE(predicted, actual) # MSE
lm3_mae <- MAE(predicted, actual) # MAE
lm3_r2 <- R2_Score(predicted, actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm3_mse, lm3_mae, lm3_r2)
## [1] "MSE: 778.545205, MAE: 22.419178, R2 Score: 0.101004"

Según el \(R^2\) obtenido, claramente debemos adicionar otra variable, ya que la potencia fue muy baja

Modelo lineal generalizado + variable Clase 

datos_lm4 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   CLASE) %>% dplyr::count(name = "NRO_ACCID")
lm4 <- glm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA+CLASE, family = "poisson", 
           data = datos_lm4)

Predicción y evaluación para los datos 

datos_lm4_p <- datos_lm4[,-5]
y_train <- round(predict(lm4, newdata= datos_lm4_p, type="response"))
y_actual <- datos_lm4$NRO_ACCID
lm4_tmse <- MSE(y_train, y_actual)
lm4_tmae <-  MAE(y_train, y_actual)
lm4_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm4_tmse, lm4_tmae, lm4_r2)
## [1] "MSE: 98.608264, MAE: 5.762025, R2 Score: 0.886066"

Predicción y evaluación para los datos de validación en año 2018 

datos_lm4_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm4_v2 <- datos_lm4_v1[,-5]
y_train <- round(predict(lm4, newdata= datos_lm4_v2, type="response"))
y_actual <- datos_lm4_v1$NRO_ACCID
lm4_tmse <- MSE(y_train, y_actual)
lm4_tmae <-  MAE(y_train, y_actual)
lm4_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm4_tmse, lm4_tmae, lm4_r2)
## [1] "MSE: 88.650194, MAE: 5.470850, R2 Score: 0.895079"

Este modelo es mucho más acertivo, pues su potencia es mucho mayor, obteniendo un \(R^2\) de 88% para la etapa de entrenamiento y un \(R^2\) de 89% para validación en el año 2018, por ende, el modelo es útil para predecir accidentalidad en Medellín.

Predicción y evaluación para los datos de validación en año 2019 

datos_lm4_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm4_v2 <- datos_lm4_v1[,-5]
y_train <- round(predict(lm4, newdata= datos_lm4_v2, type="response"))
y_actual <- datos_lm4_v1$NRO_ACCID
lm4_tmse <- MSE(y_train, y_actual)
lm4_tmae <-  MAE(y_train, y_actual)
lm4_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm4_tmse, lm4_tmae, lm4_r2)
## [1] "MSE: 99.620330, MAE: 5.971978, R2 Score: 0.886356"

Se obtiene un \(R^2\) para los datos de validación de 88.6% evidenciando que el modelo tiene potencia para realizar predicciones sobre la accidentalidad en Medellín.

Modelo agregando variable DISEÑO 

datos_lm5 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   CLASE, DISENO) %>% dplyr::count(name = "NRO_ACCID")
lm5 <- glm(NRO_ACCID ~ FESTIVIDAD+DIA_SEMANA+CLASE+DISENO, family = "poisson", 
           data = datos_lm5)

Predicción y evaluación para los datos de entrenamiento 

datos_lm5_p <- datos_lm5[,-6]
y_train <- round(predict(lm5, newdata= datos_lm5_p, type="response"))
y_actual <- datos_lm5$NRO_ACCID
lm5_tmse01 <- MSE(y_train, y_actual)
lm5_tmae <-  MAE(y_train, y_actual)
lm5_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm5_tmse01, lm5_tmae, lm5_r2)
## [1] "MSE: 32.500873, MAE: 3.238339, R2 Score: 0.865619"

Predicción y evaluación para los datos de validación en año 2018 

datos_lm5_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE, DISENO) %>% dplyr::count(name = "NRO_ACCID")
datos_lm5_v2 <- datos_lm5_v1[,-6]
y_train <- round(predict(lm5, newdata= datos_lm5_v2, type="response"))
y_actual <- datos_lm5_v1$NRO_ACCID
lm5_tmse02 <- MSE(y_train, y_actual)
lm5_tmae <-  MAE(y_train, y_actual)
lm5_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm5_tmse02, lm5_tmae, lm5_r2)
## [1] "MSE: 27.591440, MAE: 3.236468, R2 Score: 0.819347"

Utilizando este modelo lineal generalizado poisson, teniendo en cuenta que se cuenta también con la variable DISEÑO, se obtiene un \(R^2\) en las etapas de entrenamiento y validación, de 86.5% y 82% respectivamente,lo que indica una potencia adecuada.

Predicción y evaluación para los datos de validación en año 2019 

datos_lm5_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                      CLASE, DISENO) %>% dplyr::count(name = "NRO_ACCID")
datos_lm5_v2 <- datos_lm5_v1[,-6]
y_train <- round(predict(lm5, newdata= datos_lm5_v2, type="response"))
y_actual <- datos_lm5_v1$NRO_ACCID
lm5_tmse02 <- MSE(y_train, y_actual)
lm5_tmae <-  MAE(y_train, y_actual)
lm5_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm5_tmse02, lm5_tmae, lm5_r2)
## [1] "MSE: 29.108963, MAE: 3.385074, R2 Score: 0.829171"

Para concluir, se decide trabajar con el modelo lineal generalizado, que contiene la variable CLASE, ya que según los resultados anteriores, son mucho más óptimos para la realización de predicciones

Predicción a nivel mensual

Modelo lineal generalizado Poisson agregando la variable Clase 

datos_lm7 <- base_final03 %>% group_by(FECHA, FESTIVIDAD, MES, 
                                   CLASE) %>% dplyr::count(name = "NRO_ACCID")
lm7 <- glm(NRO_ACCID ~ FESTIVIDAD+MES+CLASE, family = "poisson", 
           data = datos_lm7)

Predicción y evaluación para los datos de entrenamiento 

datos_lm7_p <- datos_lm7[,-5]
y_train <- round(predict(lm7, newdata= datos_lm7_p, type="response"))
y_actual <- datos_lm7$NRO_ACCID
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm7_tmse, lm7_tmae, lm7_r2)
## [1] "MSE: 96.534080, MAE: 5.724992, R2 Score: 0.888488"

Predicción y evaluación para los datos de validación en año 2018 

datos_lm7_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, MES, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm7_v2 <- datos_lm7_v1[,-5]
y_train <- round(predict(lm7, newdata= datos_lm7_v2, type="response"))
y_actual <- datos_lm7_v1$NRO_ACCID
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm7_tmse, lm7_tmae, lm7_r2)
## [1] "MSE: 94.413068, MAE: 5.552602, R2 Score: 0.888105"

Se obtiene un \(R^2\) en la etapa de entrenamiento y validación del 88.84% y de 88.81% respectivamente, lo cual indica una buena potencia para predecir accidentalidad en Medellín

Predicción y evaluación para los datos de validación en año 2019 

datos_lm7_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, MES, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm7_v2 <- datos_lm7_v1[,-5]
y_train <- round(predict(lm7, newdata= datos_lm7_v2, type="response"))
y_actual <- datos_lm7_v1$NRO_ACCID
lm7_tmse <- MSE(y_train, y_actual)
lm7_tmae <-  MAE(y_train, y_actual)
lm7_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm7_tmse, lm7_tmae, lm7_r2)
## [1] "MSE: 97.898352, MAE: 5.923626, R2 Score: 0.888320"

En este modelo la diferencia del MSE entre los datos de entrenamiento y validaci?n del 2019 es del 1.3%, que fue menor al valor obtenido con la validaci?n del 2018 (2.19%), que indica claramente que no hay problemas de sobreentrenamiento. Adem?s se puede observar que el \(R^2\) de los datos de validaci?n para el AÑO 2019 predice un 88.83%, evidenciando as? que este modelo lineal generalizado con la adici?n de la variable Clase es un buen candidato para predecir la accidentalidad en Medell?n.

Se concluye que el modelo tiene buena potencia para predecir accidentes en la ciudad de Medellin.

  1. Predicciones 

datos_real_pred <- read_excel("Base_2020_real_predict.xlsx")

Predicciones para el año 2020

Encabezado de datos con predicciones diarias para el año 2020 

Base_prediccion <- read.csv("prediccion.csv", sep = ",", encoding = "UTF-8")
Base_prediccion <- Base_prediccion[,-1]
Base_prediccion_2020 <- subset(Base_prediccion, (AÑO != '2020'))
Base_prediccion_2020$FECHA <- as.Date(Base_prediccion_2020$FECHA)
Base_prediccion_2020$CLASE <- as.factor(Base_prediccion_2020$CLASE)
Base_prediccion_2020$DIA_SEMANA <- as.factor(Base_prediccion_2020$DIA_SEMANA)
Base_prediccion_2020$AÑO <- as.integer(Base_prediccion_2020$AÑO)
Base_prediccion_2020$FESTIVIDAD <- as.factor(Base_prediccion_2020$FESTIVIDAD)
prediccion_2020 <- predict(object = lm4, newdata = Base_prediccion_2020,
                          type = "response")
prediccion_diaria2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020,0))
diario_20_02 <- prediccion_diaria2020 %>%
  group_by(FECHA, DIA_SEMANA, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=NRO_ACCID)
head(diario_20_02, 6)

Encabezado de datos con predicciones mensual para el año 2020 

Base_prediccion <- read.csv("prediccion.csv", sep = ",", encoding = "UTF-8")
Base_prediccion_2020 <- subset(Base_prediccion, (AÑO != '2020'))
Base_prediccion_2020$FECHA <- as.Date(Base_prediccion_2020$FECHA)
Base_prediccion_2020$CLASE <- as.factor(Base_prediccion_2020$CLASE)
Base_prediccion_2020$DIA_SEMANA <- as.factor(Base_prediccion_2020$DIA_SEMANA)
Base_prediccion_2020$AÑO <- as.integer(Base_prediccion_2020$AÑO)
Base_prediccion_2020$FESTIVIDAD <- as.factor(Base_prediccion_2020$FESTIVIDAD)
prediccion_2020 <- predict(object = lm4, newdata = Base_prediccion_2020,
                          type = "response")
prediccion_mensual2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020,0))
#Se borraron columnas no necesarias
prediccion_mensual2020 <-  prediccion_mensual2020[,c(-1,-2,-3,-5,-7)]
#Agrupamiento por mes 2020
#Se agrup? por mes y si fue en d?a festivo o no
mensual_20 <- prediccion_mensual2020 %>% group_by(CLASE, MES, NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
mensual_20 <- mutate(mensual_20, NRO_ACCID_TOTAL=NRO_ACCID*total)
mensual_20_02 <- mensual_20 %>%
  group_by(MES, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))
head(mensual_20_02, 6)

Predicciones para el año 2021

Encabezado de datos con predicciones diarias para el año 2021 

Base_prediccion_2021 <- subset(Base_prediccion, (AÑO != '2021'))
Base_prediccion_2021$FECHA <- as.Date(Base_prediccion_2021$FECHA)
Base_prediccion_2021$CLASE <- as.factor(Base_prediccion_2021$CLASE)
Base_prediccion_2021$DIA_SEMANA <- as.factor(Base_prediccion_2021$DIA_SEMANA)
Base_prediccion_2021$AÑO <- as.integer(Base_prediccion_2021$AÑO)
Base_prediccion_2021$FESTIVIDAD <- as.factor(Base_prediccion_2021$FESTIVIDAD)
prediccion_2021 <- predict(object = lm4, newdata = Base_prediccion_2021,
                          type = "response")
prediccion_diaria2021 <- Base_prediccion_2021 %>% 
  mutate(NRO_ACCID = round(prediccion_2021,0))
diario_21_02 <- prediccion_diaria2021 %>%
  group_by(FECHA, DIA_SEMANA, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=NRO_ACCID)
head(diario_21_02, 6)

Encabezado de datos con predicciones mensual para el año 2021 

Base_prediccion <- read.csv("prediccion.csv", sep = ",", encoding = "UTF-8")
Base_prediccion_2020 <- subset(Base_prediccion, (AÑO != '2021'))
Base_prediccion_2020$FECHA <- as.Date(Base_prediccion_2020$FECHA)
Base_prediccion_2020$CLASE <- as.factor(Base_prediccion_2020$CLASE)
Base_prediccion_2020$DIA_SEMANA <- as.factor(Base_prediccion_2020$DIA_SEMANA)
Base_prediccion_2020$AÑO <- as.integer(Base_prediccion_2020$AÑO)
Base_prediccion_2020$FESTIVIDAD <- as.factor(Base_prediccion_2020$FESTIVIDAD)
prediccion_2020 <- predict(object = lm4, newdata = Base_prediccion_2020,
                          type = "response")
prediccion_mensual2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020,0))
#Se borraron columnas no necesarias
prediccion_mensual2020 <-  prediccion_mensual2020[,c(-1,-2,-3,-5,-7)]
#Agrupamiento por mes 2020
#Se agrup? por mes y si fue en d?a festivo o no
mensual_20 <- prediccion_mensual2020 %>% group_by(CLASE, MES, NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
mensual_20 <- mutate(mensual_20, NRO_ACCID_TOTAL=NRO_ACCID*total)
mensual_20_02 <- mensual_20 %>%
  group_by(MES, CLASE, FESTIVIDAD) %>%
  summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))
head(mensual_20_02, 6)

  1. Clusters Barrios

Mapa de calor de accidentalidad en la ciudad de Medellín entre los años 2014 y 2019 (Histórico) 

#lectura de .csv y .shp
catastral <- read.csv("Limite_Barrio_Vereda_Catastral.csv", encoding="UTF-8")
catastro <- read_sf("Limite_Barrio_Vereda_Catastral.shp")
barrio_vereda <- read.csv("Barrio_Vereda_2014.csv", encoding="UTF-8")
#Mapa para todos los barrios, usando 'innerjoin' con el .shp de Limite_Barrio_Vereda_Catastral
Unido <- inner_join(catastral, base_final, by = c("COMUNA" = "NUMCOMUNA"))
nueva_base <- Unido %>% filter(AÑO >= 2014 & AÑO <= 2019) %>% 
  group_by(CODIGO) %>%
  dplyr::summarise(accidentes = n()) %>%
  dplyr::ungroup()
#Se realiz? la conversi?n de CODIGO a formato num?rico
catastro$CODIGO <- as.numeric(as.character(catastro$CODIGO))
#Se utiliz? 'inner join' para unir dos bases y para luego generar mapa
mapa <- inner_join(catastro, nueva_base, by = c("CODIGO" = "CODIGO"))
mypal <- colorNumeric(palette = c("#000000","#280100","#3D0201","#630201","#890100","#B00100","#DD0100","#F50201",
                                   "#FF5F5E","#FF7A79","#FF9796","#FEB1B0","#FDC9C8", "#FFE5E4"), domain = mapa$accidentes, reverse = T)
# Creaci?n del mapa
leaflet() %>% addPolygons(data = mapa, color = "#0A0A0A", opacity = 0.6, weight = 1, fillColor = ~mypal(mapa$accidentes),
                          fillOpacity = 0.6, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "black", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", mapa$NOMBRE_BAR, "<br>", "Accidentes: ", mapa$accidentes, "<br>")) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", pal = mypal, values = mapa$accidentes, title = "Accidentes", opacity = 0.6)
medellin_comuna <- base_final %>% filter(AÑO >= 2014 & AÑO <= 2019) %>% 
  group_by(COMUNA) %>% 
  dplyr::summarize(accidentes = n())
ggplot(data = medellin_comuna, aes(x = reorder(COMUNA,+accidentes), y = accidentes)) +
  geom_bar(stat = "identity",  fill = "red", color = "black", alpha = 0.6) +
  geom_text(aes(y = accidentes, label = accidentes),
                position = position_dodge(width = 0.7), size = 3.5, vjust = 0.5, hjust = -0.1, col = "black") +
  xlab("Comuna") + 
  ylab("Total de Accidentes") +
  ggtitle("Total de Accidentes por Comuna") +
  ylim(c(0,50000)) +
  theme_minimal() +
  coord_flip()

Se puede evidenciar de los gráficos anteriores, que la comuna la candelaria, lidera significativamente en índices de accidentalidad en la ciudad de Medellín, siguiéndole laureles, posteriormente, Palmitas es la que manifiesta menor cantidad de accidentalidad. 

geo.dist = function(df) {
  require(geosphere)
  d <- function(i,z){         # z[1:2] contain long, lat
    dist <- rep(0,nrow(z))
    dist[i:nrow(z)] <- distHaversine(z[i:nrow(z),1:2],z[i,1:2])
    return(dist)
  }
  dm <- do.call(cbind,lapply(1:nrow(df),d,df))
  return(as.dist(dm))
}
base_final03$LATITUD <- as.numeric(as.character(base_final03$LATITUD))
base_final03$LONGITUD<- as.numeric(as.character(base_final03$LONGITUD))
df <- data.frame(long = base_final03$LONGITUD, lat = base_final03$LATITUD, barrios = base_final03$BARRIO)
library(geosphere)
library(shapefiles)
library(foreign)
#Se creó con la función 'geo.dist', una matriz de distancias
df1 <- df[1:1000, ]
d <- geo.dist(df1)
hc <- hclust(d)
df1$clust <- cutree(hc, k = 6)
head(df1,10)

Mapa según Latitud y Longitud

s <- shapefile("Limite_Catastral_de__Comunas_y_Corregimientos.shp")
map.df1 <- (s)
ggplot(map.df1)+
  geom_path(aes(x=long, y=lat, group=group))+
  geom_point(data=df1, aes(x=long, y=lat, color=factor(clust)), size=4)+
  scale_color_discrete("Cluster")+
  ggtitle("Cluster Medellín")+
  coord_fixed()

El mapa anterior fue realizado según medidas geoespaciales, no obstante, la agrupación fue realizada sin ningún método de elección de un k óptimo, para su solución, utilizamos el algoritmo “K means” con el fin de allar tal K que optimice el agrupamiento.

Clusterización por gravedad de accidente

#Numero de accidentes por Barrio
datos_cluster <- base_final03 %>% group_by(BARRIO) %>% dplyr::count(name = "TOTAL_ACCIDENTES")
#Número de accidentes por barrio, seg?n gravedad almacenado en 'df'
df <- as.matrix(table(base_final03$BARRIO, base_final03$GRAVEDAD))
df <- data.frame(Con_heridos = df[,1], Con_muertos = df[,2], Solo_danos = df[,3])
#Escalamiento y centrado de la base de datos.
scaled_data = as.matrix(scale(df))
head(scaled_data, 10)
##                          Con_heridos Con_muertos  Solo_danos
## Aguas Frias               -0.7826375 -0.62264869 -0.64131796
## Aldea Pablo VI            -0.8244313 -0.34908854 -0.67245041
## Alejandria                -0.5945652 -0.34908854  0.13388001
## Alejandro Echavarria       0.5373509 -0.07552839 -0.16810474
## Alfonso Lopez              0.8995640  0.19803176 -0.02800872
## Altamira                   0.1925518 -0.62264869 -0.11829282
## Altavista                 -0.4099758 -0.34908854 -0.38603188
## Altavista Sector Central  -0.4587353 -0.07552839 -0.59150604
## Altos del Poblado         -0.5667027 -0.62264869 -0.35801268
## Andalucia                 -0.5492886 -0.07552839 -0.55103386
kmm = kmeans(scaled_data, 5, nstart = 50, iter.max = 15 )

Búsqueda de un K óptimo

Método del codo 

#Se fij? una semilla y se realiz? el calculo y se grafico el WSS(total within - cluster sum of square) para k = 2 hasta k = 10
set.seed(2021021)
k.max <- 10
datos <- scaled_data
wss <- sapply(2:k.max, 
              function(k){kmeans(datos, k, nstart = 50, iter.max = 15 )$tot.withinss})
plot(2:k.max, wss, 
     type = "b", pch = 19, frame = FALSE,
     xlab = "Número de Clusters (k)",
     ylab = "WSS Total", 
     main = "Método del Codo", col="forestgreen")

kn <- kmeans(datos, 5)

Según el método anterior, los candidatos más marcados para ser K óptimo serían k=4 y k=5, dado que presentan un delta mucho más bajo en sus pendientes en comparación con las demás.

Método de la silueta

Método Gap statistic 

set.seed(1152)
fviz_nbclust(scaled_data, kmeans, nstart = 25,  method = "gap_stat", nboot = 50)+
  labs(subtitle = "Gap statistic method")

En conclusión de todos los métodos, tener \(k = 4\) es el más común para tener en cuenta como K óptimo para la realización de la clusterización, ahora tenemos en cuenta lo siguiente: 

kmm = kmeans(scaled_data, 4, nstart = 50, iter.max = 15 )
df_clust <- data.frame(Con_heridos = df[,1], Con_muertos = df[,2], Solo_danos = df[,3], kmm$cluster)
head(df_clust, 10)

Tenemos en cuenta la siguiente clasificación:

  • Grupo 1: Accidentalidad Moderada
  • Grupo 2: Accidentalidad Baja
  • Grupo 3: Accidentalidad Alta
  • Grupo 4: Accidentalidad Media-Alta

Información para cada grupo:

GRUPO 1: Accidentalidad Moderada 

#Accidentalidad Moderada
dfclust_clust1 <- df_clust[df_clust$kmm.cluster == 1, ]
dfclust_clust1$total <- rowSums(dfclust_clust1[,1:3])
sum(dfclust_clust1$Con_heridos)
## [1] 13145
sum(dfclust_clust1$Con_muertos)
## [1] 111
sum(dfclust_clust1$Solo_danos)
## [1] 14639
sum(dfclust_clust1$total)
## [1] 27895

Donde obtenemos una cantidad total de accidentes de 27895, de los cuales se tienen resultados a nivel de gravedad de 14639 accidentes de solo daños, 111 accidentes con presencia de muertos y 13145 accidentes con heridos.

GRUPO 2: Accidentalidad Baja 

#Accidentalidad Baja
dfclust_clust2 <- df_clust[df_clust$kmm.cluster == 2, ]
dfclust_clust2$total <- rowSums(dfclust_clust2[,1:3])
sum(dfclust_clust2$Con_heridos)
## [1] 17817
sum(dfclust_clust2$Con_muertos)
## [1] 106
sum(dfclust_clust2$Solo_danos)
## [1] 10699
sum(dfclust_clust2$total)
## [1] 28622

Donde obtenemos una cantidad total de accidentes de 28622, de los cuales se tienen resultados a nivel de gravedad de 10699 accidentes de solo daños, 106 accidentes con presencia de muertos y 17817 accidentes con heridos.

GRUPO 3: Accidentalidad Alta 

#Accidentalidad Alta
dfclust_clust3 <- df_clust[df_clust$kmm.cluster == 3, ]
dfclust_clust3$total <- rowSums(dfclust_clust3[,1:3])
sum(dfclust_clust3$Con_heridos)
## [1] 17147
sum(dfclust_clust3$Con_muertos)
## [1] 237
sum(dfclust_clust3$Solo_danos)
## [1] 18863
sum(dfclust_clust3$total)
## [1] 36247

Donde obtenemos una cantidad total de accidentes de 36247, de los cuales se tienen resultados a nivel de gravedad de 18863 accidentes de solo daños, 237 accidentes con presencia de muertos y 17147 accidentes con heridos.

GRUPO 4: Accidentalidad Media-Alta 

#Accidentalidad Media-Alta
dfclust_clust4 <- df_clust[df_clust$kmm.cluster == 4, ]
dfclust_clust4$total <- rowSums(dfclust_clust4[,1:3])
sum(dfclust_clust4$Con_heridos)
## [1] 32887
sum(dfclust_clust4$Con_muertos)
## [1] 222
sum(dfclust_clust4$Solo_danos)
## [1] 21732
sum(dfclust_clust4$total)
## [1] 54841

Donde obtenemos una cantidad total de accidentes de 54841, de los cuales se tienen resultados a nivel de gravedad de 21732 accidentes de solo daños, 222 accidentes con presencia de muertos y 32887 accidentes con heridos.

Luego, se obtiene el siguiente mapa clúster según gravedad del accidente con K= 4 

#Se vuelve a utlizar catastro para este mapa
#Se import? el archivo .xlsx basemapa
basemapa <- read_excel("basemapa.xlsx")
base_mapa <- data_frame(basemapa)
catastro$CODIGO <- as.numeric(as.character(catastro$CODIGO))
base_mapa$Codigo <- as.numeric(as.character(base_mapa$Codigo))
#Se utiliz? 'inner join' de nuevo para unir dos bases y para as? luego generar mapa
mapa02 <- inner_join(catastro, base_mapa, by = c("CODIGO" = "Codigo"))
colorgrupos <- c("#CCFF00", "#00FF66", "#A52A2A", "#FF7F00")
    mapa02$colores <- ifelse(mapa02$kmm.cluster == "1", "#CCFF00",
                             ifelse(mapa02$kmm.cluster == "2", "#00FF66",
                                    ifelse(mapa02$kmm.cluster == "3", "#A52A2A",
                                           ifelse(mapa02$kmm.cluster == "4", "#FF7F00",0))))
leaflet() %>% addPolygons(data = mapa02, opacity = 0.4, color = "#545454",weight = 1, fillColor = mapa02$colores,
                          fillOpacity = 0.4, label = ~NOMBRE_BAR,
                          highlightOptions = highlightOptions(color = "#262626", weight = 3, bringToFront = T, opacity = 1),
                          popup = paste("Barrio: ", mapa02$NOMBRE_BAR, "<br>", "Grupo: ", mapa02$kmm.cluster, "<br>", "Número de Accidentes con heridos: ", mapa02$Con_heridos, "<br>", "Número de Accidentes con muertos: ", mapa02$Con_muertos, "<br>", "Número de Accidentes con solo dAÑOs: ", mapa02$Solo_danos)) %>%
  addProviderTiles(providers$OpenStreetMap) %>%
  addLegend(position = "bottomright", colors = colorgrupos, labels = c("Grupo 1: Accidentalidad Moderada", "Grupo 2: Accidentalidad Baja", "Grupo 3: Accidentalidad Alta", "Grupo 4: Accidentalidad Media-Alta"))

En el anterior mapa, podemos analizar como se presentan los niveles de gravedad enunciados alrededor del mapa de Medellín, donde se puede evidenciar Accidentalidad alta más que todo cerca al centro de Medellín, aunque se da atípicamente accidentalidad alta en las afueras de la ciudad, más específicamente en la Cabecera de San Antonio de Prado. Estos resultados, pueden impulsar a las entidades gubernamentales, a realizar su debido proceso, ya sea con más señalización o implementar medidas que aporten significativamente a la reducción de accidentes de éstas zonas, más que todo en La candelaria, que fue la zona que más accidentes ha reportado en términos históricos.

8. Bibliografia

  • James, G., Witten, D., Hastie, T., & Tibshirani, R. (2013). An introduction to statistical learning (Vol. 112, p. 18). New York: springer.
  • Ospina J. (2019) Introducción a la analítica predictiva, Vehículos registrados en el RUNT. Material de clase especialización en Analítica Unal.
  • Concejo de Medellín. La Alta accidentalidad en Medellín se convirtió en un problema de Salud Pública. Boletin 070. 2018
  • Página Shiny Rstudio. https://shiny.rstudio.com/ +Dudas R. https://es.stackoverflow.com/

–