TAE - Accidentalidad en Medellín - Reporte Técnico

1 Introducción

Los accidentes viales son una problemática que ha afectado de gran manera a la ciudad de Medellín. El crecimiento económico de la ciudad y la rápida urbanización causan un gran aumento en el parque automotor, la densidad vehícular y las dificultades de transitar en las vías que conforman la ciudad. En este trabajo se plantea crear un sistema que pueda ayudar a mitigar los accidentes en las vías, además de servir como guía para tomar medidas frente a la situación actual.

Durante el periodo entre 2014 y 2018 se logró tomar una cantidad significativa de datos sobre la accidentalidad en Medellín y algunas de sus zonas aledañas, para evidenciar así, el problema de movilidad que presenta el área. Estos datos se encuentran en el portal GeoMedellín, donde cada uno de estos incidentes se encuentra georreferenciado, dividido por barrio y comuna. Cuenta, además, con la información sobre la vía donde el accidente ocurrió y también sobre el tipo de accidente que se dio, la fecha de éste e incluso los daños perjudiciales causados a los ciudadanos -muertos, heridos o solo daños materiales-.

Ahora bien, al tener más de 200.000 datos durante este intervalo, el equipo se planteó el objetivo de darle un buen uso a estos, no sólo para analizar lo que sucedió en este lapso, sino también para ayudar a predecir accidentes a futuro y así poder esgrimir un modelo prevencionista. En ese orden de ideas, para acceder a la información recolectada y al presente análisis se habilitará un enlace a aplicación TransitApp, donde la información a la que se accede estará presentada de una forma interactiva e intuitiva, para obtener un mejor entendimiento del usuario.

En el presente reporte se encuentra información detallada respecto a la depuración, imputación, análisis descriptivo, creación de clusters y modelos predictivos realizados para crear la aplicación TransitApp. Se verá que los modelos mixtos (modelos jerárquicos) presentan una buena estimación de la accidentalidad (número de accidentes) en la ciudad de Medellín.

Enlace aplicativo TransitApp: https://valentina-garcia-velasquez.shinyapps.io/TransitAPP_2020/

Enlace repositorio GitHub: https://github.com/vagarciave/TransitApp_1.0/

2 Justificación

La ciudad de Medellín presenta un índice de accidentalidad vial bastante alto comparada con otras ciudades de la región e incluso con la capital del país, Bogotá. Según el “Ránking (sic) Latinoamericano de Ciudades Fatales” de 2018, informe realizado por el centro de pensamiento LA Network, Medellín estaba ranqueada en el puesto 29 de las ciudades con más muertes en las vías de Latinoamérica, con 10.9 muertos por cada 100 mil habitantes. Además, durante el año 2017 hubo un total 41.667 accidentes viales en Medellín, los cuales ocasionaron 258 muertes.

También es necesario tener en cuenta al sistema de salud en los accidentes viales, debido a que estos afectan gravemente la capacidad de cobertura para otro tipo de emergencias. Solo en 2017, 14000 de los 30000 heridos en los accidentes viales requirieron asistencia médica, y de este grupo, 3800 presentaron heridas graves. Aunque, durante los años 2019 y 2020 se ha notado una disminución de estas estadísticas sobre la accidentalidad vial, se puede evidenciar que existe la necesidad de tener un mayor entendimiento de los datos, especialmente para los conductores, motociclistas y hasta ciclistas y peatones. Con la implementación del sistema planteado en este trabajo, se pretende crear una política de prevención de accidentalidad y además una conciencia ciudadana para manejar con mayor precaución en ciertos lugares y poder planear mejor sus viajes.

3 Objetivos

El análisis y los modelos que se pretenden construir sólo abarca el área urbana de Medellín correspondiente a sus 16 comunas y sus respectivos barrios.

3.1 Objetivo general

Construir un sistema que le permita a la Secretaría de Movilidad de Medellín tomar decisiones en la creación de estrategias que ayuden a reducir la accidentalidad en zona específicas de Medellín, observando cuales son los tipos de accidentes más frecuentes en las diferentes comunas y sus barrios.

También se busca que el sistema les permita a todos los ciudadanos que se movilizan en un vehículo tomar precauciones respecto a los accidentes más frecuentes que suceden en lugares específicos.

Lo anterior creando una interfaz amigable para los usuarios y fácil de utilizar.

3.2 Objetivos específicos

• Construir modelos que predigan el número de accidentes de cada tipo, tomando como entradas una ventana y una resolución temporal específicas, además de una zona espacial que puede ser un barrio o comuna de Medellín.

• Elaborar un sistema que agrupe barrios con características similares por los tipos de accidentes (clase) y la gravedad que tienen lugar en estos.

• Crear una interfaz gráfica que le permita al usuario visualizar con mayor comodidad la accidentalidad que hay en Medellín y las predicciones de esta a futuro.

4 Descripción y Limpieza de los datos

Los datos utilizados se encuentran en las bases de datos de Geomedellín (Portal Geográfico del Municipio de Medellín) y contienen información referente a múltiples accidentes de tránsito, en los que se detalla el tipo de accidente, dónde y cuándo ocurrió. La base de datos utilizada en este trabajo es una unión de las bases de incidentes en 2014, 2015, 2016, 2017 y 2018, se usará este último año para validar los modelos planteados.

La base de datos original cuenta con las siguientes variables:

• X: Componente de coordenada.
• Y: Componente de coordenada.
• OBJECTID: Id de cada incidente.
• RADICADO: Código emitido por la secretaría de movilidad de Medellín.
• FECHA: Fecha del incidente.
• HORA: Hora del incidente.
• DIA: Día del mes en el que ocurre el incidente
• PERIODO: Año del incidente.
• CLASE: Tipo de accidente
• DIRECCION.
• DIRECCION_ENC.
• CBML: Código de ubicación del predio en la ciudad.
• TIPO_GEOCOD: Tipo de
• GRAVEDAD: repercusiones del accidente
• BARRIO.
• COMUNA.
• DISENO: clasificación del lugar del accidente.
• DIA_NOMBRE: Día de la semana en el que ocurre el incidente
• MES: Mes del incidente en número (del 1 al 12)
• MES_NOMBRE: Mes del incidente
• X_MAGNAMED: Componente de coordenada
• Y_MAGNAMED: Componente de coordenada
• LONGITUD: Componente de coordenada
• LATITUD: Componente de coordenada

Se deciden eliminar las variables “X”, “Y”, “X_MAGNAMED” y “Y_MAGNAMED” ya que cumplen la misma función que las variables “LONGITUD” y “LATITUD”. También se elimina la variable “RADICADO” ya que sirve para identificar a un respectivo incidente, al igual que la variable “OBJECTID”. Se elimina la variable “MES_NOMBRE” ya que es redundante el la base de datos, debido a que se encuentra la misma información en la variable “MES”. Además ponerle los nombres a partir de la variable “MES” traería un aumento al coste computacional innecesario. También se decide eliminar las variables “DIRECCION”, “DIRECCION_ENC”, “CBML”, “TIPO_GEOCOD” y “HORA” ya no serán de utilidad en el modelo.

Además se crea otra variable llamada “DIA_FESTIVO”. Para poder analizar el número de accidentes en los días festivos y demás fechas especiales. Por lo tanto, las variables de interés que se usan en este proyecto son las siguientes:

• OBJECTID: Id de cada incidente.
• CLASE: Tipo de accidente
• GRAVEDAD: repercusiones del accidente
• COMUNA.
• BARRIO.
• DISENO: clasificación del lugar del accidente.
• LATITUD: Componente de coordenada
• LONGITUD: Componente de coordenada
• FECHA: Fecha del incidente.
• DIA: Día del mes en el que ocurre el incidente
• MES: Mes del incidente en número (del 1 al 12)
• PERIODO: Año del incidente.
• DIA_NOMBRE: Día de la semana en el que ocurre el incidente
• DIA_FESTIVO.

Definido lo anterior, se realiza una limpieza en la base de datos, ya que hay registros u observaciones que están mal escritos o simplemente no deberían estar. A continuación se explican los cambios realizados.

4.1 Depuración

Se decide quitar los acentos a cada una de las variables categóricas debido a que hay valores repetidos los cuales representan el mismo nivel o categoria, pero con acento y sin acento. Un ejemplo de esto es el barrio Berlpín de la comuna Aranjuez, el cual aparece 647 veces con el nombre de “Berlín” y 15 veces con el nombre de “Berlín”. Por lo tanto se decide arreglar este problema.

4.1.1 Para COMUNA

La variable COMUNA debe tener las comunas urbanas de Medellín y sus corregimientos.

Se sabe que en Medellín hay 16 comunas urbanas y 5 corregimientos, ambos están compuestos por barrios, por lo que la variable debería tener 16 + 5 = 21 niveles. Pero al observar los niveles de la variable “COMUNA” en el conjunto de entrenamiento, se evidencian 84 comunas, por lo cual se decide buscar las razones de esto y tratar de corregirlo.

Comunas: Popular, Santa Cruz, Manrique, Aranjuez, Castilla, Doce de Octubre, Robledo, Villa Hermosa, Buenos Aires, La Candelaria, Laureles-Estadio, La América, San Javier, El Poblado, Guayabal, Belén.

Corregimientos: Corregimiento de San Cristóbal, Corregimiento de San Antonio de Prado, Corregimiento de Santa Elena, Corregimiento de Altavista, Corregimiento de San Sebastián de Palmitas.

Uno de los problemas que se encuentra es que la persona encargada de digitar los datos se confundió entre las variables “COMUNA” y “BARRIO”. Un ejemplo de esto es el barrio Boston que aparece como comuna en 2 observaciones, con su respectiva comuna (La candelaria) en barrio. Así como el ejemplo anterior, hay varios casos. Por lo cual, estas comunas y barrios se ponen en su posición correcta.

Luego de resolver el problema anterior, se encuentra que aún hay comuna que no pertenecen a las 21 mencionadas anteriormente, estas son: “In”, “AU”, “SN”, “0”. Al no encontrar información de estos valores en la página oficial o en internet, se decide reemplazar estos valores con sus respectivos barrios como datos faltantes o NA.

4.1.2 Para BARRIO

Se detecta que los barrios contienen algunos problemas de espacios, algunos ejemplos son:

• “Asomadera No. 1” y “Asomadera No.1”
• “Aures No. 2” y “Aures No.2”
• “Bombona No. 1” y “Bombona No.1”
• “B. Cerro El Volador” y “B. Cerro El Volador”

Por lo tanto se arregla este problema quitando el espacio que hay entre “No.” y su número correspondiente en los barrios que tienen esta característica. También se eliminan los espacios que hay antes de la primera letra y después de la última letra.

También se reemplaza “Barrios de Jesús” por “Barrio de Jesús”, “Nueva Villa de La Iguana” por “Nueva Villa de la Iguana”, “Santa María de Los Ángeles” por “Santa María de los Ángeles”, “Villa Lilliam” por “Villa Liliam”.

4.1.3 Para CLASE

“Caida de Ocupante” se reemplaza por “Caida Ocupante” ya que se consideran equivalentes. También se elimina la observación con CLASE = “Choque y Atropello” ya que solo hay una.

4.1.4 Para DIA_NOMBRE

Se eliminan los espacios que hay antes de la primera letra y después de la última letra en todas las observaciones.

4.2 Imputación

En los datos faltantes de la variable DISENO, se decide usar el nombre de “No Registrado” para reemplazarlos.

En la base de datos se tiene la siguiente cantidad de datos faltantes por variable: la variable BARRIO tiene 19957 de datos faltantes, COMUNA tiene 19957 y CLASE tiene 7. Las demás variables no contienen datos faltantes.

Como la variable CLASE solo tiene 7 datos faltantes, se decide eliminar estas observaciones. Por el contrario, como BARRIO y COMUNA presentan una cantidad importante de los datos, se decide realizar una imputación de los mismos.

4.2.1 Para COMUNA

Primero se decide realizar la imputación en la variable COMUNA, se utiliza un modelo knn (k-Nearest Neighbour) con variable respuesta COMUNA y como predictoras las variables LATITUD y LONGITUD escaladas. Dado que este es un modelo que usa la distancia euclidiana para clasificar observaciones, se piensa que funcionaría bien para clasificar a las comunas teniendo su ubicación en latitud y longitud.

Para ajustar el modelo se usó la función train del paquete caret, usando K-Fold CV con K=10 y el número de vecinos igual a 1, 4, 7, 10, 13, 16 y 19. Se entrenó el modelo usando los años de 2014 a 2017 y se validó con el año 2018. Se usó la precisión de prueba del modelo como función de optimización. Los resultados obtenidos son los siguientes:

Accuracy	Kappa	AccuracyLower	AccuracyUpper	AccuracyNull	AccuracyPValue
0.9985	0.9984	0.9981	0.9989	0.2062	0

Los resultados se consideran buenos y se decide realizar el reemplazo de los valores faltantes de COMUNA por medio de este modelo.

4.2.2 Para BARRIO

Se decide usar un modelo knn (k-Nearest Neighbour) con variable respuesta BARRIO y como predictoras las variables LATITUD y LONGITUD escaladas y COMUNA.

Dado que habían valores de barrios que están en el año 2018 pero no en los demás, no se pudo hacer uso del paquete caret ya que presentaba un error al entrenar el modelo, por lo que se decide usar la función knn.cv del paquete class, para medir el comportamiento del modelo knn con k = 1 mediante LOOCV, entrenando el modelo con los años de 2014 hasta 2017 y validando el resultado con el año 2018. Se obtiene una precisión de prueba de 0.9975, por lo que se decide usar este modelo para imputar los datos faltantes en BARRIO.

4.3 Base final

Finalmente se decide usar solo las comunas del área urbana de Medellín y dejar a un lado los corregimientos, es decir, solo se tendrán las 16 comunas mencionadas anteriormente en este reporte.

Tras las modificaciones mencionadas, de las 228693 observaciones originales entre las bases de datos del 2014 al 2018, se utilizarán 204002 observaciones en el análisis y construcción de los modelos. La base de datos luce de esta manera:

A continuación se encuentra el acceso al código utilizado en todo este proceso

Enlace: https://github.com/vagarciave/TransitApp_1.0/blob/master/limpieza/Depuracion.Rmd

5 Análisis descriptivo

5.1 Sobre la base de datos

Luego de la depuración e imputación de los datos, se cuenta con un total de 204002 incidentes de tránsito, los cuales tienen 14 características que se pueden llamar “variables”. De las 14 variables 7 son categorícas (OBJECTID, CLASE , GRAVEDAD, COMUNA, BARRIO, DISENO, DIA_NOMBRE). Hay un total de 16 comunas, 268 barrios, 6 clases de accidente, 3 niveles de gravedad, 13 tipos de diseño de vía.

Para este análisis no se considerarán las variables FECHA, LATITUD y LONGITUD.

5.2 Periodo

La cantidad total de accidentes por periodo es bastante simular, se nota una leve reducción para el 2018.

5.3 Clase

## List of 1
##  $ legend.position: chr "bottom"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

Para clase, es claro que en cada periodo lo que más se presentan son choques, de manera análoga es poco común que hayan incendios, de igual forma los volcamientos también son poco comunes. Para las clases atropello, caída de ocupante y otro parece haber cantidades muy similares, estos comportamientos son parecidos en cada periodo pero hay diferencias que podrían ser significativas en los modelos.

5.4 Gravedad

## List of 1
##  $ legend.position: chr "bottom"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

En los accidentes de tránsito lo más común es que se presenten heridos, seguido por solo daños, por otra parte el número de muertos para cada periodo es muy bajo, comparado con los otros tipos de gravedad. Para el año 2018 se nota una leve disminución en la cantidad de heridos.

5.5 Comuna

En los accidentes por comuna se evidencia que la comuna “La Candelaria” es la que mas accidentes ha registrado durante los periodos analizados, la de menor cantidad de accidentes es la comuna “Santa Cruz”, esto es con el número de accidentes observado, del 2014 al 2018.

5.6 Comunas por periodo

Ahora discriminando el número de accidentes en las comunas por periodo, se observa que se presenta una tendencia similar para cada periodo, por lo que se piensa que cada año van a existir tendencias similares de la accidentalidad por comuna.

5.7 Barrios por periodo

Para los barrios se tiene el cuenta el top 10 con más accidentes por cada periodo, el barrio que presentó mayor número de accidentes de 2014 a 2017 fue la Candelaria, para 2018 este barrio fue superado por una cantidad no muy grande por el barrio campo amor, se evidencia como a través del tiempo algunos barrios ingresan o salen del top 10, por ejemplo el Carlos E. Restrepo empezó a tener una accidentalidad más alta a partir del 2016, por su parte Guayaquil redujo su accidentalidad para 2018, el barrio Naranjal solo estuvo en el top 10 en el año 2015.

Los primeros 5 barrios con mas cantidad de accidentes por gravedad

Estos son los 5 barrios con mas cantidad de accidentes, los cuales se clasifican por la gravedad, que son:“solo daños”,“Herido” y “Muerto”, se identificaron los primeros barrios con mas accidentalidad por año, notar que la gravedad menos concurrente es “Muerto”, aunque la Candelaria es el barrio con mas accidentes por año, tamibien note que Castilla tiene mas accidentes con gravedad de muertes mas alta que la de los demas barrios.

5.8 Diseño

## List of 1
##  $ legend.position: chr "bottom"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

El diseño que más accidentalidad presenta es el tramo de vía, esto para todos los periodos. De 2014 a 2018 hubo una disminución en este tipo de accidentes, por su lado los accidentes para intersección y Lote o predio aumentaron levemente de 2014 a 2018.

5.9 Accidentalidad por mes

Se puede apreciar que la accidentalidad para 2018 es menor comparada con otros años, esto para casi todos los meses, se observa como hay diferencias para cada periodo y cada año, por lo cual una variable temporal en los modelos predictivos podría sr significativa. Para los meses de Enero a Marzo la accidentalidad parece ser más baja, para el mes de Junio los accidentes disminuyen. De Agosto a Octubre la accidentalidad en más alta.

5.10 Día de la semana

Es claro que los días festivos son un factor realmente influyente en el número de accidentes, pues para estos se presentan menor cantidad de accidentes, para el día domingo el total de accidentes es menor, para los demás días no se observa una diferencia significativa.

5.11 Conclusión del análisis

Del análisis descriptivo se concluye que las variables temporales no son influyentes en la accidentalidad, para el mes hay variaciones significativas en el número total de accidentes, por su parte para la clase de accidente, gravedad y diseño también varía el total de accidentes. Principalmente la variabilidad a nivel de accidentalidad depende en mayor medida de barrio y comuna

6 Modelos Predictivos

Se define la accidentalidad como el número de accidentes que hay a nivel diario, semanal y mensual, discriminando por tipo de accidente (CLASE) y barrio o comuna donde éste ocurre. Lo anterior debido a que existe una variabilidad notable a nivel de accidentalidad para cada columna y para cada barrio como se observó en el análisis descriptivo.

Dado que el objetivo es predecir la accidentalidad a nivel diario, mensual o semanal, se busca construir modelos que predigan el número de accidentes en cada uno de los rangos temporales es decir se crearán modelos específicos para cada comuna/barrio según la unidad temporal.

Se piensa que los modelos mixtos o modelos jerárquicos pueden representar bien esta definición de accidentalidad, debido a la agrupación que existe al momento de definir la accidentalidad, es decir, existe una agrupación entre comuna y clase; y otra agrupación entre barrio y clase. Por lo que se decide crear estos modelos con la ayuda del paquete lme4 y medir su ajuste. En cada uno de estos casos se crea un modelo que prediga el número de accidentes, por lo que se considera que la accidentalidad puede seguir una distribución Poisson.

Para evaluar los modelos se crean conjuntos de entrenamiento y conjuntos de validación. Los conjuntos de entrenamiento constan de todas las combinaciones posibles de tipos de accidentes según si es comuna o barrio y su unidad de tiempo con los años (periodos) del 2014 al 2017, el conjunto de validación será con el año 2018.

La medida para evaluar el ajuste de los modelos será el error cuadrático medio en las predicciones tanto para los conjuntos de entrenamiento como para los de validación.

\[MSE = \frac{\sum_{i=1}^N (accidentalidad_i - \widehat{accidentalidad_i)^2}}{N}\]

Con la ayuda de R se crea la siguiente función para obtenerlo:

MSE <- function(y, y_est) mean((y-y_est)**2)

Se usarán los siguientes paquetes para ajustar los modelos

library(lme4)      # Paquete para la creación de modelos mixtos poisson
library(tidyverse) # Paquete para la creación de los conjuntos de datos

6.1 Ajuste de la base de datos

Se debe modificar la base de datos depurada para poder realizar los modelos. Se crean variables de tiempo a nivel diario, semanal y mensual, además de la variable accidentalidad y días festivos.

#  Vector con fechas de dias festivos
  # Definir dias festivos
  festivos <- ymd(c(
    #2014
    '2014-01-01','2014-01-06','2014-03-24','2014-04-17','2014-04-18', '2014-05-01',
    '2014-06-02','2014-06-23','2014-06-30','2014-07-20','2014-08-07', '2014-08-18',
    '2014-10-13','2014-11-03','2014-11-17','2014-12-08','2014-12-25',
    #2015
    '2015-01-01','2015-01-12','2015-03-23','2015-03-29','2015-04-02','2015-04-03',
    '2015-04-05','2015-05-01','2015-05-18','2015-06-08','2015-06-15','2015-06-29',
    '2015-07-20','2015-08-07', '2015-08-17','2015-10-12','2015-11-02','2015-11-16',
    '2015-12-08','2015-12-25',
    #2016
    '2016-01-01','2016-01-11','2016-03-20','2016-03-21','2016-03-24','2016-03-25',
    '2016-03-27','2016-05-01','2016-05-09','2016-05-30','2016-06-06','2016-07-04',
    '2016-07-20','2016-08-07', '2016-08-15','2016-10-17','2016-11-07','2016-11-14',
    '2016-12-08','2016-12-25',
    #2017
    '2017-01-01','2017-01-09','2017-03-20','2017-04-09','2017-04-13','2017-04-14',
    '2017-04-16','2017-05-01','2017-05-29','2017-06-19','2017-06-26','2017-07-03',
    '2017-07-20','2017-08-07', '2017-08-21','2017-10-16','2017-11-06','2017-11-13',
    '2017-12-08','2017-12-25',
    #2018
    '2018-01-01','2018-01-08','2018-03-19','2018-03-25','2018-03-29','2018-03-30',
    '2018-04-01','2018-05-01','2018-05-14','2018-06-04','2018-06-11','2018-07-02',
    '2018-07-20','2018-08-07', '2018-08-20','2018-10-15','2018-11-05','2018-11-11',
    '2018-12-08','2018-12-25',
    #2019
    '2019-01-01','2019-01-07','2019-03-25','2019-04-18','2019-04-19','2019-05-1',
    '2019-06-03','2019-06-24','2019-07-01','2019-07-20','2019-08-07','2019-08-19',
    '2019-10-14','2019-11-04','2019-11-11','2019-12-08','2019-12-25'
    ))

df <- df %>% select(COMUNA, CLASE, FECHA, PERIODO, MES, DIA_NOMBRE)

df$FECHA <- as.Date(df$FECHA)

# Se agrega la variable TIEMPO y MES_NOMBRE
df$MES_NOMBRE <- paste(df$PERIODO, df$MES, sep="-") %>% as.yearmon("%Y-%m")

# Para obtener la inversa se usaría: zoo::as.Date(df$FECHA, origin="2014-01-01")
df$TIEMPO_DIA <- as.numeric(as.Date(df$FECHA)) - as.numeric(as.Date("2014-01-01")) + 1

# Se crea la variable SEMANA
df <- df %>% mutate(SEMANA = strftime(FECHA, format = "%Y-%V"),
                          TIEMPO_SEMANA = match(SEMANA, sort(unique(SEMANA))))

# Se crea la variable MES
df <- df %>% mutate(MES = strftime(FECHA, format = "%Y-%m"),
                          TIEMPO_MES = match(MES, sort(unique(MES))))

# días festivos
df <- df %>% mutate(DIA_FESTIVO = ifelse(ymd(FECHA) %in% festivos,1,0))

# accidentalidad
df <- df %>% mutate(ACCIDENTALIDAD = 1)

6.2 Modelos predictivos para las Comunas

Para obtener la base datos en la cual se encuentre el número de accidentes para cada evento considerado, se crea una base de datos con todos los eventos posibles con ayuda de la función expand.grid, en la cual se tiene en cuenta todas la combinaciones de comuna, clase (tipo de accidente) y fecha.

fecha_vector <- as.Date(as.Date("2014-01-01"):as.Date("2018-12-31"))
base <- expand.grid(COMUNA = levels(df$COMUNA), CLASE = levels(df$CLASE),
                             FECHA = fecha_vector)
base <- base %>% mutate(TIEMPO_DIA = as.numeric(FECHA) -
                                            as.numeric(as.Date("2014-01-01")) + 1)

# PERIODO
base <- base %>% mutate(PERIODO = as.numeric(format(FECHA,'%Y')))

# Se crea la variable SEMANA
base <- base %>% mutate(SEMANA = strftime(FECHA, format = "%Y-%V"),
                          TIEMPO_SEMANA = match(SEMANA, sort(unique(SEMANA))))

# Se crea la variable MES
base <- base %>% mutate(MES = strftime(FECHA, format = "%Y-%m"),
                          TIEMPO_MES = match(MES, sort(unique(MES))))

# días festivos
base <- base %>% mutate(DIA_FESTIVO = ifelse(ymd(FECHA) %in% festivos,1,0))

Luefo se realiza un left join con la base que contiene todos los posibles eventos y la base de datos que contiene la información de los accidentes (base depurada).

base <- left_join(base, subset(df, select = -DIA_NOMBRE),
                  by = c("COMUNA", "CLASE", "FECHA", "TIEMPO_DIA",
                                       "PERIODO", "SEMANA", "TIEMPO_SEMANA",
                                       "TIEMPO_MES", "DIA_FESTIVO"))
base[is.na(base)] <- 0

accidentes_dia_comuna <- left_join(base, distinct(df[, c("FECHA", "DIA_NOMBRE")]), by = "FECHA")

6.2.1 Modelos - unidad de tiempo Día

La base de datos con el número de accidentes por comuna, clase y día se obtuvo anteriormente, ésta contiene 175296 observaciones.

dim(accidentes_dia_comuna)

## [1] 175296     12

La base luce de la siguiente manera:

Se realiza la partición de los datos para seleccionar el conjunto de entrenamiento y el conjunto de prueba.

test_dia_comuna <- accidentes_dia_comuna[accidentes_dia_comuna$PERIODO == 2018, ]
train_dia_comuna <- accidentes_dia_comuna[accidentes_dia_comuna$PERIODO %in%
                                            c(2014, 2015, 2016, 2017), ]

Se procede a la creación de modelos:

Se consideran distintos modelos partiendo de un modelo jerárquico con 2 niveles (CLASE y COMUNA) con solo un intercepto, ambos niveles o variables con efectos aleatorios se consideran correlacionados y se pondrá la variable CLASE (tipo de accidente) dentro de la variable comuna.

mod_dia_comuna0 <- glmer(ACCIDENTALIDAD ~ 1 + (1 | COMUNA/CLASE),
                         data = train_dia_comuna, family= poisson())
mod_dia_comuna1 <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | COMUNA/CLASE),
                         data = train_dia_comuna, family= poisson())
mod_dia_comuna2 <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + TIEMPO_DIA + (1 | COMUNA/CLASE),
                         data = train_dia_comuna, family= poisson())
mod_dia_comuna3 <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + DIA_NOMBRE + (1 | COMUNA/CLASE),
                         data = train_dia_comuna, family= poisson())
mod_dia_comuna4 <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + DIA_NOMBRE +
                           (1 + DIA_FESTIVO| COMUNA/CLASE),
                         data = train_dia_comuna, family= poisson())

Se comparan los modelos por un analísis de varianza ANOVA:

anova(mod_dia_comuna0, mod_dia_comuna1, mod_dia_comuna2, mod_dia_comuna3, mod_dia_comuna4)

Al parecer cada modelo explica mejor la accidentalidad que el anterior, pero se debe recordar que el Valor-P es sensible al número de datos. En la base de entrenamiento hay 140256 observaciones, esto se considera grande, por lo que no hay que fijarse mucho en el Valor-P dado que es más fácil que rechace las hipótesis nulas.

Se procede a calcular el MSE de entrenamiento y de prueba en cada uno de los modelos, cabe resaltar que los valores predichos se redondearán al entero más cercano debido a que son datos de conteo y la función predict devuelve valores en la escala original pero con décimales.

# Modelo 0
y_est_train_dia_comuna0 <- round(predict(mod_dia_comuna0, newdata = train_dia_comuna,
                                         type = "response"),0)
y_est_test_dia_comuna0 <- round(predict(mod_dia_comuna0, newdata = test_dia_comuna,
                                        type = "response"),0)
mse_train_dia_comuna0 <- round(MSE(train_dia_comuna$ACCIDENTALIDAD,
                                   y_est_train_dia_comuna0),4)
mse_test_dia_comuna0 <- round(MSE(test_dia_comuna$ACCIDENTALIDAD,
                                 y_est_test_dia_comuna0),4)

# Modelo 1
y_est_train_dia_comuna1 <- round(predict(mod_dia_comuna1, newdata = train_dia_comuna,
                                         type = "response"),0)
y_est_test_dia_comuna1 <- round(predict(mod_dia_comuna1, newdata = test_dia_comuna,
                                        type = "response"),0)
mse_train_dia_comuna1 <- round(MSE(train_dia_comuna$ACCIDENTALIDAD,
                                   y_est_train_dia_comuna1),4)
mse_test_dia_comuna1 <- round(MSE(test_dia_comuna$ACCIDENTALIDAD,
                                  y_est_test_dia_comuna1),4)

# Modelo 2
y_est_train_dia_comuna2 <- round(predict(mod_dia_comuna2, newdata = train_dia_comuna,
                                         type = "response"),0)
y_est_test_dia_comuna2 <- round(predict(mod_dia_comuna2, newdata = test_dia_comuna,
                                        type = "response"),0)
mse_train_dia_comuna2 <- round(MSE(train_dia_comuna$ACCIDENTALIDAD,
                                   y_est_train_dia_comuna2),4)
mse_test_dia_comuna2 <- round(MSE(test_dia_comuna$ACCIDENTALIDAD,
                                  y_est_test_dia_comuna2),4)

# Modelo 3
y_est_train_dia_comuna3 <- round(predict(mod_dia_comuna3, newdata = train_dia_comuna,
                                         type = "response"),0)
y_est_test_dia_comuna3 <- round(predict(mod_dia_comuna3, newdata = test_dia_comuna,
                                        type = "response"),0)
mse_train_dia_comuna3 <- round(MSE(train_dia_comuna$ACCIDENTALIDAD,
                                   y_est_train_dia_comuna3),4)
mse_test_dia_comuna3 <- round(MSE(test_dia_comuna$ACCIDENTALIDAD,
                                  y_est_test_dia_comuna3),4)

# Modelo 4
y_est_train_dia_comuna4 <- round(predict(mod_dia_comuna4, newdata = train_dia_comuna,
                                         type = "response"),0)
y_est_test_dia_comuna4 <- round(predict(mod_dia_comuna4, newdata = test_dia_comuna,
                                        type = "response"),0)
mse_train_dia_comuna4 <- round(MSE(train_dia_comuna$ACCIDENTALIDAD,
                                   y_est_train_dia_comuna4),4)
mse_test_dia_comuna4 <- round(MSE(test_dia_comuna$ACCIDENTALIDAD,
                                  y_est_test_dia_comuna4),4)

Modelo_Mixto_Poisson	Train_MSE	Test_MSE	Porcentaje_Variacion
ACCIDENTALIDAD ~ 1 + (1 | COMUNA/CLASE)	2.0222	1.9352	2.20
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | COMUNA/CLASE)	1.9038	1.8241	2.14
ACCIDENTALIDAD ~ DIA_FESTIVO + TIEMPO_DIA + (1 | COMUNA/CLASE)	1.9010	1.8598	1.10
ACCIDENTALIDAD ~ DIA_FESTIVO + DIA_NOMBRE + (1 | COMUNA/CLASE)	1.5791	1.5311	1.54
ACCIDENTALIDAD ~ DIA_FESTIVO + DIA_NOMBRE + (1 + DIA_FESTIVO | COMUNA/CLASE)	1.5400	1.5118	0.92

En los modelos mod_dia_comuna3 y mod_dia_comuna4 se descarta usar la variable TIEMPO_DIA debido a que en lugar en disminuir el Test_MSE, lo aumentó, y tampoco mejoró el Train_MSE significativamente, por lo que se considera que no es significativa para explicar la accidentalidad en las comunas de Medellín por tipo de accidente.

Finalmente se observa que el modelo con intercepto y pendiente aleatoria mod_dia_comuna4:

\[\widehat{\text{ACCIDENTALIDAD}} = \text{DIA_FESTIVO} + \text{DIA_NOMBRE} + (1 + \text{DIA_FESTIVO} | \text{COMUNA/CLASE})\]

Es el que obtiene el menor Train_MSE y Test_MSE, por lo que se escoge este modelo para la explicación de la accidentalidad en comuna por día.

A continuación se muestran algunos de los valores observados de la accidentalidad diaria por clase y por comuna, con su respectivo valor predicho última columna) en el conjunto de prueba.

6.2.2 Modelos - unidad de tiempo Semana

Para este modelo se agrupa la unidad de tiempo por semanas, en las cuales se suma la cantidad de accidentes que suceden el la misma, también se suma el número de días festivos. La base luce como sigue:

accidentes_semana_comuna <- accidentes_dia_comuna %>%
  group_by(COMUNA, CLASE, PERIODO, SEMANA, TIEMPO_SEMANA) %>%
  summarise(DIA_FESTIVO = sum(DIA_FESTIVO), ACCIDENTALIDAD = sum(ACCIDENTALIDAD))

accidentes_semana_comuna

Se procede a realizar la partición del conjunto de entrenamiento y el conjunto de prueba:

test_semana_comuna <- accidentes_semana_comuna[accidentes_semana_comuna$PERIODO == 2018, ]
train_semana_comuna <- accidentes_semana_comuna[accidentes_semana_comuna$PERIODO %in%
                                                  c(2014, 2015, 2016, 2017), ]

Ahora se ajustan los modelos mixtos con efectos aleatorios en COMUNA y CLASE al igual que se hizo a nivel diario, es decir, CLASE dentro de COMUNA.

mod_semana_comuna0 <- glmer(ACCIDENTALIDAD ~  1 + (1 | COMUNA/CLASE),
               data = train_semana_comuna, family= poisson())
mod_semana_comuna1 <- glmer(ACCIDENTALIDAD ~  DIA_FESTIVO + (1 | COMUNA/CLASE),
               data = train_semana_comuna, family= poisson())
mod_semana_comuna2 <- glmer(ACCIDENTALIDAD ~  DIA_FESTIVO + TIEMPO_SEMANA +
                              (1 | COMUNA/CLASE),
               data = train_semana_comuna, family= poisson())
mod_semana_comuna3 <- glmer(ACCIDENTALIDAD ~  DIA_FESTIVO +
                              (1 + DIA_FESTIVO | COMUNA/CLASE),
               data = train_semana_comuna, family= poisson())

Se comparan los modelos por un analísis de varianza ANOVA:

anova(mod_semana_comuna0, mod_semana_comuna1, mod_semana_comuna2, mod_semana_comuna3)

Se tiene el mismo problema del tamaño de muestra que en el modelo diario. Ahora, se procede a calcular el MSE de entrenamiento y de prueba en cada uno de los modelos, cabe resaltar que los valores predichos se redondearán al entero más cercano debido a que son datos de conteo y la función predict devuelve valores en la escala original pero con décimales.

# Modelo 0
y_est_train_semana_comuna0 <- round(predict(mod_semana_comuna0, newdata = train_semana_comuna,
                                         type = "response"),0)
y_est_test_semana_comuna0 <- round(predict(mod_semana_comuna0, newdata = test_semana_comuna,
                                        type = "response"),0)
mse_train_semana_comuna0 <- round(MSE(train_semana_comuna$ACCIDENTALIDAD,
                                   y_est_train_semana_comuna0),4)
mse_test_semana_comuna0 <- round(MSE(test_semana_comuna$ACCIDENTALIDAD,
                                 y_est_test_semana_comuna0),4)

# Modelo 1
y_est_train_semana_comuna1 <- round(predict(mod_semana_comuna1, newdata = train_semana_comuna,
                                         type = "response"),0)
y_est_test_semana_comuna1 <- round(predict(mod_semana_comuna1, newdata = test_semana_comuna,
                                        type = "response"),0)
mse_train_semana_comuna1 <- round(MSE(train_semana_comuna$ACCIDENTALIDAD,
                                   y_est_train_semana_comuna1),4)
mse_test_semana_comuna1 <- round(MSE(test_semana_comuna$ACCIDENTALIDAD,
                                  y_est_test_semana_comuna1),4)

# Modelo 2
y_est_train_semana_comuna2 <- round(predict(mod_semana_comuna2, newdata = train_semana_comuna,
                                         type = "response"),0)
y_est_test_semana_comuna2 <- round(predict(mod_semana_comuna2, newdata = test_semana_comuna,
                                        type = "response"),0)
mse_train_semana_comuna2 <- round(MSE(train_semana_comuna$ACCIDENTALIDAD,
                                   y_est_train_semana_comuna2),4)
mse_test_semana_comuna2 <- round(MSE(test_semana_comuna$ACCIDENTALIDAD,
                                  y_est_test_semana_comuna2),4)

# Modelo 3
y_est_train_semana_comuna3 <- round(predict(mod_semana_comuna3, newdata = train_semana_comuna,
                                         type = "response"),0)
y_est_test_semana_comuna3 <- round(predict(mod_semana_comuna3, newdata = test_semana_comuna,
                                        type = "response"),0)
mse_train_semana_comuna3 <- round(MSE(train_semana_comuna$ACCIDENTALIDAD,
                                   y_est_train_semana_comuna3),4)
mse_test_semana_comuna3 <- round(MSE(test_semana_comuna$ACCIDENTALIDAD,
                                  y_est_test_semana_comuna3),4)

Modelo_Mixto_Poisson	Train_MSE	Test_MSE	Porcentaje_Variacion
ACCIDENTALIDAD ~ 1 + (1 | COMUNA/CLASE)	16.7066	15.1418	4.91
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | COMUNA/CLASE)	15.0649	13.2131	6.55
ACCIDENTALIDAD ~ DIA_FESTIVO + TIEMPO_SEMANA + (1 | COMUNA/CLASE)	14.9353	13.7151	4.26
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 + DIA_FESTIVO | COMUNA/CLASE)	14.9111	13.1024	6.46

En el modelo mod_semana_comuna3 se descarta usar la variable TIEMPO_SEMANA debido a que en lugar en disminuir el Test_MSE, lo aumentó, y tampoco mejoró el Train_MSE significativamente, por lo que se considera que no es significativa para explicar la accidentalidad en las comunas de Medellín por tipo de accidente.

Finalmente se observa que el modelo con intercepto y pendiente aleatoria mod_semana_comuna3:

\[\widehat{\text{ACCIDENTALIDAD}} = \text{DIA_FESTIVO} + (1 + \text{DIA_FESTIVO} | \text{COMUNA/CLASE})\]

Es el que obtiene el menor Train_MSE y Test_MSE, por lo que se escoge este modelo para la explicación de la accidentalidad en comuna por semana.

A continuación se muestran algunos de los valores observados de la accidentalidad semanal por clase y por comuna, con su respectivo valor predicho (última columna) en el conjunto de prueba.

6.2.3 Modelos - unidad de tiempo Mes

Para este modelo se agrupa la unidad de tiempo por meses, en los cuales se suma la cantidad de accidentes que suceden el el mismo, también se suma el número de días festivos que hay en dicho mes. La base luce como sigue:

accidentes_mes_comuna <- accidentes_dia_comuna %>% group_by(COMUNA, CLASE, PERIODO, MES, TIEMPO_MES) %>%
  summarise(DIA_FESTIVO = sum(DIA_FESTIVO), ACCIDENTALIDAD = sum(ACCIDENTALIDAD))

accidentes_mes_comuna

Se realiza la partición de la base de entrenamiento y de prueba

test_mes_comuna <- accidentes_mes_comuna[accidentes_mes_comuna$PERIODO == 2018, ]
train_mes_comuna <- accidentes_mes_comuna[accidentes_mes_comuna$PERIODO %in%
                                            c(2014, 2015, 2016, 2017), ]

Ahora se ajustan los modelos mixtos con efectos aleatorios en COMUNA y CLASE al igual que se hizo a nivel diario y semanal, es decir, CLASE dentro de COMUNA.

mod_mes_comuna0 <- glmer(ACCIDENTALIDAD ~  1 + (1 | COMUNA/CLASE),
               data = train_mes_comuna, family= poisson())
mod_mes_comuna1 <- glmer(ACCIDENTALIDAD ~  DIA_FESTIVO + (1 | COMUNA/CLASE),
               data = train_mes_comuna, family= poisson())
mod_mes_comuna2 <- glmer(ACCIDENTALIDAD ~  DIA_FESTIVO + (1 + DIA_FESTIVO | COMUNA/CLASE),
               data = train_mes_comuna, family= poisson())

Análisis de varianza ANOVA

anova(mod_mes_comuna0, mod_mes_comuna1, mod_mes_comuna2)

Se observa que el modelo mod_mes_comuna2, el cual tiene intercepto y pendiente aleatoria con la característica DIA_FESTIVO no aporta significativamente a la explicación de la accidentalidad según el análisis de varianza ANOVA. Sin embargo, se procede a calcular los MSE de entrenamiento y de prueba.

# Modelo 0
y_est_train_mes_comuna0 <- round(predict(mod_mes_comuna0, newdata = train_mes_comuna,
                                         type = "response"),0)
y_est_test_mes_comuna0 <- round(predict(mod_mes_comuna0, newdata = test_mes_comuna,
                                        type = "response"),0)
mse_train_mes_comuna0 <- round(MSE(train_mes_comuna$ACCIDENTALIDAD,
                                   y_est_train_mes_comuna0),4)
mse_test_mes_comuna0 <- round(MSE(test_mes_comuna$ACCIDENTALIDAD,
                                 y_est_test_mes_comuna0),4)

# Modelo 1
y_est_train_mes_comuna1 <- round(predict(mod_mes_comuna1, newdata = train_mes_comuna,
                                         type = "response"),0)
y_est_test_mes_comuna1 <- round(predict(mod_mes_comuna1, newdata = test_mes_comuna,
                                        type = "response"),0)
mse_train_mes_comuna1 <- round(MSE(train_mes_comuna$ACCIDENTALIDAD,
                                   y_est_train_mes_comuna1),4)
mse_test_mes_comuna1 <- round(MSE(test_mes_comuna$ACCIDENTALIDAD,
                                  y_est_test_mes_comuna1),4)

# Modelo 2
y_est_train_mes_comuna2 <- round(predict(mod_mes_comuna2, newdata = train_mes_comuna,
                                         type = "response"),0)
y_est_test_mes_comuna2 <- round(predict(mod_mes_comuna2, newdata = test_mes_comuna,
                                        type = "response"),0)
mse_train_mes_comuna2 <- round(MSE(train_mes_comuna$ACCIDENTALIDAD,
                                   y_est_train_mes_comuna2),4)
mse_test_mes_comuna2 <- round(MSE(test_mes_comuna$ACCIDENTALIDAD,
                                  y_est_test_mes_comuna2),4)

Modelo_Mixto_Poisson	Train_MSE	Test_MSE	Porcentaje_Variacion
ACCIDENTALIDAD ~ 1 + (1 | COMUNA/CLASE)	88.9306	94.6510	3.12
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | COMUNA/CLASE)	87.6868	91.3359	2.04
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 + DIA_FESTIVO | COMUNA/CLASE)	87.6793	91.0191	1.87

Aunque el modelo mod_mes_comuna2 no tiene una mejora muy significativa respecto al modelo mod_mes_comuna1, sigue siendo mejor con el críterio del Test_MSE, por lo que el modelo elegido es el modelo mod_mes_comuna2

\[\widehat{\text{ACCIDENTALIDAD}} = \text{DIA_FESTIVO} + (1 + \text{DIA_FESTIVO} | \text{COMUNA/CLASE})\] A continuación se enseña la comparación entre la accidentalidad observada y la accidentalidad predicha mensualmente para la base de prueba (columna final).

data.frame(test_mes_comuna, A.PREDICHA = y_est_test_mes_comuna2)

6.3 Modelos predictivos para los Barrios

La base de datos para realizar el ajuste de los barrios se obtuvo igual a como se realizó para las comunas, es decir, con un expand.grid pero esta vez agrupando la base de datos por barrios.

Realizando el mismo método visto anteriormente para la selección de modelos en las comunas, se encuentran los siguiente modelos para los barrios.

6.3.1 Modelo - unidad de tiempo Día

Base de datos agrupada para la accidentalidad (número de accidentes) diaria.

accidentes_dia_barrio # Base de datos diaria

Conjunto de prueba y conjunto de validación:

test_dia_barrio <- accidentes_dia_barrio[accidentes_dia_barrio$PERIODO == 2018, ]
train_dia_barrio <- accidentes_dia_barrio[accidentes_dia_barrio$PERIODO %in%
                                            c(2014, 2015, 2016, 2017), ]

Se encontró que el modelo que mejor explica esta accidentalidad es:

mod_dia_barrio <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | BARRIO/CLASE),
               data = train_dia_barrio, family= poisson())

Con los siguientes Train_MSE y Test_MSE

y_est_train_dia_barrio <- round(predict(mod_dia_barrio, newdata = train_dia_barrio,
                                         type = "response"),0)
y_est_test_dia_barrio <- round(predict(mod_dia_barrio, newdata = test_dia_barrio,
                                        type = "response"),0)
mse_train_dia_barrio <- round(MSE(train_dia_barrio$ACCIDENTALIDAD,
                                   y_est_train_dia_barrio),4)
mse_test_dia_barrio <- round(MSE(test_dia_barrio$ACCIDENTALIDAD,
                                 y_est_test_dia_barrio),4)

Modelo_Mixto_Poisson	Train_MSE	Test_MSE	Porcentaje_Variacion
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | BARRIO/CLASE)	0.0855	0.0825	1.79

Modelo:

\[\widehat{\text{ACCIDENTALIDAD}} = \text{DIA_FESTIVO} + (1 ~|~ \text{COMUNA/CLASE})\]

6.3.2 Modelo - unidad de tiempo Semana

Base para el modelo semanal

accidentes_semana_barrio <- accidentes_dia_barrio %>% group_by(BARRIO, CLASE, PERIODO, SEMANA, TIEMPO_SEMANA) %>%
  summarise(DIA_FESTIVO = sum(DIA_FESTIVO), ACCIDENTALIDAD = sum(ACCIDENTALIDAD))

accidentes_semana_barrio # Base de datos semanaria

Conjunto de prueba y conjunto de validación:

test_semana_barrio <- accidentes_semana_barrio[accidentes_semana_barrio$PERIODO == 2018, ]
train_semana_barrio <- accidentes_semana_barrio[accidentes_semana_barrio$PERIODO %in% c(2014, 2015, 2016, 2017), ]

Se encontró que el modelo que mejor explica esta accidentalidad es:

mod_semana_barrio <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | BARRIO/CLASE),
               data = train_semana_barrio, family= poisson())

y_est_train_semana_barrio <- round(predict(mod_semana_barrio, newdata = train_semana_barrio,
                                         type = "response"),0)
y_est_test_semana_barrio <- round(predict(mod_semana_barrio, newdata = test_semana_barrio,
                                        type = "response"),0)
mse_train_semana_barrio <- round(MSE(train_semana_barrio$ACCIDENTALIDAD,
                                   y_est_train_semana_barrio),4)
mse_test_semana_barrio <- round(MSE(test_semana_barrio$ACCIDENTALIDAD,
                                 y_est_test_semana_barrio),4)

Modelo_Mixto_Poisson	Train_MSE	Test_MSE	Porcentaje_Variacion
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | BARRIO/CLASE)	0.5967	0.6149	1.5

Modelo:

\[\widehat{\text{ACCIDENTALIDAD}} = \text{DIA_FESTIVO} + (1 ~|~ \text{COMUNA/CLASE})\]

6.3.3 Modelo - unidad de tiempo Mes

Base para el modelo mensual

accidentes_mes_barrio <- accidentes_dia_barrio %>% group_by(BARRIO, CLASE, PERIODO, MES, TIEMPO_MES) %>%
  summarise(DIA_FESTIVO = sum(DIA_FESTIVO), ACCIDENTALIDAD = sum(ACCIDENTALIDAD))

accidentes_mes_barrio # Base de datos mes

Conjunto de prueba y conjunto de validación:

test_mes_barrio <- accidentes_mes_barrio[accidentes_mes_barrio$PERIODO == 2018, ]
train_mes_barrio <- accidentes_mes_barrio[accidentes_mes_barrio$PERIODO %in% c(2014, 2015, 2016, 2017), ]

Se encontró que el modelo que mejor explica esta accidentalidad es:

mod_mes_barrio <- glmer(ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | BARRIO/CLASE),
               data = train_mes_barrio, family= poisson())

y_est_train_mes_barrio <- round(predict(mod_mes_barrio, newdata = train_mes_barrio,
                                         type = "response"),0)
y_est_test_mes_barrio <- round(predict(mod_mes_barrio, newdata = test_mes_barrio,
                                        type = "response"),0)
mse_train_mes_barrio <- round(MSE(train_mes_barrio$ACCIDENTALIDAD,
                                   y_est_train_mes_barrio),4)
mse_test_mes_barrio <- round(MSE(test_mes_barrio$ACCIDENTALIDAD,
                                 y_est_test_mes_barrio),4)

Modelo_Mixto_Poisson	Train_MSE	Test_MSE	Porcentaje_Variacion
ACCIDENTALIDAD ~ DIA_FESTIVO + (1 | BARRIO/CLASE)	2.8034	3.6059	12.52

Modelo:

\[\widehat{\text{ACCIDENTALIDAD}} = \text{DIA_FESTIVO} + (1 ~|~ \text{COMUNA/CLASE})\]

A continuación se enseña la comparación entre la accidentalidad observada y la accidentalidad predicha mensualmente para la base de prueba (columna final).

data.frame(test_mes_barrio, A.PREDICHA = y_est_test_mes_barrio)

remove(mod_dia_barrio)

7 Clustering

Según lo descrito en la introducción, es necesario centrarse en los métodos preventivos más que en las sanciones, debido a que, si la población conoce sus cifras y logra entender qué es lo que esta pasando dentro de la ciudad, se logrará no sólo reducir los incidentes de tránsito, sino también crear conciencia en las generaciones futuras.

Asimismo, para reducir los accidentes de transito en la ciudad es necesario comprender qué esta sucediendo a profundidad dentro de ella, y, gracias al acceso a las bases de datos anuales de Geomedellin, es posible clasificar las zonas donde más ocurren, para que estas sean intervenidas de una forma más rigurosa. En esta sección del proyecto fue necesario plantearnos la siguiente pregunta: ¿Cómo podríamos clasificar los barrios de Medellín de acuerdo a su accidentalidad?.

Considerando las caracteristicas de Medellín y del valle de Aburrá en general, pueden existir muchos factores que catalogan un barrio, en cuanto a seguridad vial, como peligroso o seguro. Día a día en la ciudad transitan todo tipo de vehículos, que dadas las condiciones de la zona, las calles angostas en algunos sectores, las vías en estado de deterioro y la imprudencia de la gente, hacen que sea muy complicado ir de un lugar a otro sin tener algún tipo de altercado en algún momento de la vida de una persona. Por este motivo, el criterio escogido para clasificar un barrio en vista de su accidentalidad, con los datos brindados de geomedellin, fue la gravedad del accidente. Adicionalmente, al ser solo 3 tipos de gravedad es facil de representar el análisis de forma más gráfica.

En primer lugar, la gravedad del accidente habla mucho de las condiciones en las que se dio el mismo, ya que muchas veces un mayor daño es causado por un mayor grado de irresponsabilidad. En consecuencia, es posible inferir que los barrios en los que se muestran los peores accidentes tengan la mayor cantidad de estos en general (aplica para el caso contrario). De manera anólaga, en el caso de que lo anterior no se cumpla, con este dato por lo menos se podrá identificar los barrios con los accidentes más graves, que del mismo modo, servirá para prevenir a la ciudadanía y tratar de identificar de manera profunda lo que está pasando en el lugar. Al fin y al cabo reducir los daños materiales es importante, pero reducir la cantidad de muertos en Medellín debe ser la prioridad.

En cuanto al clustering desarrollado para la clasificación de los barrios, se utilizó el método de K-Means. Aquí, se tomó la decisión de basarnos en la columna de la gravedad, ya que, no solo mostraría los barrios con más accidentalidad, sino el peligro que estos representan en cuanto a seguridad vial.

Primero, fue necesario extraer las columnas BARRIO y GRAVEDAD de la base de datos, para así solo centrarnos en lo importante durante el proceso.

Posteriormente, fue necesario hacer 2 pasos: primero se crearon 3 columnas correspondientes a los tipos de gravedad, asignando un valor de uno en la fila correspondiente al mismo tipo y un valor de cero para las demás; despues, fue necesario agrupar las columnas por barrio, para obtener la suma total de accidentes por barrio, dando como resultado la tabla que se enseña a continuación

Más adelante y ahora con sólo los datos necesarios, se realizó el método del codo para encontrar el número de clusters, al mismo tiempos que se escogía la mejor división utilizando el método k-means++ que posee python. Todo esto arrojó que lo mejor sería dividir los datos en 4 clusters.

Al realizar la respectiva clasificación se arrojó la siguiente gráfica, donde se puede observar la distribución de los barrios y permite aceptar la siguiente premisa: los 3 tipos de gravedad se relacionan entre sí, evidenciandose que a mayor numero de heridos, también hay mayor número de solo daños y muertos en promedio. En ese orden de ideas, se muestra otra información importante: Unos cúantos barrios son los que presentan la mayor concentración de accidentes, con una frecuencia mayor con respecto a los más seguros.

Teniendo en cuenta lo anterior, se analizó la media de accidentes por cluster, lo que arrojó una diferencia significativa entre ellos a la hora de analizar los heridos, muertos y daños ocasionados. El heatmap que se encuentra a continuación describe lo visto en las gráficas anteriores de una manera mas precisa. Por lo tanto, los números muestran una diferencia mucho mayor cuando se comparan los 2 grupos más seguros con los 2 más peligrosos, donde visualmente se puede notar al igual que en los gráficos anteriores. De ahí que, es importnte resaltar que en las 3 componentes, el promedio del cluster con más accidentes es mayor a la suma de los otros 3

Estos datos se muestran en la aplicación de shiny, en la parte de abajo de la sección de agrupamiento.

Al contar el número de barrios, encontramos otro dato muy interesante: Además de que el número de barrios cambia mucho de un cluster a otro, el cluster que tiene una mayor accidentalidad, que en promedio equivale a casi la mitad de todos los accidentes, solo se compone de 9 barrios.

Al tener esta información y con el fin de reducir considerablemente estos incidentes viales en Medellín, es necesario prevenir a los cuidadanos de los siguientes barrios, a su vez la alcaldía plantea posibles soluciones en la movilidad de estos:

Teniendo en cuenta la información anterior y teniendolos de manera gráfica en el mapa, la mayoría de los problemas de movilidad de la ciudad, como era de esperarse, ocurren en zonas centrales, especificamente en barrios como Caribe y La Candelaria.

Al igual que la información del heatmap, estos datos se muestran en la aplicación de shiny, debajo de la sección de agrupamiento.

Para finalizar esta sección, se agregó la tabla que clasifica el riesgo del barrio en la base de datos pricipal, para proximamente mostrar esta distribución en el mapa de Medellín asignando colores a los clusters. La columna quedó con el nombre SEGURIDAD_BARRIO.

8 Conclusiones

El modelo mixto o modelo jerárquico demuestra ser una buena alternativa para predecir la accidentalidad en las comuna y barrios de Medellín en las diferentes unidades de tiempo planteadas en este trabajo (día, semana, mes), ya que presenta un Train_MSE y un Test_MSE relativamente pequeño.

Al incluir todos los casos posibles en la construcción de los modelos, es decir, incluso los eventos donde no ocurren accidentes, los modelos mixtos se hacen muy pesados computacionalmente, y las bases de datos consideras crecen en observaciones significativamente, por lo que se recomienda crear una alternativa de optimización para el uso de los mismos.

El clustering realizado para el agrupamiento de mapas y el criterio elegido para agruparlos (número de accidentes) demuestran ser bien elegidos debido a que la forma en que se agrupan los datos tiene sentido con la realidad de la ciudad (La Candelaria y Laureles son barrios con mucho flujo vehícular y se encuentran en zonas donde hay mucha afluencia de personas de todo el Valle de Aburrá)

En los meses de agosto se puede observar que la accidentalidad es mucho mas alta que en los demás meses, algo a destacar es que el mes de agosto en los años analizados se festejaba la feria de flores, esto lleva a que la cantidad de personas en la ciudad de Medellín sea mas alta y por ende mas probabilidad de accidentes abra, ya que la ciudad estará mas congestionada, en cambio que en los meses de diciembre y enero la ciudad esta un poco mas descongestionada esto causa una posible disminución en la accidentalidad de Medellín.

9 Recomendaciones

Crear un plan de acción inmediato en los 9 barrios que se encuentran con la categoría de “Accidentalidad alta” ya que presentan índices muy altos que pueden tener como consecuencia un incremento en la accidentalidad y muertes en las vías.

Invertir en una campaña de conciencia y responsabilidad automovilista en el centro ya que en este lugar es donde mas se presentan accidentes en Medellín.

Cuando hay eventos como la feria de flores tener planes de movilidad mas estrictas esto con el fin de disminuir la alta accidentalidad en esa festividad.

10 Referencias

Andrew Gelman, Jennifer Hill (2007). Data Analysis Using Regression and Multilevel/Hierarchical Models.

Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.

Wickham et al., (2019). Welcome to the tidyverse. Journal of Open Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686

Achim Zeileis and Gabor Grothendieck (2005). zoo: S3 Infrastructure for Regular and Irregular Time Series. Journal of Statistical Software, 14(6), 1-27. doi:10.18637/jss.v014.i06

C. Sievert. Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC Florida, 2020.

Rstudio.github.io. 2020. Leaflet For R - Introduction. [online] Available at: https://rstudio.github.io/leaflet/ [Accessed 20 November 2020].

Tiempo, C., 2020. Estos Son Los 5 Puntos Con Más Accidentes De Tránsito En Medellín. [online] El Tiempo. Available at: https://www.eltiempo.com/colombia/medellin/estos-son-los-5-untos-con-mas-accidentes-de-transito-en-medellin-229534 [Accessed 20 November 2020].

Ctumedellin.com. 2020. Medellín, Una De Las Ciudades Más Letales Para Manejar – CTU COLOMBIA. [online] Available at: http://ctumedellin.com/sitio/2018/01/16/medellin-una-de-las-ciudades-mas-letales-para-manejar/ [Accessed 20 November 2020].

LA.Network. 2020. Estas Son Las 10 Ciudades Latinoamericanas Con Más Muertes En Las Vías. [online] Available at: https://la.network/ciudades-latinoamericanas-muertes-vias/ [Accessed 20 November 2020].