Análisis, predicción y pronóstico de la accidentalidad en la ciudad de Medellín (2014-2021)

1. Introducción

La accidentalidad vial en la ciudad de Medellín, es un problema estructural que limita el crecimiento económico y social, que debe ser atacado con politicas y estrategias de corto, mediano y largo plazo, en donde el gobierno, los entes privados y la academia, deben aunar esfuerzos por encontrar soluciones que beneficien a la sociedad en general.

En el presente reporte, se busca contribuir a esas estrategías que permitan dar ideas o soluciones, para afrontar la problemática de la accidentalidad vial en la ciudad de Medellín, a través del análisis y la predicción de la accidentalidad en la ciudad, utilizando como insumo los datos proporcionados por la Alcaldía de Medellín en el portal MeData, en donde se reportan los incidentes viales en la ciudad entre los años 2014 y 2020.

---

2. Datos

La base de datos inicial, fue consultada el día 18/11/2021, en formato .CSV. Se usaron las herramientas de Excel y R para el procesamiento de la base de datos. Inicialmente, esta contenía 270.765 observaciones y 17 variables. Directamente en el excel, antes de iniciar con el procesamiento en R, se crean dos variables nuevas, las cuales son; DIA_SEMANA y HORA. Luego, se generan 19 variables en total. Una vista previa de la primeras observaciones de la base de datos inicial es la siguiente:

head(base,n=5)

#coercionando fecha accidente a tipo date
base$FECHA_ACCIDENTE <- as.Date(base$FECHA_ACCIDENTE, format="%d/%m/%Y") 

#class(base$FECHA_ACCIDENTE)


#colSums(is.na(base)|base=="") #cuenta los vacios

#summary(base)
#str(base)

Lo siguiente, será analizar detalladamente variable por variable en busca de observaciones faltantes, erróneas, particulares, entre otras y con esto se quiere eliminar las variables que no se consideren pertinentes para este análisis.

AÑO: año de ocurrencia del incidente. (2014 hasta 2016)

CBML: código catastral que corresponde al código comuna, barrio, manzana, lote catastral de un predio. Hallazgo: 18.156 vacíos, se tienen 962 registros más con caracteres extraños como: AUC1, AUC2, Inst_14, Inst_16, Inst_18, Inst_19, Sin Inf SN01, para un total de 19.118 registros mal estructurados o vacíos.

CLASE_ACCIDENTE: clasificación del IPAT sobre la clase de accidente de transito: choque, atropello, volcamiento, caida de ocupante, incendio, u otro (que no corresponde a las anteriores 5 clasificaciones, p. ej: sumersión). Hallazgo: 6 datos vacíos, se cambiarán por “otro”, pues no fueron reportados o se perdieron al momento de tomar la medición.

DISEÑO: en la variable DISEÑO que es el sitio donde ocurrió el accidente( Cicloruta, Glorieta, Intersección, Lote o Predio, Paso a Nivel, Paso Elevado, Paso Inferior, Pontón, Puente, Tramo de via, Tunel, Via peatonal). Hallazgo: se encuentran 1.148 vacíos, se reemplazarán los vacíos por “otro”.

BARRIO, COMUNA, NUMCOMUNA:

BARRIO: barrio de ocurrencia del incidente vial.

Hallazgos: 19.006 vacíos. Además de los vacíos, se tienen 1.822 registros adicionales con carácteres como: números entre 0 y 9.086, AUC1, AUC2, Inst, Sin Inf, Sin nombre. Para un total de 20.828 registros mal estructurados o vacíos.

COMUNA: denominación con la cual se identifica cada Comuna o Corregimiento.

Hallazgos: 12.798 vacíos. Se tienen 7.064 registros adicionales con carácteres como: No Georef, 0, In, AU, Sin Inf, SN. Para un total de 19862 registros mal estructurados o vacíos.

NUMCOMUNA: número de la comununa en la que ocurrió incidente vial. Hallazgos: se tienen 20.116 registros adicionales con carácteres como: AU, In, Sin Inf, SN.

LOCATION: fuente de información con la cual se realizó la geocodificación. Contiene la latitud y longitud. Posteriormente será separada en dos variables.

X: coordenada X en metros del accidente, en sistema de coordenadas MAGNA Medellin Local

Y: coordenada Y en metros del accidente, en sistema de coordenadas MAGNA Medellin Local

NRO_RADICADO: consecutivo que asigna UNE, según el orden de llegada de los expedientes para su diligenciamiento.

MES: mes de ocurrencia del incidente vial. Esta variable no se modifica.

GRAVEDAD_ACCIDENTE: clasificación del IPAT - Informe Policial de Accidentes de Tránsito, sobre la gravedad del accidente, corresponde al resultado más grave presentado en el accidente. Daños materiales “Sólo daños”, accidente con heridos “Herido”, accidente con muertos “Muerto”. No indica cantidad. Hallazgo: esta variable contenía algunos carácteres extraños, estos se modifican y se restauran.

FECHA_ACCIDENTES: fecha de los accidente (formato YYYY-MM-DD hh:mi:ss), proviene del IPAT - Informe Policial de accidentes de Tránsito

FECHA_ACCIDENTE: fecha del accidente, proviene del IPAT - Informe Policial de accidente de Tránsito. Esta variable posteriormente se elimina debido a que porporciona menos información que la variable FECHA_ACCIDENTES.

EXPEDIENTE: consecutivo que asigna UNE, según el orden de llegada de los expedientes para su diligenciamiento. Esta variable posteriormente se elimina.

DIRECCION ENCASILLADA: dirección encasillada que entrega el geocodificador. Esta variable se elimina.

DIRECCION: dirección donde ocurrió el incidente. Esta variable no se modifica.

NRO_RADICADO: consecutivo que asigna UNE, según el orden de llegada de los expedientes para su diligenciamiento.

2.1. Integración con los datos de Geo medellín y depuración de vacíos

Con la variable CBML y creando una nueva variable con los 4 primeros dígitos de esta, podemos obtener la comuna y el barrio de las filas en donde estos datos se encuentren vacíos. Para eso usaremos un archivo proporcionado por GEO medellín, donde se registran los nombres de barrios y comunas con su código (CBML). Lo siguiente, será cruzar los datos de incidentes viales y los datos de Geo medellín teniendo como variable común el código y CBML.

Una vez realizado este proceso, se eliminan 23.386 observaciones, de las que no se pudo recuperar datos.

Luego se procedió a crear/eliminar/modificar, las variables que no se consideran relevantes para el análisis.

Se eliminan las variables:

BARRIO,COMUNA,DIRECCION.ENCASILLADA,CBML,CB,JCB,FECHA_ACCIDENTES,EXPEDIENTE,DISEÑO.
Algunas de estas variables se encontraban duplicadas.
* Se quitan tildes
* Se separar la variable LOCATION para obtener las variables LATITUD y LONGITUD.
* Se renombran algunas columnas.
* Se crea la variable SEMANA.

2.2. Fechas especiales

Para las fechas especiales se crean dos nuevas variables; FESTIVIDAD y TIPO_FESTIVIDAD. Estas variables provienen de una base de datos externa que se adiciona a la base de análisis y abarca los días feriados en colombia desde 2014 hasta 2021.

FESTIVIDAD: contiene dos etiquetas (SI/NO). SI: cuando hay festividad para ese día. NO: cuando no hay festividad para ese día,

TIPO_FESTIVIDAD: contiene seis tipos de festividad:
* FESTIVO: día feriado.
* NAVIDAD: 24,25 y 31 de diciembre.
* SEM_SANTA: toda la semana santa, desde el lunes hasta el domingo.
* BRUJAS: 31 de octubre.
* MADRES: el día d emadres designado para el año respectivo.
* A_NUEVO: primero de enero de cada año.

Una vez agregada la festividad, la nueva base de datos queda conformada por 246.417 observaciones y 19 variables.

head(base_prueba, n=6)

base_final <- base_prueba

---

3. Análisis descriptivo

3.1. Promedio accidentes mensuales por año

aggregate(numero_de_accidentes~ano*mes, data = acc, FUN = mean) %>%
  plot_ly(x = ~mes,
          y = ~numero_de_accidentes, type = "scatter", mode = "lines",
          split = ~ano, line = list(width = 1.5)) %>%
  layout(title = 'Promedio accidentes mensuales mensual por año',
         xaxis = list(title = "Mes"),
         yaxis = list(title = "Número de accidentes"))

3.2. Registros de accidentalidad diarios 2014-2020

plot_ly(data = acc, x = ~FECHA, y = ~numero_de_accidentes,
        type = "scatter", mode = "lines", split = ~ano,
        line = list(width = 1)) %>%
  layout(title = 'Registros de accidentalidad diarios 2014-2020',
        xaxis = list(title = "Día"),
        yaxis = list(title = "Número de accidentes"))

De las gráficas anteriores, se puede observar que:

• Después de febrero el promedio de accidentes en el 2020 disminuyen a gran medida, a tal punto que en mayo solo hay en promedio 26 accidentes.
• El día con menos accidentes es el 23 de mayo del 2020, con solo 9 accidentes.
• El día con más accidentes es el 4 de agosto de 2017, con 213 accidentes.
• En el año 2014 no hay registrados datos de accidentes antes del 4 de julio, lo cual no significa que en estas fechas no ocurrieran accidentes. -Al igual que en 2014, en el año 2020 no se han registrado accidentes en la base de datos abiertos de Medellín, después del 31 de agosto hasta la fecha actual.

3.3. Accidentes por día de la semana

ggarrange(gdia,gsem)

El día que presenta mayor cantidad de personas accidentadas, es el lunes seguido del domingo, con una diferencia de 660 personas. Tras estos dos sigue el martes, mientras que los días miércoles, jueves, viernes y sábado, presentan accidentes más cercamos, con la mayor diferencia entre jueves y sábado de 300 accidentes. El sábado, es el día de la semana con menor cantidad de accidentes.

3.4. Accidentes por mes

ggplot(data = mes, aes(fill = MES, x = MES, y = numero_de_accidentes)) + 
  geom_bar(stat = "identity") +
  scale_fill_viridis_d( option = "D") +
  geom_text(aes(y = numero_de_accidentes, label = numero_de_accidentes), vjust = -0.5) +
  xlab("Mes") +
  ylab("Número de accidentes") +
  ylim(c(0,26000)) +
  ggtitle("Número de accidentes por mes")+
  theme(legend.position = "none")

3.5. Accidentes por año

ggplot(data = año, aes(fill = Var1, x = Var1, y = Freq)) +
  geom_bar(stat = "identity") +
  geom_text(aes(y = Freq, label = Freq), vjust = -0.5) +
  scale_fill_viridis_d( option = "H") + 
  xlab("Año") +
  ylab("Número de accidentes") +
  ylim(c(0,45000)) +
  ggtitle("Número de accidentes por año")+
  theme(legend.position = "none")

Se debe tener en cuenta, que en los años 2014 y 2020 faltan algunos meses, por ello se presenta una cantidad baja de accidentes, si tener en cuenta ambos años, año con menor cantidad de accidentes es 2018 y 2016 fue el año con más accidentes sólo por unos cientos más que 2015 y 2017.

3.6. Festividad

ggarrange(gfes1,gfes2)

3.7. Acidentes por comuna

ggplot(data = acc_comuna, aes( fill = COMUNA, x = reorder(COMUNA,+accidentes), y = accidentes)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis_d( option = "B") +
  geom_text(aes(y = accidentes, label = accidentes), hjust = -0.1) +
  xlab("Comuna") + 
  ylab("Total accidentes") +
  ggtitle("Total accidentes por comuna 2014-2018") +
  ylim(c(0,60000)) +
  coord_flip() +
  theme(legend.position = "none")

Se puede apreciar en la gráfica, que la mayor cantidad de accidentes ocurren en La Candelaria, incluso la cantidad de accidentes casi que duplica los accidentes en Laureles, que es la comuna que le sigue en cantidad de accidentes. Por otro lado, la comuna con menor cantidad de accidentes es Palmitas.

3.8. Número de accidentes por su gravedad

ggplot(data = gravedad, aes(fill = GRAVEDAD, x = GRAVEDAD, y = numero_de_accidentes)) +
  geom_bar(stat = "identity") + 
  scale_fill_manual(values = c("green", "red", "blue")) +
  geom_text(aes(y = numero_de_accidentes, label = numero_de_accidentes), vjust = -0.5) +
  xlab("Gravedad") + 
  ylab("Número de acccidentes") + 
  ggtitle("Total Accidentes por gravedad") +
  ylim(c(0,150000))+
  theme(legend.position = "none")

De los accidentes registrados de 2014 a 2020, visualizando la gravedad en esta gráfica podemos apreciar que la cantidad de muertos respecto a la cantidad de accidentes total es poca, y con heridos en la mayoria de accidentes aunque hablemos de muertes, es un alivio que en la mayor parte de los accidentes sólo sean daños o personas heridas.

---

4. Entrenamiento de un modelo predictivo

En la etapa de entrenamiento del modelo predictivo, se utilizó el registro histórico de accidentes desde el año 2014 hasta el año 2017, y para la etapa de validación se hizo uso de los registros de los años 2018 y 2019. Ésto debido a las especificaciones del trabajo para predecir la accidentalidad en Medellín.

4.1. Modelo Lineal

Inicialmente se utiliza un modelo lineal con las variables Festividad, Día Semana y Diseño.

En este modelo se va a observar el MSE (Error Cuadrático Medio) y el \(R^2\) (Coeficiente de Determinación) para determinar la potencia del modelo para predecir.

4.1.1. Predicción y Evaluación para los datos de Entrenamiento

lm1_data <- base_final03 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                   DISENO) %>% dplyr::count(name = "NRO_ACCID") 
lm1_tr <- lm1_data[,-c(5)]

predicted <- round(predict(lm1, newdata=lm1_tr))
actual <- lm1_data$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: 116.151546, MAE: 5.882622, R2: 0.896806"

4.1.2. Predicción y Evaluación para los datos de Validación en el año 2018

lm1_2018 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                  DISENO) %>% dplyr::count(name = "NRO_ACCID")

predicted <- round(predict(lm1, newdata=lm1_2018))
actual <- lm1_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"

La diferencia del MSE entre los datos de entrenamiento y validación del 2018 es del 34.48%, que al ser mayor que el 15% indica un posible sobreajuste. Además se puede apreciar que el \(R^2\) de los datos de validación para el año 2018 predice un 72.86%, sin embargo disminuyó un 16.82% en cuanto al \(R^2\) para los datos de entrenamiento. Así que luego de ésto se decide igualmente ver qué sucede con el mismo modelo validando con el año 2019.

4.1.3. Predicción y Evaluación para los datos de Validación en el año 2019

datos_vl02 <- subset(base_final, (AÑO == '2019'))

lm1_2019 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, DIA_SEMANA, 
                                  DISENO) %>% dplyr::count(name = "NRO_ACCID")

predicted <- round(predict(lm1, newdata=lm1_2019))
actual <- lm1_2019$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: 157.439013, MAE: 7.161181, R2: 0.765329"

La diferencia del MSE entre los datos de entrenamiento y validación del 2019 es del 41.29%, que al ser mayor que el 15% indica un posible sobreajuste. Además se puede evidenciar que el \(R^2\) de los datos de validación para el año 2019 predice un 76.53%, sin embargo aunque mejoró con respecto a la validación del año 2018, disminuyó un 13.15% en cuanto al \(R^2\) para los datos de entrenamiento. Por tanto, al obtener estos resultados validando con los años 2018 y 2019 se decide buscar un nuevo modelo cambiando las variables.

4.2. Modelo Lineal con disminución de Variables

Para este nuevo modelo se decide utilizar un modelo lineal únicamente con las variables “Festividad” y “Día Semana”. Es decir, se omite en este caso “Diseño” para observar qué cambios pueden ocurrir en el modelo.

4.2.1. 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"

4.2.2. Predicción y Evaluación para los datos de Validación en el 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"

Con este nuevo modelo disminuyendo variables no se obtuvieron resultados positivos, ya que el R2 disminuyó notablemente tanto en el entrenamiento, como en la validación con el año 2018 y la diferencia entre el MSE de datos de entrenamiento y validación para el año 2018 fue de 57.89%, lo cual indica que aumentó, indicando así una alta variabilidad en cuanto a las predicciones y de esa forma, una baja variabilidad explicada por este nuevo modelo.

4.2.3. Predicción y Evaluación para los datos de Validación en el 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 este modelo la diferencia del MSE entre los datos de entrenamiento y validación del 2019 es del 91.02%, que claramente indica un sobreajuste. Y se puede observar que el \(R^2\) de los datos de validación para el año 2019 predice un 10.07%. Así que se concluye que este modelo no sirve para predecir según los resultados obtenidos tanto en el entrenamiento, como para la validación en los años de 2018 y 2019. Así que se decide utilizar un modelo lineal generalizado.

4.3. Modelo Lineal Generalizado

4.3.1. 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"

4.3.2. Predicción y Evaluación para los datos de Validación en el 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 este modelo lineal generalizado con la familia de distribución Poisson se obtuvo un \(R^2\) en la etapa de entrenamiento de 7.21%, y para la etapa de validación para el año 2018 el \(R^2\) fue de 2.88%,lo cual indica que dicho modelo no sirve para predecir según los resultados obtenidos. Además la diferencia entre el MSE de entrenamiento y validación para el año 2018 fue de 57.89%, que indica un sobreajuste. Sin embargo se procede a validar igualmente con el año 2019.

4.3.3. Predicción y Evaluación para los datos de Validación en el 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"

Para este modelo la diferencia del MSE entre los datos de entrenamiento y validación del 2019 es del 90.91%, que claramente indica un sobreajuste. Y se puede observar que el \(R^2\) de los datos de validación para el año 2019 predice un 10.1%. Así que se evidencia que este modelo no sirve para predecir según los resultados obtenidos tanto en el entrenamiento, como para la validación en los años de 2018 y 2019. Así que se decide utilizar el modelo lineal generalizado adicinando otra variable.

4.4. Modelo lineal generalizado con Adición de la variable Clase

4.4.1. Predicción y Evaluación para los datos de Entrenamiento

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.678030, MAE: 5.765783, R2 Score: 0.886021"

4.4.2. Predicción y Evaluación para los datos de Validación en el 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"

En este modelo lineal generalizado con la familia de distribución Poisson con adición de la variable Clase se obtuvo un \(R^2\) en la etapa de entrenamiento de 88.6%,y para la etapa de validación para el año 2018 el \(R^2\) fue de 89.51%,lo cual indica que dicho modelo sirve para predecir según los resultados obtenidos. Además la diferencia entre el MSE de entrenamiento y validación para el año 2018 fue de 10.03%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. Luego se procede a validar con el año 2019.

4.4.3. Predicción y Evaluación para los datos de Validación en el 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"

En este modelo la diferencia del MSE entre los datos de entrenamiento y validación del 2019 es del 0.94%, que fue menor al valor obtenido con la validación del 2018 (10.03%), 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.63%, 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. Sin embargo, se adicionará otra variable para ver si se obtiene un mejor modelo.

4.5. Modelo lineal generalizado con Adición de la variable Diseño

4.5.1. 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.509575, MAE: 3.239037, R2 Score: 0.865613"

4.5.2. Predicción y Evaluación para los datos de Validación en el 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"

En este modelo lineal generalizado con la familia de distribución Poisson con adición de la variable Diseño se obtuvo un \(R^2\) en la etapa de entrenamiento de 86.56%, y para la etapa de validación para el año 2018 el \(R^2\) fue de 81.93%,lo cual indica que dicho modelo sirve para predecir según los resultados obtenidos. La diferencia entre el MSE de entrenamiento y validación para el año 2018 fue de 4.92%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. Así, se procede a validar con el año 2019.

4.5.3. Predicción y Evaluación para los datos de Validación en el 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"

En este modelo la diferencia del MSE entre los datos de entrenamiento y validación del 2019 es del 3.4%, que indica que no hay problemas de sobreentrenamiento. También se puede observar que el \(R^2\) de los datos de validación para el año 2019 predice un 82.91%, que evidencia que este modelo lineal generalizado con la adición de la variable Diseño es buen candidato para predecir la accidentalidad en Medellín.

Este modelo presenta el mejor MSE, pero no es un modelo viable ya que no es posible obtener la variable DISENO para realizar predicciones, según las pautas dadas.

Finalmente, se decide trabajar con el modelo lineal generalizado que posee la variable CLASE, ya que como se analizó anteriormente, se observa que tanto para los datos de entrenamiento, como para los de validación (2018, 2019) se pueden obtener buenas predicciones.

Luego, con el modelo elegido se procedió a realizar la etapa de entrenamiento y validación de forma semanal y mensual ya que dicho modelo para la predicción diaria tuvo buenos resultados, tanto para el MSE, como para el \(R^2\), para así observar si también con dicho modelo se pueden obtener buenas predicciones no solamente de manera diaria.

4.6. Semanal

4.6.1. Modelo lineal generalizado con adición de la variable Clase

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

datos_lm6_p <- datos_lm6[,-5]
y_train <- round(predict(lm6, newdata= datos_lm6_p, type="response"))
y_actual <- datos_lm6$NRO_ACCID
lm6_tmse <- MSE(y_train, y_actual)
lm6_tmae <-  MAE(y_train, y_actual)
lm6_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", 
        lm6_tmse, lm6_tmae, lm6_r2)

## [1] "MSE: 96.627665, MAE: 5.730144, R2 Score: 0.888415"

4.6.1.2. Predicción y Evaluación para los datos de Validación en el año 2018

datos_lm6_v1 <- datos_vl %>% group_by(FECHA, FESTIVIDAD, SEMANA, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm6_v2 <- datos_lm6_v1[,-5]

y_train <- round(predict(lm6, newdata= datos_lm6_v2, type="response"))
y_actual <- datos_lm6_v1$NRO_ACCID
lm6_tmse <- MSE(y_train, y_actual)
lm6_tmae <-  MAE(y_train, y_actual)
lm6_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm6_tmse, lm6_tmae, lm6_r2)

## [1] "MSE: 87.475847, MAE: 5.460300, R2 Score: 0.896469"

En este modelo lineal generalizado con la familia de distribución Poisson con adición de la variable Clase se obtuvo un \(R^2\) en la etapa de entrenamiento de 88.84%,y para la etapa de validación para el año 2018 el \(R^2\) fue de 89.65%,lo cual indica que dicho modelo sirve para predecir según los resultados obtenidos. Además la diferencia entre el MSE de entrenamiento y validación para el año 2018 fue de 9.15%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. Luego se procede a validar con el año 2019.

4.6.1.3. Predicción y Evaluación para los datos de Validación en el año 2019

datos_lm6_v1 <- datos_vl02 %>% group_by(FECHA, FESTIVIDAD, SEMANA, 
                                      CLASE) %>% dplyr::count(name = "NRO_ACCID")
datos_lm6_v2 <- datos_lm6_v1[,-5]

y_train <- round(predict(lm6, newdata= datos_lm6_v2, type="response"))
y_actual <- datos_lm6_v1$NRO_ACCID
lm6_tmse <- MSE(y_train, y_actual)
lm6_tmae <-  MAE(y_train, y_actual)
lm6_r2 <- R2_Score(y_train, y_actual)
sprintf("MSE: %f, MAE: %f, R2 Score: %f", lm6_tmse, lm6_tmae, lm6_r2)

## [1] "MSE: 97.902198, MAE: 5.926374, R2 Score: 0.888316"

En este modelo la diferencia del MSE entre los datos de entrenamiento y validación del 2019 es del 1.27%, que fue menor al valor obtenido con la validación del 2018 (9.15%), 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.

4.7. Mensual

4.7.1. Modelo lineal generalizado con adición de la variable Clase

4.7.1.1. 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.602242, MAE: 5.728723, R2 Score: 0.888444"

4.7.1.2. Predicción y evaluación para los datos de validación en el 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"

En este modelo lineal generalizado con la familia de distribución Poisson con adición de la variable Clase se obtuvo un \(R^2\) en la etapa de entrenamiento de 88.84%,y para la etapa de validación para el año 2018 el \(R^2\) fue de 88.81%,lo cual indica que dicho modelo sirve para predecir según los resultados obtenidos. Además la diferencia entre el MSE de entrenamiento y validación para el año 2018 fue de 2.19%, lo que indica que al ser menor del 15% no hay problemas de sobreentrenamiento. Luego se procede a validar con el año 2019.

4.7.1.3. Predicción y evaluación para los datos de validación en el 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.

5. Modelamiento - Predicción

Ahora, se muestra gráficamente la potencia de la predicción para el año 2020, según los datos reales de 2020 que se poseen desde la fecha 2020-01-01, hasta la fecha 2020-08-31.

Para ello se creó una base que contuviese los datos del año 2020 real y el del año 2020 con las predicciones; como en este año se tomó desde 2020-01-01 hasta 2020-12-31 porque es el predictivo, entonces para esta comparación se decide que va hasta 2020-08-31, para que esta posea la misma longitud de los datos que se tienen reales.

Igualmente, se diseñó una función para contener la relación entre la desviación estándar y la raíz cuadrada de la longitud. Para luego con ayuda de esto y junto a las especificaciones del gráfico, observar la cantidad de accidentes reales y predichos mediante boxplots, desde el mes 01 hasta el mes 08.

#Se cre? una funci?n que contuviera la relaci?n entre desviaci?n est?ndar y la ra?z cuadrada de la longitud
sem <- function(x, na.rm = FALSE) {
  out <- sd(x, na.rm = na.rm)/sqrt(length(x))
  return(out)}

#Se a?ade la base de datos especificada
datos_real_pred <- read_excel("Base_2020_real_predict.xlsx")

#Se hacen los agrupamientos necesarios y se omite la varible Incendio para que en ambos a?os est?n las mismas 'Clases'

datos_real_pred_001 <- datos_real_pred[datos_real_pred$CLASE!="Incendio", ]

datos_real_pred_001$AÑO<-as.integer(datos_real_pred_001$AÑO)

datos_real_pred_01 <- datos_real_pred_001 %>% group_by(AÑO, MES, CLASE, NRO_ACCID) %>% 
  dplyr::count() %>% mutate(NRO_ACCID=as.integer(NRO_ACCID))

#Gr?fico Comparativo de a?os reales y predictivos del 2020
datos_real_pred_01 %>% 
  group_by(MES, AÑO) %>% 
  dplyr::summarize(Número_de_Accidentes = mean(NRO_ACCID),
            se = sem(NRO_ACCID)) %>% 
  ggplot(aes(x = MES, 
             y = Número_de_Accidentes, 
             group = AÑO,  
             color = AÑO)) +
  geom_point() +
  geom_errorbar(aes(ymin = Número_de_Accidentes - se, 
                    ymax = Número_de_Accidentes + se)) +
  geom_line()

Nota: se aclara que en esta gráfica los años que se comparan son dos (2020 con datos reales y 2020 con los datos predictivos). En ese orden, se hace referencia a ‘año 1’ como el año 2020 con los datos reales que se poseen, y de la misma forma, se hace referencia a ‘año 2’, como el año 2020 con las predicciones. Por tanto, se aconseja hacer caso omiso a los valores de años 1.25, 1.5 y 1.75.

En la gráfica se observa que el año 1 (representado con la línea y boxplots color negro) tiene menos variabilidad ya que son los valores reales que se poseen de la cantidad de accidentes durante el tiempo que se especifica anteriormente, e igualmente se puede ver su promedio de números de accidentes a través de los 8 meses, representado por el punto de los diagramas de Cajas y Bigotes.

Luego, para el año 2 (representado con la línea y boxplots color azul) se nota más variabilidad, lo cual tiene sentido ya que hace referencia a los valores predictivos según el modelo empleado, e igualmente se puede ver su promedio de números de accidentes a través de los 8 meses, representado por el punto que se observa para cada Boxplot.

Se nota claramente que para los dos primeros meses del año 2020 el promedio de número de accidentes del ‘año 2’, que contiene los valores predictivos, son muy similares a los que se observan del año 1, con los valores reales. Luego, para los demás meses del año 2020 se puede apreciar que el año con los datos reales (año 1) su promedio de accidentes bajó, observándose así que el año con las predicciones quedó por encima de éste. Sin embargo, ésto no significa que el año 2 tenga malas predicciones, ni tampoco quiere decir que el modelo no tenía potencia predictiva, ya que anteriormente se había notado su capacidad predictiva en el estudio de la Accidentalidad de Medellín según su MSE y \(R^2\), lo que realmente sucede es que se debe poner en contexto lo que sucedía en el país en dichos meses.

Para el año 2020 en Colombia, El Ministerio de Salud y Protección Social confirma el primer caso de COVID-19 (enfermedad infecciosa, causada por el coronavirus) en el mes 3, a causa de ello el 8 de marzo de 2020, el presidente Iván Duque Márquez dio a conocer las acciones de contención en el país. Por consiguiente, los ciudadanos adoptando dichas medidas y con el propósito de que no aumentara la tasa de contagios por coronavirus, empezaron a tener conciencia social y a cuidarse y así mismo a cuidar a los demás. De hecho, desde febrero que ya se estaba especulando que el COVID-19 iba a llegar al país, muchas personas ya se estaban cuidando.

Así, entre las medidas que se dieron, se encuentra la de cierre de vías en marzo, que es por eso que gráficamente se notó un declive para los datos del promedio de accidentes a partir del mes 3 para los datos reales, ‘año 1’. Además de que luego se implementó cuarentena obligatoria por el causal de pandemia, que fue conocido como ‘El Aislamiento Preventivo Obligatorio en Colombia o confinamiento de Colombia de 2020’ que se inició el 25 de marzo de 2020, y que en un principio se había decretado por 19 días, pero luego fue variando y extendiéndose dicha cuarentena, según la tasa de contagios.

Para el año 2021 no se poseen valores reales, puesto que no tenemos datos con respecto a dicha información.

A continuación, se presentan las predicciones diarias, semanales y mensuales que se obtuvieron para los años 2020 y 2021.

5.1. Predicción diaria del 2020

Se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones diarias obtenidas 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 != '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_diaria2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020,0))

diario_20_02 <- prediccion_diaria2020 %>%
  group_by(FECHA, DIA_SEMANA, CLASE, FESTIVIDAD) %>%
  dplyr::summarise(NRO_TOTAL_ACCID=NRO_ACCID)

head(diario_20_02, 10)

5.2. Predicción diaria del 2021

A continuación, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones diarias obtenidas para el año 2021:

Base_prediccion_2021 <- subset(Base_prediccion, (AÑO != '2020'))

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) %>%
  dplyr::summarise(NRO_TOTAL_ACCID=NRO_ACCID)

head(diario_21_02, 10)

5.3. Predicción semanal del 2020

De igual forma, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones semanales que se obtuvieron para el año 2020:

Base_prediccion03 <- Base_prediccion[,-4]
Base_prediccion04 <- Base_prediccion03[,-4]

Base_prediccion_2020 <- subset(Base_prediccion, (AÑO != '2021'))

Base_prediccion_2020 <- subset(Base_prediccion04, (AÑO != '2021'))


Base_prediccion_2020$SEMANA <- as.integer(Base_prediccion_2020$SEMANA)


prediccion_2020_02 <- predict(object = lm6, newdata = Base_prediccion_2020,
                          type = "response")
prediccion_semanal2020 <- Base_prediccion_2020 %>% 
  mutate(NRO_ACCID = round(prediccion_2020_02,0))

semanal <- prediccion_semanal2020 %>% group_by(CLASE, SEMANA = week(FECHA), NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
semanal <- mutate(semanal, NRO_ACCID_TOTAL=NRO_ACCID*total)

semanal_20_02 <- semanal %>%
  group_by(SEMANA, CLASE, FESTIVIDAD) %>%
  dplyr::summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))

head(semanal_20_02, 10)

5.4. Predicción semanal del 2021

También, acá se puede observar una tabla con el encabezado de las primeras 10 observaciones para las predicciones semanales que se obtuvieron para el año 2021:

base_final03 <- subset(base_final02, (AÑO != '2020'))
Base_prediccion_2021 <- subset(Base_prediccion, (AÑO != '2020'))

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_semanal2021 <- Base_prediccion_2021 %>% 
  mutate(NRO_ACCID = round(prediccion_2021,0))

#borrando columnas no necesarias
prediccion_semanal2021 <-  prediccion_semanal2021[,c(-1,-2,-4,-5)]

#Se agrupó por semana 2021
semanal_21 <- prediccion_semanal2021 %>% group_by(CLASE, SEMANA, NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
semanal_21 <- mutate(semanal_21, NRO_ACCID_TOTAL=NRO_ACCID*total)

semanal_21_02 <- semanal_21 %>%
  group_by(SEMANA, CLASE, FESTIVIDAD) %>%
  dplyr::summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))

head(semanal_21_02, 10)

5.5. Predicción mensual del 2020

De manera resumida, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones mensuales que se obtuvieron para el año 2020:

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) %>%
  dplyr::summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))

head(mensual_20_02, 10)

5.6. Predicción mensual del 2021

Igualmente, se presenta una tabla con el encabezado de las primeras 10 observaciones para las predicciones mensuales que se obtuvieron para el año 2021:

Base_prediccion_2021 <- subset(Base_prediccion, (AÑO != '2020'))

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_mensual2021 <- Base_prediccion_2021 %>% 
  mutate(NRO_ACCID = round(prediccion_2021,0))

#Se borranron columnas no necesarias
prediccion_mensual2021 <-  prediccion_mensual2021[,c(-1,-2,-3,-5,-7)]

#Agrupamiento por mes 2021

#Se agrupó por mes y si fue en día festivo o no
mensual_21 <- prediccion_mensual2021 %>% group_by(CLASE, MES, NRO_ACCID, FESTIVIDAD) %>% dplyr::summarize(total = n())
mensual_21 <- mutate(mensual_21, NRO_ACCID_TOTAL=NRO_ACCID*total)

mensual_21_02 <- mensual_21 %>%
  group_by(MES, CLASE, FESTIVIDAD) %>%
  dplyr::summarise(NRO_TOTAL_ACCID=sum(NRO_ACCID_TOTAL))

head(mensual_21_02, 10)

---

6. Agrupamiento de los barrios de Medellín de acuerdo a su accidentalidad

Antes de realizar la agrupación de barrios según la accidentalidad se decide mostrar primero un mapa de calor de accidentalidad en la ciudad de Medellín entre los años de 2014 y 2019, ya que este mapa a su vez representaría la historia de la accidentalidad. Lo anterior, teniendo presente que para este proyecto se pretende predecir para los años de 2020 y 2021.

6.1. Mapa de accidentes en la ciudad de Medellín

Para la elaboración del mapa calor según la accidentalidad, se descargó un archivo .shp del Límite Barrio Vereda Catastral

#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()) %>%
  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)

# Cantidad de Accidentes por Comuna

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", position = "dodge", fill = "blue", 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 entre los años 2014 y 2019") +
  ylim(c(0,50000)) +
  theme_minimal() +
  coord_flip()

Al analizar la accidentalidad por Barrio y Comuna, tanto en el mapa de accidentalidad como en el gráfico, se observó que La Candelaria es el lugar donde más se presentan accidentes entre los años 2014 y 2019, seguida por Laureles y Castilla. También, se evidencia que Palmitas es el lugar donde ocurre menos accidentalidad en Medellín durante los años 2014-2019. Además en el mapa de calor se puede notar que la zona centro y las vías principales son las más afectadas por los accidentes.

Luego se procede a crear una función para calular distancias para datos geoespaciales y además se crea un dendograma para así emprender la búsqueda del \(K\) óptimo para el agrupamiento (tomando los datos de 2014 hasta 2017, ya que éstos fueron los años que se utilizaron para el modelo de predicción).

#Se realizó la conversión de la latitud y longitud al formato numérico

base_final03$LATITUD <- as.numeric(as.character(base_final03$LATITUD))
base_final03$LONGITUD<- as.numeric(as.character(base_final03$LONGITUD))

#Se creó con la función 'geo.dist', una matriz de distancias
df1 <- df[1:1000, ]
d <- geo.dist(df1)
hc <- hclust(d)
plot(hc, main = "Dendograma", col = "#00AFBB")

df1$clust <- cutree(hc, k = 6)
head(df1,10)

6.2. Mapa de agrupamiento según latitud y longitud.

Para la realización del mapa de agrupamiento según latitud y longitud, se descargó un archivo .shp del Límite Catastral de Comunas y Corregimientos.

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")+
  coord_fixed()

El anterior mapa muestra una posible agrupación, según las medidas geoespaciales de la latitud y la longitud de los accidentes. Sin embargo, este agrupamiento se utiliza como referencia porque para su creación no se utilizó ningún método para la elección del \(K\) óptimo.

Así que se procede a realizar una clusterización del número de accidentes por Gravedad y Barrio haciendo uso del algoritmo “k means” para la búsqueda del \(K\) óptimo.

6.3. Clusterización con número de accidentes por gravedad y barrio.

Con los datos preprocesados y el subconjunto de datos que se seleccionaron, Se les realizó un escalamiento y centrado de la base de datos.

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

Para este caso, con un k=5, se evidencia que el valor de Suma de Cuadrados entre grupos (SS between) sobre la Suma de Cuadrados Totales fue de aproximadamente 85.1% (0.851), que indica un buen ajuste porque es cercano a 1. Sin embargo, es mejor graficar el WSS contra el número de clúster, ya que este número se debe especificar de antemano.

Luego se procede a hallar el k óptimo.

6.4. El método del codo

#Se fijó una semilla y se realizó el cálculo y se gráficó el WSS(total within - cluster sum of square) para k = 2 hasta k = 10
set.seed(2021022)
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")

Según la gráfica del Método del Codo posiblemente el k=4 o k=5 serían buenos candidatos para el k óptimo, ya que presentan un cambio más suave en las pendientes en comparación con k=2 o k=3. Igualmente al observar el between SS / total SS para k=5, mencionado anteriomente, se evidencia un 85.1 %, lo cual indica un buen ajuste. Además como se graficó el WSS contra el número de clústeres, se refleja que es un buen candidato.

Después, se busca el \(K\) óptimo haciendo uso del paquete NbClust, de la siguiente forma:

nb <- NbClust(scaled_data, diss=NULL, distance = "euclidean", 
              min.nc=4, max.nc=8, method = "kmeans", 
              index = "all", alphaBeale = 0.1)

## *** : The Hubert index is a graphical method of determining the number of clusters.
##                 In the plot of Hubert index, we seek a significant knee that corresponds to a 
##                 significant increase of the value of the measure i.e the significant peak in Hubert
##                 index second differences plot. 
## 

## *** : The D index is a graphical method of determining the number of clusters. 
##                 In the plot of D index, we seek a significant knee (the significant peak in Dindex
##                 second differences plot) that corresponds to a significant increase of the value of
##                 the measure. 
##  
## ******************************************************************* 
## * Among all indices:                                                
## * 13 proposed 4 as the best number of clusters 
## * 2 proposed 5 as the best number of clusters 
## * 7 proposed 6 as the best number of clusters 
## * 2 proposed 8 as the best number of clusters 
## 
##                    ***** Conclusion *****                            
##  
## * According to the majority rule, the best number of clusters is  4 
##  
##  
## *******************************************************************

El cual sugiere que el \(K\) óptimo es 4.

hist(nb$Best.nc[1,], breaks = max(na.omit(nb$Best.nc[1,])), main = "Histograma del K óptimo ", xlab = "K", ylab = "Frecuencia", col="darkorchid4")

Según el histograma, también indica que el \(K\) óptimo es el 4.

6.5. Resúmen de métodos

En este resumen se encuentra el método de la Silueta, el Método del Codo y Gap Statistic. Donde se obtuvieron los siguientes resultados:

6.5.1. Método de la silueta.

fviz_nbclust(scaled_data, kmeans, method = c("silhouette"))

6.5.2. Método del codo.

fviz_nbclust(scaled_data, kmeans, method = "wss") +
    geom_vline(xintercept = 4, linetype = 2)+
  labs(subtitle = "Método del Codo")

6.5.3. Gap statistic

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

6.6. Generación de clusterización, según el k óptimo.

Según los diferentes métodos \(k = 4\), parecía ser muy óptimo para la generación de la clusterización.

Luego, se muestran las primeras 10 observaciones de los barrios ordenados alfabéticamente, donde se observa el tipo de gravedad y el grupo al cual pertenecen.

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)

Con lo realizado anteriormente poseemos el agrupamiento según el k óptimo igual a 4, es decir, ahora los barrios pertenecen a un tipo de cluster (enumaredado del 1 al 4). Para el nombramiento de cada grupo se tomó como referencia el mapa de calor, quedando de la siguiente forma:

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

La información que se obtuvo de cada grupo es la siguiente:

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

Para este grupo se observa que la cantidad total de accidentes son de 27895, la cantidad de accidentes con heridos es de 13145, la cantidad de accidentes con muertos 111 y la cantidad de accidentes con solo daños es de 14639.

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

Para este grupo se observa que la cantidad total de accidentes son de 28622, la cantidad de accidentes con heridos es de 17817, la cantidad de accidentes con muertos 106 y la cantidad de accidentes con solo daños es de 10699.

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

Para este grupo se observa que la cantidad total de accidentes son de 36247, la cantidad de accidentes con heridos es de 17147, la cantidad de accidentes con muertos 237 y la cantidad de accidentes con solo daños es de 18863.

6.6.4. Grupo 4 - accidentalidad media-alta

Para este grupo se observa que la cantidad total de accidentes son de 54841, la cantidad de accidentes con heridos es de 32887, la cantidad de accidentes con muertos 222 y la cantidad de accidentes con solo daños es de 21732.

Finalmente, se realizó el mapa de agrupamiento, según lo mostrado anteriormente.

6.7. Mapa de accidentalidad en la ciudad de Medellín por agrupamiento

#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("#00FF66", "#CCFF00", "#FF0000", "#0066FF")
mapa02$colores <- ifelse(mapa02$kmm.cluster == "1", "#00FF66",
                            ifelse(mapa02$kmm.cluster == "2", "#CCFF00",
                                   ifelse(mapa02$kmm.cluster == "3", "#FF0000",
                                          ifelse(mapa02$kmm.cluster == "4", "#0066FF",0))))

#Mapa final
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"))

---

7. Conclusiones

• El lunes, es el día de la semana que más accidentes se presentan en la ciudad de Medellín.

• El modelo lineal generalizado con la familia de distribución Poisson y con la adición de la variable Clase, demostró eficiencia y es un buen candidato para predecir la accidentalidad en Medellín. Con los resultados obtenidos, no se observaron problemas de sobreentrenamiento para este modelo.

• La mayor cantidad de accidentes ocurren en la comuna La Candelaria. La zona centro y las vías principales son las más afectadas por los accidentes.

• Con el k óptimo obtenido se puede notar que la agrupación sí representa el nivel de accidentalidad por clusterización en la ciudad de Medellín. Lo cual indica que los cuatro grupos muestran de manera eficiente lo que ocurre en la ciudad, en cuanto a la accidentalidad.

---

8. Comentarios

• Con un poco más de tiempo, se podrían implementar modelos con redes neuronales o también con potenciación del gradiente.

• Sería interesante analizar los datos actuales con una nueva variable que nos indique el estado del clima del día en que ocurrió el accidente. Podríamos definir una variable que puede tener las etiquetas día soleado o día lluvioso, y de esta forma verificar el comportamiento de los accidentes dependiendo del clima.

• Para las predicciones del año 2020, se obtienen resultados por encima de los valores reales, esto es debido a las restricciones de movilidad en el país, por la emergencia sanitaria ocasionada por el COVID-19, más no porque el modelo utilizado no sea eficiente.

• Las herramientas virtuales que se utilizaron (Google Meet, el software estadístico de R, R-Studio y WhatsApp, principalmente), así como nuestro compromiso con el proyecto y nuestros conocimientos en Estadística, hizo posible que desarrolláramos este trabajo de forma eficaz, a pesar de no estar presentes de forma presencial. Igualmente, el acompañamiento del docente del curso y su disposición a escuchar las dudas en las clases programadas, hizo posible que tuviéramos éxito en esta investigación.

---

9. Enlaces

Es importante resaltar que este reporte trae consigo dos elementos adicionales complementarios:

• Un Repositorio en Github, donde usted puede acceder a las bases de datos utilizadas a lo largo del reporte, y también a todos los scripts utilizados para obtener los resultados.
• Una Aplicación web elaborada con R Shiny Web App, en donde usted puede seleccionar una ventana de tiempo y obtener los datos históricos de accidentalidad por tipo de accidente. También podrá predecir la accidentalidad por tipo de accidente, visualizar los grupos de barrios, entre otras.

---

10. Referencias

[1]Incidentes viales. (2021). URL: http://medata.gov.co/dataset/incidentes-viales

[2]Alcaldía de Medellín OpenData, URL:https://geomedellin-m-medellin.opendata.arcgis.com/datasets/M-Medellin::limite-barrio-vereda-catastral/about

[3]R Markdown Cookbook, Yihui Xie, Christophe Dervieux, Emily Riederer (2021). URL:https://bookdown.org/yihui/rmarkdown-cookbook/

[4]R Markdown Theme Gallery. Andrew Zieffler (2018). URL:https://www.datadreaming.org/post/r-markdown-theme-gallery/

[5] Frequently asked questions. Tal Galili, (2021). URL: https://cran.r-project.org/web/packages/dendextend/vignettes/FAQ.html

[6] ¿Qué significa ss total y entre ss en el agrupamiento k-means? (s.f.) URL:https://qastack.mx/stats/82776/what-does-total-ss-and-between-ss-mean-in-k-means-clustering

[7]Modelos Estadísticos Avanzados. Pedro A. Torres-Saavedra, (2018). URL:https://pegasus.uprm.edu/~pedro.torres/book/

[8]Basic Aesthetics. URL:https://bookdown.org/yih_huynh/Guide-to-R-Book/basic-aesthetics.html