Número de acidentes por dia nos últimos 10 anos

Veja a seguir o número de acidentes, por id, diariamente e duas curvas de médias móveis com k = 30 e com k = 30*6 :

g1 <- 
  base %>%
  map_df(~.x%>% select(id, data_inversa, causa_acidente, tipo_acidente) ) %>% 
  dplyr::group_by(data_inversa)%>% 
  dplyr::summarise(n=n())%>%
  mutate(rMM=rollmean(n,30, na.pad=TRUE, align="left")) %>%
  mutate(rMY=rollmean(n,30*6, na.pad=TRUE, align="left")) %>%
  tidyr::gather(medida,valor,n:rMY) %>% 
  ggplot(aes(x=data_inversa, y=valor,col=medida)) +
  geom_line() +
  scale_colour_manual(name="Legenda:",values=c("black","red","blue"), 
                      labels = c('Acidentes','Médias móveis k=30', 'Médias móveis k=30*6')) +
  theme_bw() +
  labs(y= 'Número de acidentes',
       x = 'Data',
       subtitle = "Por dia", 
       title= 'Número de acidentes na Ponte Rio-Niterói', 
       caption = "Fonte dos dados:\nhttps://www.prf.gov.br/portal/dados-abertos/acidentes\nhttps://gomesfellipe.github.io/",col="") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(breaks=seq(0,15,1))

g1 +
  scale_x_date(date_breaks = '6 months',limits = as.Date(c("2008-02-01","2018-12-01")),
               expand = c(0.01,0.01)) + 
  theme(legend.position = c(0.85,0.7), 
        legend.background = element_rect(fill=alpha('lightgrey', 0.2)))

Este gráfico esta um pouco carregado com a informação levando em conta o número de acidentes por dia mas mesmo assim já é possível notar algumas informações como:

• Existem diversos dias que não houveram acidentes registrados pela PRF (diversos pontos tocam o eixo x)
• No ano de 2015 houve uma grotesca diminuição no registros do número de acidentes (o que será que aconteceu?)

A primeira pergunta parece ter sido respondida: O número de acidentes não parece ter sido diário nos últimos 10 anos segundo esses dados (e supondo que não houve subnotificação).

Já a hipótese levantada com a segunda pergunta parece ser verdadeira pois o comportamento dos dados sugerem ter havido a diminuição no número de acidentes.

Veja a seguir outra rápida análise exploratória para checar como se distribuem os número de acidentes ano longo dos anos e dos meses com alguns box-plots:

bp1 <- 
  base %>%
  map_df(~.x%>% select(id, data_inversa, causa_acidente, tipo_acidente) ) %>% 
  dplyr::group_by(data_inversa)%>% 
  dplyr::summarise(n=n()) %>%
  tidyr::gather(medida,valor,n) %>% 
  ggplot(aes(x=factor(year(data_inversa)), y=valor,col=medida)) +
  geom_jitter(col="darkgrey",alpha=0.3)+
  geom_boxplot(col="black",alpha=0.3)+
  theme_bw() +
  labs(y= 'Número de acidentes',x="") + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1)) +
  scale_y_continuous(breaks=seq(0,15,1))+ 
  geom_vline(xintercept = 8.5, linetype="dotted", 
             color = "red", size=1.5) +
  annotate("label",y = 7, x = 8.5,
           label = "Comportamentos \n     diferentes",
           size = 4, colour = "red",hjust=0.1)

bp2 <- 
  base %>%
  map_df(~.x%>% select(id, data_inversa, causa_acidente, tipo_acidente)) %>% 
  filter(year(data_inversa)<2016) %>% 
  dplyr::group_by(data_inversa) %>% 
  dplyr::summarise(n=n()) %>%
  tidyr::gather(medida,valor,n) %>% 
  ggplot(aes(x=factor(months(data_inversa),
                      levels = months(base[[1]]$data_inversa) %>% unique),
             y=valor,col=medida)) +
  geom_jitter(col="darkgrey",alpha=0.3)+
  geom_boxplot(color="black",alpha=0.3)+
  theme_bw() +
  labs(y= 'Número de acidentes',
       title = 'Dados de 2008 até 2015',x="", col="") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(breaks=seq(0,15,1))+ 
  theme(legend.position = c(0.85,0.7))

bp3 <- 
  base %>%
  map_df(~.x%>% select(id, data_inversa, causa_acidente, tipo_acidente) ) %>% 
  filter(year(data_inversa)>=2016) %>% 
  dplyr::group_by(data_inversa) %>% 
  dplyr::summarise(n=n()) %>%
  tidyr::gather(medida,valor,n) %>% 
  ggplot(aes(x=factor(months(data_inversa),
                      levels = months(base[[1]]$data_inversa) %>% unique),
             y=valor,col=medida)) +
  geom_jitter(col="darkgrey",alpha=0.3)+
  geom_boxplot(color="black",alpha=0.3)+
  theme_bw() +
  labs(y= '',
       title = 'Dados de 2016 até 2018',x="",col="") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(breaks=seq(0,15,1))+ 
  theme(legend.position = c(0.85,0.7))

bp2 <- arrangeGrob(bp2,bp3,ncol=2)

grid.arrange(bp1,bp2,
             top = "Boxplots do número de acidentes na Ponte Rio-Niterói",
             bottom = textGrob("Fonte dos dados:\nhttps://www.prf.gov.br/portal/dados-abertos/acidentes\nhttps://gomesfellipe.github.io/",
                               gp = gpar(fontface = 3, fontsize = 9),
                               hjust = 1,x = 1 )) 

É possível notar que antes de 2016 todos aos boxplots apresentaram medianas parecidas enquanto que após o ano de 2016 parece haver mais acidentes apenas nos primeiros meses do ano.

Para responder a última pergunta vamos analisar com mais detalhes a série para entender melhor seu comportamento.

g1 +
  facet_wrap( ~ year(data_inversa), scales = "free_x") +
  ggtitle("Número de acidentes na Ponte Rio-Niterói de acordo com o ano, por dia") +
  scale_x_date(date_breaks = '1 month') +
  theme(legend.position = c(0.85, 0.1))

Mas o que aconteceu em 2015?

Já que em 2015 houve este comportamento tão brusco resolvi procurar no Google por notícias de acidentes na Ponte Rio-Niterói neste período e em uma breve pesquisa me deparei com a seguinte notícia:

Essa afirmação na manchete, que houve “redução de quase 60% após o videomonitoramento” é intrigante e seria interessante conseguir achar algum resultado concordante.

Na notícia " Monitoramento com câmeras ajuda a reduzir acidentes na Ponte Rio-Niterói " um inspetor da PRF faz a seguinte declaração:

“(…) Ajuda a amplificar o nosso poder de visualização das infrações que são cometidas. Concequentemente reduz o numero de acidentes por conta da presença quase que em todo trecho, ao mesmo tempo, da PRF”

Número de acidentes por mês nos últimos 10 anos agrupados

Como o gráfico ficou carregado ao levar em consideração o número de acidentes ao longo dos dias, a mesma informação será agrupada por mês.

Vejamos visualmente os dados fornecidos pela PRF e o dia em que a notícia foi atualizada:

g2 <- 
  base %>%
  map_df(~.x
         %>% select(id, data_inversa, causa_acidente, tipo_acidente) %>%
           mutate_if(is.character,~ dmy(.x) )) %>%
  group_by(month=floor_date(data_inversa, "1 month")) %>% 
  summarise(n = n()) %>%  
  ggplot(aes(x=month, y=n)) +
  geom_line() +
  geom_point() +
  theme_bw() + 
  labs(y= 'Número de acidentes',
       x = 'Data',
       title= 'Número de acidentes na Ponte Rio-Niterói', 
       subtitle = "Por mês", 
       caption = "Fonte dos dados: \nhttp://glo.bo/210j8gF \nhttps://www.prf.gov.br/portal/dados-abertos/acidentes\nhttps://gomesfellipe.github.io/") +
  scale_y_continuous(breaks=seq(0,150,25))

g2 +
  geom_segment(x = as.Date("2015-11-13"), xend = as.Date("2015-11-13"),
               y=0,yend=100, linetype="dashed", col="red") +
  geom_text(aes(y = 100, x = as.Date("2015-11-13"),
                label = "G1 em 17/11/2015 :\n\"Monitoramento com câmeras\n ajuda a reduzir acidentes\n  na Ponte Rio-Niterói\"") ,hjust=-0.1) +
  scale_x_date(date_breaks = '6 month') + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

De fato, visualmente é possível notar uma aparente diminuição no ano do qual a notícia foi divulgada mas será que é possível identificar algum tipo de tendência a partir destes dados?

g2 +
  # geom_smooth(method = "lm", se = F, aes(col = "blue")) +
  facet_wrap( ~ year(month), scales = "free_x") +
  scale_colour_manual(name=" ", values=c("blue"),labels=c("Regressão Linear")) + 
  theme(legend.position = c(0.88,0.20))

Ajuste do modelo descritivo

Para entender melhor essa tendência apresentada em 2015 será ajustado modelo de regressão linear simples conforme metodologia apresentada em @tend, que utilizam o tempo como variável explicativa e entender qual foi a variação do número de acidentes médios por mês.

O método utilizado para avaliar a tendência consiste em um ajuste de tendência em função do tempo (\(t\)) com um modelo de regressão linear de maneira que:

\[ Y_{i} = \beta_0 + \beta_1 t_i + \epsilon_{i} \]

em que:

• \(Y_i\): Tendência para o número de acidentes médios por mês no instante \(i\)
• \(t_i\): i-ésimo instante de tempo
• \(\beta_0\): Intercepto do modelo (valor de \(Y_i\) quando \(t_i=0\))
• \(\beta_1\): Inclinação da reta ajustada (efeito da variável \(t_i\))
• \(\epsilon_i \sim N(0,\sigma^2)\) onde \(\sigma^2\) é constante

Este modelo de regressão linear mede a variação da tendência para diferentes valores de intervalo de tempo na faixa de interesse. O ajuste é interpretado de acordo com a inclinação da reta formada pelos diferentes valores de referência em relação a respectiva tendência pode não ser o melhor método para a previsão mas pode ajudar a entender o comportamento da série como análise descritiva.

to_lm <- 
  base %>%
  map_df(~.x%>% select(id, data_inversa, causa_acidente, tipo_acidente)%>%
           mutate_if(is.character,~ dmy(.x) )) %>% 
  group_by(data_inversa=floor_date(data_inversa, "1 month")) %>% 
  summarise(n = n()) %>% 
  filter(data_inversa >= "2015-01-01" & data_inversa < "2016-01-01")

y   <- to_lm %$% n %>%  ts(start = c(2015, 01, 01))

g2 +
  geom_smooth(method='lm',se=F,aes(col=I("darkred"))) +
  geom_smooth(data=to_lm, aes(x=data_inversa, y=n,col=I("blue")), method = "lm",se=F) +
  scale_colour_manual(name=NULL, values=c("red","blue","darkred"),labels=c(
    "Regressão Linear para 2015",
    "Regressão Linear para todo período",
    "Notícia"
  )) + 
  theme(legend.position = c(0.2,0.18)) +
  geom_segment(x = as.Date("2015-11-13"), xend = as.Date("2015-11-13"),
               y=0, yend=100, linetype="dashed", aes(col=I("red"))) +
  geom_text(aes(y = 100, x = as.Date("2015-11-13"),
                label = "G1 em 17/11/2015 :\n\"Monitoramento com câmeras\n ajuda a reduzir acidentes\n  na Ponte Rio-Niterói\" ") ,hjust=-0.1) +
  scale_x_date(date_breaks = '6 month') + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1), 
        legend.background = element_rect(fill=alpha('lightgrey', 0.2)))

Antes de conferir os resultados do ajuste do modelo, veja o resultado para o teste de sequências (Wald-Wolfowitz) que tem como hipótese alternativa \(H_1:\) existe tendência na série

randtests::runs.test(y)

## 
##  Runs Test
## 
## data:  y
## statistic = -3.0277, runs = 2, n1 = 6, n2 = 6, n = 12, p-value =
## 0.002465
## alternative hypothesis: nonrandomness

Ao nível de significância de \(\alpha :\) 0.05 rejeita-se a hipótese nula de que esta série não apresenta tendência no período selecionado (2015).

Outra hipótese importante para conferir é sobre a normalidade da variável resposta, a seguir será aplicado o teste de Shapiro-wilk que tem como hipótese nula de que os dados possuem distribuição normal:

shapiro.test(y)

## 
##  Shapiro-Wilk normality test
## 
## data:  y
## W = 0.9702, p-value = 0.9129

De acordo com este teste, não existem evidências estatísticas para se rejeitar a hipótese da normalidade dos dados, portanto, para uma análise mais detalhada o modelo ajustado no período de 2015 pode ser obtido da seguinte maneira:

fit <- tslm(y ~ trend, data = y) 

Com o modelo está salvo no objeto fit, sua informação será incrementada no gráfico da regressão de 2015 para facilitar o entendimento dos resultados:

g1 <- 
  base %>%
  map_df(~.x%>% select(id, data_inversa, causa_acidente, tipo_acidente)) %>% 
  group_by(data_inversa=floor_date(data_inversa, "1 month")) %>% 
  dplyr::summarise(n=n()) %>%
  ggplot(aes(x=data_inversa, y=n)) +
  geom_line() +
  theme_bw() +
  labs(y= 'Número de acidentes',
       x = 'Data',
       subtitle = "Por dia", 
       title= 'Número de acidentes na Ponte Rio-Niterói', 
       caption = "Fonte dos dados:\nhttp://glo.bo/210j8gF\nhttps://www.prf.gov.br/portal/dados-abertos/acidentes\nhttps://gomesfellipe.github.io/") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  scale_y_continuous(breaks=seq(0,100,25)) + 
  annotate("text", x=as.Date("2015-08-13"), y = 75,
           label = paste("Modelo teórico: Yi = β^0 + β1 ti + ei:\nEstimador do Nº de acidentes =",round(fit$coefficients[1],2),round(fit$coefficients[2],2),"*  Mês"))

t <- 
  fit %>%
  tidy() %>% 
  bind_cols(r2=c(glance(fit)[[1]] %>% round(4),"-"),
            confint_tidy(fit)) %>% 
  mutate_if(is.numeric,~ round(.x,4)) %>% 
  mutate(term = c("Intercepto", "Tendência")) %>% 
  set_names(c("Coef.", "Estimativa", "Erro Padrão", "Estatística", "Valor P", "R²","Lim. Inf.", "Lim. Sup.")) 

# Funcao auxiliar para colorir a celula
g <- tableGrob(t,rows = NULL)
find_cell <- function(table, row, col, name="core-fg"){
  # Fonte: 
  # https://cran.r-project.org/web/packages/gridExtra/vignettes/tableGrob.html
  l <- table$layout
  which(l$t==row & l$l==col & l$name==name)
}

ind <- find_cell(g, 3, 2, "core-bg")
g$grobs[ind][[1]][["gp"]] <- gpar(fill="darkolivegreen1", col = "darkolivegreen4", lwd=5)

g1 +
  scale_x_date(date_breaks = '1 months',
               limits = as.Date(c("2015-01-01","2015-12-01")),
               expand = c(0.01,0.01)) +
  geom_smooth(data=to_lm, aes(x=data_inversa, y=n), method = "lm",col="red",se=F)+
  annotation_custom(grob = g  , ymax=220,xmin=as.Date("2015-04-01"))

Nota-se que o coeficiente estimado para determinar a variação do número de acidentes ao aumentar a variável explicativa trend em uma unidade se apresenta de maneira significante, isto quer dizer que parece existir relação linear estatisticamente significante entre o mês da observação e a variável resposta.

Além disso é possível notar que para cada aumento de uma unidade da tred (ou seja, no decorrer de 1 mês), o valor esperado no número de acidentes observados naquele mês diminui em aproximadamente 4.7972.

Como -4.7972 x 12 = -57.5664 então o número de acidentes esperados decresce para aproximadamente 57.57% do valor original conforme a variável tred aumentou em 12 unidades (em 1 ano).

Veja o resultado para o teste de Shapiro - wilk para avaliar a hipótese de normalidade dos resíduos do modelo, os resíduos padronizados(quociente entre o resíduo e a estimativa do seu desvio padrão) e os resíduos studentizados, que levam em conta o \(h\) (medida de alavancagem ou " leverage") em sua medida:

bind_cols(
  residuals = fit$residuals,
  rstandard = rstandard(fit),
  rstudent = rstudent(fit)
) %>%
  map_dbl(~shapiro.test(.x)$p.value)  %>%
  tidy() %>% 
  set_names(c("Medida", "Valor p")) %>% 
  kable(escape = F) %>%
  kable_styling("hover",full_width = F)

Não existem evidências para rejeitar a hipótese de que os resíduos do modelo e os resíduos padronizados seguem distribuição normal. Já para os resíduos studentizados rejeita-se a hipótese de normalidade e como o valor de R² não foi muito alto (66.51%), indicando que a reta ajustada explica parcialmente a variação dos dados então pode ser interessante estudar com mais detalhes qual foi a influência do número de acidentes em cada mês no ajuste do modelo calculando as medias de influência para cada observação do modelo.

Uma observação é considerada influente caso sua exclusão no ajuste do modelo de regressão cause uma mudança substancial na análise de regressão e para identificar essas observações serão utilizadas as mesmas medidas utilizadas no post que fiz sobre Pacotes do R para avaliar o ajuste de modelos.

De maneira visual e com medidas de influência do valor observado para cada mês no ajuste do modelo:

t1  <- influence.measures(fit) %$%
  infmat %>% 
  as_tibble() %>% 
  map_df(~.x %>% round(4)) %>% 
  mutate(mes = 1:12) %>% 
  select(mes,everything()) %>% 
  set_names(c("Mês","DF β1","DF β2","   DFFit","    CovRatio"," D.Cook","   h"))

g   <- tableGrob(t1,rows = NULL)

ind <- map_dbl(1:7, ~find_cell(g, 6, .x, "core-bg"))  

walk(ind, function(i) g$grobs[i][[1]][["gp"]] <<-
       gpar(fill = "darkolivegreen1", col = "darkolivegreen1", lwd = 5)
)

# g2  <-
#   ggplot2::autoplot(fit, print=F,suppressMessages=T) +
#   theme_bw()
# 
# gn  <- arrangeGrob(g2[[1]], g2[[2]], g2[[3]], g2[[4]], g2[[5]], g2[[6]])
# 
# grid.arrange(gn, g, ncol=2)

grid.arrange(g)

A tabela apresenta um resumo das medidas de inflencia de cada variável, segundo @gauss:

• DF Beta: Alteração no vetor estimado β ao se retirar o i-ésimo ponto da análise
  
• DF Fit: Alteração provocada no valor ajustado pela retirada da observação i
  
• CovRatio: Expressa o relação de covariância
  
• D.Cook: Medida de afastamento das estimativas ao retirar i e também considera o resíduo estudentizado internamente
  
• h: Elementos da diagonal da matriz H

É possível notar que no mês 5 (Maio) foi o mês que provocaria a maior alteração provocada no valor ajustado após sua retirada. Outra medida importante é a Distância de Cook que indica o quão atípica a observação \(i\) se apresenta no ajuste do modelo combinando os resíduos estudentizados e as medidas de alavancagem possibilitando examinar observações que influenciam bastante as estimativas dos parâmetros.

Em geral, se \(D_i>1\) o ponto é excessivamente influente mas como no geral as observações apresentaram valor muito inferior a 1 então a eliminação de qualquer uma dessas variáveis não vai alterar substancialmente as estimativas dos parâmetros porém nós podemos examinar quaisquer observações cujo \(D_i\) seja muito maior em relação aos demais valores estimados para \(D_i\) com mais cuidado para entender qual foi foi sua influência no ajuste do modelo:

# Como a funcao augment do broom nao parece mais suportar o modelo ajustado por tslm, vamos ajustar na mão o modelo com a tendencia para criar o gráfico:
to_lm <- to_lm %>% mutate(trend = 1:nrow(.))
fit_lm <- lm(y ~ trend, data = to_lm) 


fit_lm %>% 
  broom::augment() %>% 
  mutate(col = if_else(.resid==max(.resid), T, F))%>%
  mutate(trend = factor(months(base[[1]]$data_inversa) %>% unique,
                        levels = months(base[[1]]$data_inversa) %>% unique) ) %>% 
  ggplot(aes(x=trend,y=.resid,col=col,label=trend)) +
  geom_point() +
  geom_label(vjust = -0.25) +
  theme_bw() +
  scale_y_continuous(limits=c(-20,35),breaks = seq(-20,35,5)) +
  scale_color_manual(values=c("black", "red")) +
  geom_hline(yintercept = 0,linetype="dashed",col="red") +
  theme(legend.position="none") +
  labs(y="Resíduos", x="Valores ajustados", title="Resíduos vs Ajuste")

De fato o valor ajustado para o mês de maior é muito diferente do valor real observado, o que implica que houve um comportamento atípico que pode ser confirmado ao voltar no gráfico da série temporal e constatar um comportamento diferente (o número de acidente aumentou em relação ao anterior) e fugindo da tendência geral dos dados.