Modelación de la Accidentalidad vehicular para las 16 comunas del área urbana de Medellín entre los años 2014-2019 a través de regresión multinivel
Juan Felipe Múnera Vergara
3/21/2020
Introducción
Los accidentes de tráfico o accidentes de tránsito son siniestros viales que perturban el funcionamiento normal de una vía, este debe involucrar a un vehículo automotor. Estos sucesos son clasificados según sus repercusiones (su gravedad: solo daños materiales, heridos, muertos, etc.) y su tipo (Choques, atropellos).
En este documento se presenta el ajuste del número de accidentes del área Urbana de la ciudad de Medellín en función del tiempo, los datos que se presentaran cumplen con las características de una serie de tiempo, por lo que metodologías aptas para el modelamiento de estas son claras candidatas a la hora de elegir un modelo para explicar el número de accidentes. Sin embargo, se pueden encontrar diferencias significativas en el número de accidentes registrados según divisiones territoriales como barrios o comunas, considerar esto significaría ajustar múltiples modelos, uno para cada división territorial que se considere.
Los modelos mixtos, jerárquicos o multinivel, son una metodología que permite crear un único modelo mientras se diferencia por una variable de agrupamiento, que usualmente representa sujetos u objetos de los que se tienen repetidos registros, permitiendo involucrar efectos aleatorios que cambian de sujeto a sujeto.
En el presente trabajo se utilizan los modelos multinivel generalizados para ajustar dos modelos: el primero un modelo con respuesta Poisson para en número total de accidentes y el segundo, un modelo Multinomial para el número de accidentes según la clase (este último es tan solo un añadido). Ambos modelos son complementarios, ya que el modelo multinomial requiere de una estimación del total de accidentes para realizar las predicciones de la cantidad de accidentes según la categoría. Se utilizara como variable de agrupamiento las comunas.
Para mejorar la legibilidad del documento se ocultaron algunas lineas de código, pero se facilita la descarga del rmd utilizado para generar el documento en la parte superior, al lado del titulo.
Carga de librerías
#CAGRA DE LIBRERIAS
## Manejo de bases de datos
library(dplyr)
library(magrittr)
## Gráficos
library(ggplot2)
library(lubridate)
## Ajuste de modelos mixtos
library(lme4)
library(brms)
## Manejo de fechas
library("zoo")
## Tablas latex
library(kableExtra)Objetivos
- Construir un modelo predictivo para el número de accidentes ocurridos en la ciudad de Medellín a nivel comunal en intervalos mensuales, con base a los datos de años del 2014-2019 utilizando la regresión multinivel (conocida también como modelos jerárquicos o mixtos), utilizando los paquetes
lmea4ybrmsde R. - Comprender la forma en que la accidentalidad se relaciona con el tiempo para las comunas de Medellín.
Conjunto de datos:
Los datos análizados en el presente documento se obtienen de la página del portal GeoMedellín de la Alcaldía de Medellín en su sección de datos abiertos. Corresponden a las bases de “Accidentalidad georreferenciada” entre los años 2014 y 2019, los elementos de dichas bases representan un accidente registrado por la Secretaría de Movilidad de la Alcaldía de Medellín y presentan una serie de covariables que sirven para categorizar el tipo de accidente, el momento de su incidencia y su ubicación georreferenciada.
Para el ajuste de los modelos se descartaron algunas de las covariables (no interesantes para modelar la accidentalidad) y se realizaron cambios al formato de los datos (las modificaciones realizadas y los pasos para construir el conjunto de datos que se utiliza aquí pueden ser consultados en detalle aquí). En síntesis, se limitó el estudio a las 16 comunas de la ciudad de Medellín (descartándose aquellos accidentes ocurridos en corregimientos), se omitieron covariables relativas a la georreferenciación y se agruparon los registros por barrios en periodos mensuales, para obtener la cantidad de accidentes ocurridos distinguidos por su clase.
Los datos definitivos poseen los siguientes atributos:
COMUNA: las 16 comunas de Medellín con niveles ‘Aranjuez’, ‘Belén’, ‘Buenos Aires’, ‘Castilla’, ‘Doce de Octubre’, ‘El Poblado’, ‘Guayabal’, ‘La América’, ‘La Candelaria’, ‘Laureles Estadio’, ‘Manrique’, ‘Popular’, ‘Robledo’, ‘San Javier’, ‘Santa Cruz’, ‘Villa Hermosa’.
BARRIO: los barrios de cada comuna con 265 niveles.
PERIODO: año en el que ocurren los accidentes.
MES: con niveles ‘January’, ‘February’, ‘March’, ‘April’, ‘May’, ‘June’, ‘July’, ‘August’, ‘September’, ‘October’, ‘November’, ‘December’.
FECHA: en formato “YY/MM”.
t: tiempo transcurrido, una conversión numérica de la FECHA para agregar en el modelo, equivalente al número de días transcurridos desde la fecha de referencia 2014/01/01 hasta el tiempo del registro. Dado que los registros representan periodos mensuales, el tiempo del registro corresponderá a los días transcurridos a la mitad del MES.
Accidentalidad: número de accidentes según tipo de categoría (también llamada CLASE).
Atropello: número de accidentes de esta categoría de accidentes.
Caida.Ocupante: número de accidentes de esta categoría de accidentes.
Choque: número de accidentes de esta categoría de accidentes.
Incendio: número de accidentes de esta categoría de accidente.
Volcamiento: número de accidentes de esta categoría de accidentes.
Otro: número de accidentes de esta categoría de accidentes.
nro.accidentes: total de accidentes, suma de los accidentes de todas las categorías. Variable respuesta.
Carga del conjunto de datos
A continuación se presenta como cargar el conjunto de datos y se construyen t y FECHA:
## Carga de conjunto de datos a nivel barrial
#data_men_all <- read.csv("Datasets/base_mensual_barrios_y_comunas.csv", encoding='UTF-8')
url <- "https://raw.githubusercontent.com/jfmra99/Modelos_Jerarquicos/master/AplicacionRPubs/Datasets/base_mensual_barrios_y_comunas.csv"
data_men_all <- read.csv(url, encoding='UTF-8')
## Verificación de covariables como factor
data_men_all[,c(2,3,4)] %<>% lapply(function(x) factor(as.character(x)))
## Creación de covariables FECHA y t, en base a las otras covariables
data_men_all$FECHA <- paste(data_men_all$PERIODO, data_men_all$MES, sep="-") %>% as.yearmon("%Y-%m")
# Modificación de MES a su nombre en inglés
data_men_all$MES <- factor(months(data_men_all$FECHA), levels =month.name[1:12])
# Para obtener la inversa se usaría: zoo::as.Date(data_men_all$FECHA, origin="2014-01-01")
data_men_all$t <- as.numeric(as.Date(data_men_all$FECHA, frac=0.5)) - as.numeric(as.Date("2014-01-01"))Los accidentes del anterior conjunto son dados a nivel barrial, por lo que los agrupamos por comunas, el objeto data_men_comunas tiene los datos con los que se ajustaran los modelos. Mientras que tidy_data_men_all_comunas es una base de apoyo para construir gráficos.
## Construcción del conjunto de datos a utilizar, conseguido por Comunas.
data_men_comunas <- data_men_all %>% group_by(COMUNA, PERIODO, MES,FECHA, t) %>%
summarise(nro.accidentes = sum(nro.accidentes),
Atropello = sum(Atropello),
Caida.Ocupante = sum(Caida.Ocupante),
Choque=sum(Choque),
Volcamiento = sum(Volcamiento),
Otro = sum(Otro))
## Organización de bases de datos en formato apto para gráficar.
tidy_data_men_all <- tidyr::gather(data_men_all, "CLASE", "nro.accidentes", 5:9)
tidy_data_men_all_comunas <- tidy_data_men_all%>% group_by(COMUNA, t, CLASE) %>%
summarise(nro.accidentes = sum(nro.accidentes))A continuación se presenta el conjunto de datos:
Se puede ver que se tienen 12 registros por cada año para cada una de las 16 comunas, ya que los datos corresponden a un periodo de 6 años se tienen \(12*16*6=1152\) registros.
Análisis descriptivo
Nuestro objetivo primordial es predecir nro.accidentes en función de las covariables temporales, por lo que realizamos el gráfico de la serie de tiempo correspondiente a los accidentes de todas las Comunas juntas:
#Accidentalidad conjunta de TODAS las Comunas:
temp <- data_men_comunas %>% group_by(FECHA) %>%
summarise(nro.accidentes = sum(nro.accidentes))
## SERIE DE TIEMPO para la
ts_monthly <- ts(temp$nro.accidentes, start=c(2014,1),frequency=12)
plot(ts_monthly)
Las series de tiempo tienen componentes de tendencia y estacionalidad. Este último corresponde a comportamientos que se repiten cada cierto tiempo en la serie, existen gráficos que evaluan que periodicidad (anual, trimestral, semestral etc.) tomar como el periodograma, pero sin necesidad de recurrir a él, del gráfico de la serie presentado anteriormente se pueden notar comportamientos marcados para algunos meses entre el 2014 y el 2018, (por ejemplo, la accidentalidad cae al mínimo principio de cada año en el mes de Enero). En el 2019 se puede intuir que se conservan los comportamientos periódicos, pero parece que ocurre un ciclo que aumento el número de accidentes en el segundo semestre del año. Debido a esto se sospecha que los datos tienen una periodicidad anual, por los que la covariable MES es candidata a ser relevanta para el modelamiento del número de accidentes.
La tendencia en cambio, describe comportamiento general de la serie (Restando los periodos). Se puede obtener una forma aproximada de la tendencia de la seria por medio de decompose:
plot(decompose(ts_monthly)$trend,ylim=c(min(ts_monthly),max(ts_monthly)), ylab = "accidentes")
Se puede ver que la tendencia varia entre 3200 y 3600 accidentes, por lo que se puede hablar de una tendencia de crecimiento/decrecimiento marcada. Para modelarla se puede recurrir polinomio del tiempo, dado que se pueden apreciar varias concavidades, un polinomio de grado 3 o 4 debería ser adecuado.
Ahora exploraremos la cantidad de accidentes según la comuna y su CLASE, a continuación se presenta una tabla con el total de accidentes de los 6 años registrados para cada una de las comunas:
| COMUNA | nro.accidentes | Atropello | Caida.Ocupante | Choque | Volcamiento | Otro |
|---|---|---|---|---|---|---|
| La Candelaria | 52792 | 5120 | 3171 | 39342 | 1253 | 3906 |
| Laureles Estadio | 28385 | 1733 | 1915 | 21407 | 708 | 2622 |
| Castilla | 25359 | 1875 | 2619 | 16757 | 919 | 3189 |
| El Poblado | 20865 | 707 | 952 | 17240 | 532 | 1434 |
| Guayabal | 18343 | 1151 | 1205 | 13618 | 632 | 1737 |
| Belén | 16608 | 1243 | 1221 | 11862 | 585 | 1697 |
| Robledo | 16550 | 1382 | 2466 | 9551 | 621 | 2530 |
| Aranjuez | 14975 | 2038 | 1485 | 9143 | 538 | 1771 |
| Buenos Aires | 9554 | 1017 | 911 | 5907 | 465 | 1254 |
| La América | 8243 | 683 | 701 | 5790 | 221 | 848 |
| Manrique | 7791 | 1641 | 964 | 3762 | 360 | 1064 |
| Doce de Octubre | 6951 | 1355 | 1325 | 2890 | 247 | 1134 |
| Villa Hermosa | 6688 | 1096 | 792 | 3602 | 314 | 884 |
| San Javier | 4435 | 757 | 606 | 2223 | 214 | 635 |
| Popular | 3662 | 1074 | 489 | 1450 | 155 | 494 |
| Santa Cruz | 3393 | 802 | 366 | 1687 | 125 | 413 |
El total de accidentes varía mucho de una comuna a otra, por lo que se sospecha que existen diferencias marcadas entre las Comunas. Esto se puede apreciar mucho mejor en el siguiente gráfico de los accidentes dado el tiempo para cada una de las 16 comunas, en este y los demás gráficos que involucren el conteo de accidentes se utilizará la escala logarítmica, para que la escala permita comparar las diferentes comunas:
Se puede ver que efectivamente existe una diferencia en la cantidad de accidentes registrados de una comuna a otra. Esto nos dice que utilizar la COMUNA como la variable de agrupamiento, en base al anterior gráfico la adición de un intercepto aleatorio parece razonable, en cambio dado que que no hay diferencias de tendencia marcadas entre las comunas, se descarta la inclusión de una pendiente aleatoria y se considera que el comportamiento es similar entre ellas.
Se puede alcanzar la misma conclusión si se gráfica el número de accidentes haciendo distinción por su clase:
Ahora veremos la densidad muestral de la cantidad de accidentes según las comunas:
## Warning: Removed 39 rows containing non-finite values (stat_density).
Se puede ver la presencia de diferentes modas, pero esto puede ser explicado por las diferencias del número de accidentes según el mes. Más adelante en el documento, se considerará que el número de accidentes de distribuye \(Poisson\), por lo que este gráfico sirve como argumento para considerar que la media esperada dependerá de covariables como el tiempo. Además, se puede ver el cambio de forma de las densidades entre las comunas.
También se presentan las densidades según la CLASE del accidente:
## Warning: Removed 39 rows containing non-finite values (stat_density).
Nuevamente se puede ver el cambio de forma de las densidades entre las comunas. Lo cual es confirmado al realizar una gráfica de barras de probabilidad:
Se puede ver que la proporción de algunas clases de accidente como los Atropellos son notablemente mayores para algunas comunas como Santa Cruz y Popular, y en general varían considerablemente entre las comunas. Lo mismo ocurre en los Choques y los accidentes con Caída de Ocupante. Mientras que otras como los volcamientos los accidentes de clase “Otro” parecen tener una proporción constante entre las comunas. La conclusión general del gráfico es que las proporciones según la clase son afectadas por la variable de agrupamiento COMUNA
Modelación y Ajuste
Modelo con respuesta Poisson
Denotaremos el número de accidentes para el tiempo \(y\) de la comuna \(j\) como \(n_{ij}\) y la modelaremos por medio de una regresión Poisson:
\[n_{ij} \sim Poisson(\lambda_{ij})\] Entonces queremos encontrar una expresión que relacione a \(\lambda_{ij}\) con las variables temporales y las comunas.
Modelos propuestos
Con base al análisis descriptivo, se tiene que la tendencia general del modelo tiene que ser modelada a través de un polinomio de alto grado para captar sus concavidades (inicialmente se considera que el grado mínimo debe ser 3). Además, dado que realmente no se notaron tendencias marcadas respecto al número de accidentes según la comuna, se descartan modelos con pendiente aleatoria y se restringe a modelo con intercepto aleatorio.
Denótese:
\[j: Comuna\; (1:16) \\ i: observación\; correspondiente\; a\; un\;mes \; entre \; 2014 \; y \; 2019 \; (1:72)\] Adicionalmente se notó una periodicidad anual marcada, por lo que para modelarla la variable MES será incluida en los modelos de la siguiente forma: \[\sum\limits_{k=1}^{11}{\delta_{k}}MES_{ijk}\;\; donde\;\; MES_{ik} = \left\{ \begin{matrix} 1 \quad si ~ la ~observación ~ij~ está~ en~la~ estación~k \\ 0 \quad~ otro~caso\end{matrix} \right\}\]
\[k = 1: Febrero, ~ k=2: Marzo, ..., k=10: Noviembre, ~ k=11: Diciembre\]
Modelo 1
Polinomio de grado 3 en función del tiempo, con estacionalidad según el MES e intercepto aleatorio:
\[\log (\lambda_{ij}) = \beta_0 + \beta_1 t + \beta_2t^2 + \beta_3t^3+\sum\limits_{k=1}^{11}{\delta_{k}}MES_{ijk} + b_{0j}\] \[b_{0j} \sim N(0, \sigma_{b0})\]
Modelo 2
Polinomio de grado 4 en función del tiempo, con estacionalidad según el MES e intercepto aleatorio:
\[\log (\lambda_{ij}) = \beta_0 + \beta_1 t + \beta_2t^2 + \beta_3 t^3+ \beta_4 t^4 +\sum\limits_{k=1}^{11}{\delta_{k}}MES_{ijk} + b_{0j}\] \[b_{0j} \sim N(0, \sigma_{b0})\]
Modelo 3
Polinomio de grado 5 en función del tiempo, con estacionalidad según el MES e intercepto aleatorio:
\[\log (\lambda_{ij}) = \beta_0 + \beta_1 t + \beta_2t^2 + \beta_3 t^3+ \beta_4 t^4+ \beta_5 t^5 +\sum\limits_{k=1}^{11}{\delta_{k}}MES_{ijk} + b_{0j}\]
\[b_{0j} \sim N(0, \sigma_{b0})\]
Modelo 4
Polinomio de grado 6 en función del tiempo, con estacionalidad según el MES e intercepto aleatorio:
\[\log (\lambda_{ij}) = \beta_0 + \beta_1 t + \beta_2t^2 + \beta_3 t^3+ \beta_4 t^4 + \beta_5 t^5 + \beta_6 t^6 +\sum\limits_{k=1}^{11}{\delta_{k}}MES_{ijk}+ b_{0j}\]
\[b_{0j} \sim N(0, \sigma_{b0})\]
Ajuste con R
Se realiza el ajuste de los modelos por medio de la función glmer de la libreria lme4, añadiendo el atributo family= poisson() para modelar la distribución de \(n_{ij}\).
## glmer attempts
##Modelo1
mod1 <- glmer(nro.accidentes ~ poly(t, 3)+MES+ (1|COMUNA),
data = data_men_comunas, family= poisson())
##Modelo2
mod2 <- glmer(nro.accidentes ~ poly(t, 4)+MES+ (1|COMUNA),
data = data_men_comunas, family= poisson())
##Modelo3
mod3 <- glmer(nro.accidentes ~ poly(t, 5)+MES+ (1|COMUNA),
data = data_men_comunas, family= poisson())
##Modelo4
mod4 <- glmer(nro.accidentes ~ poly(t, 6)+MES+ (1|COMUNA),
data = data_men_comunas, family= poisson())
# Guardado de los modelos
save(mod1, mod2, mod3, mod4, file = "Modelos/comuna_glmer.RData")Comparación de modelos
En esta sección se elegirá uno de los modelos mostrados anteriormente, probando la significancia de los parámetros añadidos por medio del test de razón de verosimilitud. Considérese dos modelos: \(ModeloA\) y \(ModeloB\), con vectores de parámetros \(\theta^A\) y \(\theta^B\) respectivamente, que cumplen \(\theta^A \subset \theta^B\). Denotemos a los parámetros del \(ModeloB\) que no están es el \(ModeloA\) como \(\theta^c\) un vector de dimensión \(p\), entonces la prueba de hipótesis asociada a la significancia de estos parámetros será:
\[H_0: \theta^c = 0_{1\times p} \quad v.s \quad H_1: \theta^c \not= 0_{1\times p}\]
Cuyo estadístico de prueba es:
\[LR = -2 \times (loglikelihood(ModeloA)- loglikelihood(ModeloB))\]
A un nivel de significancia \(\alpha\) se rechaza \(H_0\) si \(VP =P(\chi^2_p > LR) < \alpha\)
Esa prueba se dice es anti-conservativa, es decir que su \(VP\) será más pequeños de lo normal, por lo que usualmente se utiliza en lugar de \(\chi^2_p\) la distribución muestral del estadístico, la cual puede ser estimada por métodos de remuestreo. Sin embargo, dada la naturaleza de serie de tiempo que tiene el conjunto de datos presentado, utilizar una técnica de remuestreo no tiene mucho sentido, por lo que en su lugar fijaremos un nivel de significancia más riguroso y por ende más pequeño de lo usual, haciendo \(\alpha = 0.01\).
A continuación, se presentan las pruebas de hipótesis para cada comparación:
Modelo 1 v.s Modelo 2
\[H_0: \beta_4 = 0 \quad v.s \quad H_1: \beta_4 \not=0\]
anova(mod1, mod2)Con un \(VP = 2.2e-16 < 0.01\) se rechaza la hipótesis nula y se considera que \(\beta_4\) es significativa.
Modelo 2 v.s Modelo 3
\[H_0: \beta_4 = 0 \quad v.s \quad H_1: \beta_4 \not=0\]
anova(mod2, mod3)Con un \(VP = 0.01764 > 0.01\) no se puede rechazar la hipótesis nula al nivel de significancia establecido y se considera que \(\beta_4\) es no significativa.
Modelo 2 v.s Modelo 4
En este caso se espera que se obtenga el mismo resultado que en la prueba del Modelo 2 v.s Modelo 3, pero igualmente se incluye:
\[H_0: \beta_4 = \beta_5 = 0 \quad v.s \quad H_1: \beta_4 \not=0 ~ | ~ \beta_5 \not=0 \]
anova(mod2, mod4)Con un \(VP = 0.05913 > 0.01\) no se puede rechazar la hipótesis nula al nivel de significancia establecido y se considera que \(\beta_4\) y \(\beta_5\) son no significativas.
De las pruebas anteriores se encuentra que el modelo óptimo es el numero 2, por lo que se procede a realizar el análisis de sus residuales.
Evaluación de residuales
Queremos chequear si se cumplen los típicos supuestos de varianza constante y normalidad de los residuales. Pero además, se evaluará si existe falta de ajuste respecto al tiempo y el supuesto de normalidad de los interceptos aleatorios \[b_{0j} \sim N(0, \sigma_{b0})\]
Los gráficos de la primera fila, correspondientes a los residuales estandarizados contra los valores ajustados (A) y a los residuales estandarizados contra el tiempo (B), no muestran ningún indicio de heterocedasticidad. los residuales en (B) no tienen ninguna tendencia marcada y se dispersan aleatoriamente, por lo que se considera no hay falta de ajuste del modelo.
En la segunda fila encontramos los qqplot de los residuales estandarizados (C) y de los intervalos aleatorios ajustados (D). De (C) se puede ver que los cuantiles de los residuales se desvían de los cuantiles teóricos, por lo que seguirían una distribución de cola pesada. De (D) no se encuentran evidencias gráficas contundentes contra la normalidad de los interceptos aleatorios.
El summary() del modelo se presenta a continuación:
summary(res1)## data: t, r, y, w [1152x4]
## mapping: x = ~w
## scales: x, xmin, xmax, xend, xintercept, xmin_final, xmax_final, xlower, xmiddle, xupper, x0
## faceting: <ggproto object: Class FacetNull, Facet, gg>
## compute_layout: function
## draw_back: function
## draw_front: function
## draw_labels: function
## draw_panels: function
## finish_data: function
## init_scales: function
## map_data: function
## params: list
## setup_data: function
## setup_params: function
## shrink: TRUE
## train_scales: function
## vars: function
## super: <ggproto object: Class FacetNull, Facet, gg>
## -----------------------------------
## geom_density: na.rm = FALSE
## stat_density: na.rm = FALSE
## position_identityExpresión del modelo
\[\log (\widehat{\lambda_{ij}}) = 4.9180 -0.1385t + 0.0906t^2 + 0.6405t^3+ 0.9706 t^4 +0.1074MES_{Febrero}+\\0.1680MES_{marzo}+0.1287MES_{Abril}+0.1924MES_{Mayo}+0.1126MES_{Junio}+\\ 0.1855MES_{Julio}+0.2369MES_{Agosto}+ 0.2147MES_{Septiembre}+ 0.2028MES_{Octubre}+\\0.1274MES_{Noviembre}+0.1641MES_{Diciembre} + b_{0j}\]
Los valores estimados de \(b_{0j}\) se muestran a continuación:
| j | Comuna | b0 |
|---|---|---|
| 1 | Aranjuez | 0.2636229 |
| 2 | Belén | 0.3671174 |
| 3 | Buenos Aires | -0.1857351 |
| 4 | Castilla | 0.7903509 |
| 5 | Doce de Octubre | -0.5037175 |
| 6 | El Poblado | 0.5952949 |
| 7 | Guayabal | 0.4664752 |
| 8 | La América | -0.3332942 |
| 9 | La Candelaria | 1.5235809 |
| 10 | Laureles Estadio | 0.9030769 |
| 11 | Manrique | -0.3896728 |
| 12 | Popular | -1.1441775 |
| 13 | Robledo | 0.3636193 |
| 14 | San Javier | -0.9528285 |
| 15 | Santa Cruz | -1.2203906 |
| 16 | Villa Hermosa | -0.5422732 |
Y la estimación de su varianza es: \(\widehat{\sigma^2_{b0}} = 0.5765\)
Modelo Ajustado
A continuación se muestra el modelo ajustado sobre los datos reales para cada una de las 16 comunas del área urbana:
Aranjuez
Belén
Buenos Aires
Castilla
Doce de Octubre
El Poblado
Guayabal
La América
La Candelaria
Laureles Estadio
Manrique
Popular
Robledo
San Javier
Santa Cruz
Villa Hermosa
Modelo con respuesta Multinomial
Como se notó en el análisis descriptivo, existen diferencias en la cantidad de accidentes según la CLASE, por lo que adicional al modelo Poisson, se presenta el ajuste una regresión con respuesta multinomial la cual puede diferenciarse de la más conocida regresión logística Multinomial, que es un método de clasificación, la respuesta multinomial tiene por respuesta conteo en diferentes categorías. Sean:
\[o = Atropello,\;\; p = Caida.Ocupante,\;\; q = Choque,\;\; r = Volcamiento,\;\; s = Otro \\ j: Comuna\; (1:16) \\ i: observación\; correspondiente\; a\; un\;mes \; entre \; 2014 \; y \; 2019 \; (1:72)\]
Ell vector de Accidentalidad \(y_{i,j}\):
\[y_{i,j} = (y_{i,j}^o, y_{i,j}^p, y_{i,j}^q, y_{i,j}^r, y_{i,j}^s)\]
Con un total de accidentes (nro.accidentes) dado por: \[n_{ij} = y_{i,j}^o + y_{i,j}^p + y_{i,j}^q + y_{i,j}^r + y_{i,j}^s\]
Y \(\theta_{ij}\) el vector de probabilidad según la clase del accidente.
\[\theta_{ij} = (\theta_{ij}^o, \theta_{ij}^p,\theta_{ij}^q, \theta_{ij}^r, \theta_{ij}^s)\]
Para ajustar nuestro modelo consideraremos que:
\[y_{ij}|\theta_{ij} \sim Multinomial(n_{ij}, \theta_{ij})\]
\[\theta_i = Softmax(0, \mu_j^p,\mu_j^q, \mu_j^r, \mu_j^s)\]
Entonces, nuestro objetivo es modelar \(\mu_j^p,\mu_j^q, \mu_j^r, \mu_j^s\) incluyendo covariables.
Para ajustar el modelo utilizaremos el paquete brmsde R, que permite ajustar modelo multinivel considerando una gran variedad de distribuciones para la variable respuesta y tiene un gran fundamento en la estadística Bayesiana. Sin embargo, se ha de mencionar que la familia multinomial es relativamente nueva en brms, por lo que carece de muchos estadísticos comparativos y metodologías para realizar comparaciones entre modelos o probar la significancia de efectos. Un modelo podría ser elegido en base a su capacidad predictiva, realizando métodos de validación cruzada (específicos para series de tiempo), tal y como expone [3]. Sin embargo, brms por sus metodologías bayesianas para el ajuste, consume muchos recursos computacionales y consiguientemente, tarda mucho en su ajuste, por lo que calcular múltiples modelos es una tarea larga y tediosa.
Dado que el ajuste de este modelo no es el objetivo principal de este documento y es tan solo un añadido que parte de la curiosidad de quien escribe estas palabras, se ajustará un único modelo.
Ajuste del modelo
\[\mu_j^p = \beta_0^p + \beta_1^p t + \beta_2^pt^2 + \beta^p_3 t^3 +\sum\limits_{k=1}^{11}{\delta^p_{k}}MES_{ijk} + b^p_{0j}\] \[\mu_j^q = \beta_0^q + \beta_1^q t + \beta_2^qt^2 + \beta^q_3 t^3+\sum\limits_{k=1}^{11}{\delta^q_{k}}MES_{ijk} + b^q_{0j}\] \[\mu_j^r = \beta_0^r + \beta_1^r t + \beta_2^rt^2 + \beta^r_3 t^3+\sum\limits_{k=1}^{11}{\delta^r_{k}}MES_{ijk} + b^r_{0j}\] \[\mu_j^s = \beta_0^s + \beta_1^s t + \beta_2^st^2 + \beta^s_3 t^3 +\sum\limits_{k=1}^{11}{\delta^s_{k}}MES_{ijk} + b^s_{0j}\]
Con \(b^p_{0}\sim N(0, \sigma^2_{b^p_{0}})\), \(b^q_{0}\sim N(0, \sigma^2_{b^q_{0}})\), \(b^r_{0}\sim N(0, \sigma^2_{b^r_{0}})\) y \(b^s_{0}\sim N(0, \sigma^2_{b^s_{0}})\)
Código en R:
El ajuste de este modelo se realizó basado en los ejemplos de [1] y [2]. Primero se hace necesario indicarle a brms cuales son las categorías para las que va a predecir, esto se hace de esta manera:
data_men_comunas$Accidentalidad <- with(data_men_comunas, cbind(Atropello, Caida.Ocupante, Choque, Volcamiento, Otro))Luego el modelo es ajustado utilizando:
formula1 <- bf( Accidentalidad | trials(nro.accidentes) ~ poly(t, 3) + MES + (1 | COMUNA))
fit1 <- brm(formula1, data = test_Data, family = multinomial(), seed = 314, cores= 4)Expresión del modelo estimado
\[\widehat{\mu^p} = -0.17 +1.53 t -0.45 t^2 -1.30 t^3 -0.04 MES_{February} + 0.06 MES_{March}\\ + 0.02 MES_{April} +0.10MES_{May} -0.05MES_{June} + 0.01MES_{July} +0.15MES_{August} + 0.08MES_{September}\\ + 0.03MES_{October} +0.09MES_{November} -0.04MES_{December} + b^p_{0}\]
\[\widehat{\mu^q} = 3.25 -1.81 t -1.81 t^2 -1.04 t^3 +0.07 MES_{February} - 0.00 MES_{March}\\ + 0.08 MES_{April} +0.08 MES_{May} +0.02MES_{June} + 0.05 MES_{July} +0.06MES_{August} + 0.04 MES_{September}\\ + 0.06 MES_{October} + 0.06 MES_{November} -0.04MES_{December} + b^q_{0}\]
\[\widehat{\mu^r} = -1.26 -3.82 t +3.91 t^2 +0.07 t^3 -0.04 MES_{February} -0.00 MES_{March}\\ + 0.13 MES_{April} + 0.11 MES_{May} +0.14 MES_{June} + 0.16 MES_{July} +0.16MES_{August} + 0.12MES_{September}\\ + 0.17MES_{October} +0.17MES_{November} +0.01MES_{December} + b^r_{0}\]
\[\widehat{\mu^s} = 0.03 +2.43 t -1.67 t^2 + 0.42 t^3 -0.00 MES_{February} - 0.08 MES_{March}\\ -0.05 MES_{April} + 0.06 MES_{May} - 0.00MES_{June} + 0.08MES_{July} +0.05MES_{August} + 0.04MES_{September}\\ + 0.02 MES_{October} - 0.02MES_{November} -0.09MES_{December} + b^p_{0}\]
Con estimación de efectos aleatorios:
ranef(fit1)## $COMUNA
## , , muCaidaOcupante_Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## Aranjuez -0.17633503 0.1103515 -0.39641088 0.03923586
## Belén 0.11234978 0.1100846 -0.10408913 0.32507215
## Buenos Aires 0.02353400 0.1118086 -0.20053442 0.24226193
## Castilla 0.46421326 0.1078272 0.24385820 0.67452564
## Doce de Octubre 0.11530043 0.1095012 -0.10731818 0.32724246
## El Poblado 0.40995730 0.1142855 0.18523787 0.63664518
## Guayabal 0.17442081 0.1109708 -0.04346020 0.39200274
## La América 0.15818200 0.1153638 -0.06896101 0.38222401
## La Candelaria -0.33996102 0.1062530 -0.55341313 -0.13627760
## Laureles Estadio 0.23276200 0.1089999 0.01540726 0.45050019
## Manrique -0.38530005 0.1109014 -0.60915505 -0.16963704
## Popular -0.62892691 0.1161356 -0.85824658 -0.40264120
## Robledo 0.70912014 0.1086911 0.49335820 0.91791067
## San Javier -0.07783837 0.1160330 -0.31644390 0.14968296
## Santa Cruz -0.62164818 0.1199816 -0.86470707 -0.38860883
## Villa Hermosa -0.18321080 0.1125102 -0.41102886 0.03027084
##
## , , muChoque_Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## Aranjuez -0.17964631 0.2269217 -0.623603546 0.28826907
## Belén 0.56543043 0.2275043 0.129524181 1.03688482
## Buenos Aires 0.07276403 0.2288188 -0.365839712 0.54172839
## Castilla 0.49973461 0.2263414 0.063099844 0.95623166
## Doce de Octubre -0.92548920 0.2279626 -1.369956621 -0.46044689
## El Poblado 1.48718262 0.2275926 1.051009034 1.94917304
## Guayabal 0.77844484 0.2279301 0.339581113 1.23864579
## La América 0.45123315 0.2283199 0.008402249 0.92122259
## La Candelaria 0.35785731 0.2264577 -0.082027794 0.81802482
## Laureles Estadio 0.82765424 0.2272842 0.385293218 1.28682341
## Manrique -0.84386603 0.2273367 -1.284379943 -0.37970620
## Popular -1.36875564 0.2291793 -1.807781777 -0.90459627
## Robledo 0.24454496 0.2273629 -0.186766797 0.71259579
## San Javier -0.59852315 0.2283576 -1.037116049 -0.13374706
## Santa Cruz -0.92381070 0.2306542 -1.374102641 -0.44599501
## Villa Hermosa -0.48853926 0.2278985 -0.932435497 -0.02829169
##
## , , muVolcamiento_Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## Aranjuez -0.188126739 0.1375472 -0.47319240 0.079181001
## Belén 0.373479626 0.1371022 0.09775886 0.639987413
## Buenos Aires 0.346128366 0.1386338 0.06788652 0.622129113
## Castilla 0.420207804 0.1340944 0.14279893 0.685373484
## Doce de Octubre -0.548467923 0.1429298 -0.83713319 -0.261412830
## El Poblado 0.823750674 0.1381171 0.54453138 1.095997320
## Guayabal 0.525571749 0.1369229 0.25399286 0.800126041
## La América 0.009832633 0.1461457 -0.28865850 0.295376608
## La Candelaria -0.267687045 0.1320386 -0.53401125 -0.006857978
## Laureles Estadio 0.241080256 0.1354838 -0.03913449 0.509678349
## Manrique -0.367388682 0.1389197 -0.64676351 -0.091710126
## Popular -0.761581347 0.1524242 -1.06914341 -0.456985946
## Robledo 0.336005725 0.1351578 0.07039200 0.603987230
## San Javier -0.109303704 0.1492565 -0.40065837 0.187427923
## Santa Cruz -0.680556167 0.1550127 -0.99782053 -0.381939934
## Villa Hermosa -0.106320531 0.1418453 -0.39432627 0.171114082
##
## , , muOtro_Intercept
##
## Estimate Est.Error Q2.5 Q97.5
## Aranjuez -0.1689083 0.1278364 -0.42466836 0.08429566
## Belén 0.2738271 0.1285301 0.02397624 0.53078837
## Buenos Aires 0.1746305 0.1289375 -0.08417407 0.43319601
## Castilla 0.4956057 0.1272579 0.24262191 0.74995827
## Doce de Octubre -0.2053506 0.1301302 -0.46556586 0.05235482
## El Poblado 0.6529525 0.1295005 0.40040169 0.91058603
## Guayabal 0.3708110 0.1300784 0.11532366 0.63046246
## La América 0.1815452 0.1339048 -0.08320888 0.44299547
## La Candelaria -0.2994731 0.1255367 -0.54818321 -0.05001284
## Laureles Estadio 0.3802947 0.1273853 0.12699345 0.63798496
## Manrique -0.4546204 0.1300006 -0.71787878 -0.19357927
## Popular -0.7852137 0.1328777 -1.05284059 -0.52070687
## Robledo 0.5676279 0.1282647 0.31714458 0.82401303
## San Javier -0.1980646 0.1337548 -0.46126894 0.06514869
## Santa Cruz -0.6689358 0.1367521 -0.93397112 -0.39835821
## Villa Hermosa -0.2406480 0.1306308 -0.49657523 0.01608736Y varianzas de efectos aleatorios estimadas \(\widehat{\sigma^2_{b^p_{0}}} = 0.42\) , \(\widehat{\sigma^2_{b^q_{0}}} = 0.87\) \(\widehat{\sigma^2_{b^r_{0}}} = 0.51\) , \(\widehat{\sigma^2_{b^s_{0}}} = 0.48\)
Comentarios generales
Los modelos mixtos son una herramienta que permite ajustar modelos globales que se caractericen a las particularidades de individuos (en este caso Comunas), lo que permite tener modelos versátiles y precisos para las predicciones.
Una expansión de este trabajo puede ir dirigida a ajustar un modelo en una escala temporal distinta, digamos semanal o diariamente. O utilizar los barrios como la variable de agrupamiento.
Enlazar un modelo Poisson para predecir el número total de accidentes y un modelo multinomial para predecir el número de accidentes por cada clase podría proveer de una herramienta útil para tomar planes de acción respecto a la movilidad de la ciudad. Aunque su utilidad inmediata podría verse afectada por el hecho de que la última vez que miré por mi ventana apenas y había vehículos en tránsito, producto de las medidas estatales de aislamiento social.
Referencias
[1] Ivan Ukhov (2020). Bayesian inference of the net promoter score via multilevel regression with poststratification. Recuperado de https://blog.ivanukhov.com/2020/02/03/net-promoter.html
[2] Example with Family Multinomial. Recuperado de https://discourse.mc-stan.org/t/example-with-family-multinomial/8707
[3] Rob J Hyndman (2016). Cross-validation for time series. Recuperado de https://robjhyndman.com/hyndsight/tscv/
[4] Freddy Hernández Barajas y Jorge Leonardo López Martínez (2020). Modelos Mixtos con R. https://fhernanb.github.io/libro_modelos_mixtos/
[5] Paul-Christian Bürkner (2017). brms: An R Package for Bayesian Multilevel Models Using Stan. Journal of Statistical Software, 80(1), 1-28. doi:10.18637/jss.v080.i01
[6] 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.