trafficdeaths.Rmd

---
title: "Bayesian data analysis - traffic deaths in Finland"
author: "Aki Vehtari"
date: '`r format(Sys.Date())`'
output:
  html_document: default
  html_notebook: default
---

License: CC-BY

This notebook demonstrates time series analysis for traffic deaths per year in Finland. Currently when the the number of traffic deaths during previous year are reported, the press release claims that the the traffic safety in Finland has improved or worsened depending whether the number is smaller or larger than the year before. Time series analysis can be used to separate random fluctuation from the slowly changing traffic safety. 

Load some libraries:
```{r, comment=NA}
library(ggplot2)
library(tidyr)
library(gridExtra)
library(rstanarm)
library(rstan)
library(bayesplot)
library(loo)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
source("stan_utility.R")
```

Read the data (there would data for earlier years, too, but this is sufficient for the demonstration)
```{r}
# file preview shows a header row
deaths <- read.csv("trafficdeaths.csv", header = TRUE)
head(deaths)
```

First plot just the data.
```{r}
ggplot() +
  geom_point(aes(year, deaths), data = deaths, size = 1) +
  labs(y = 'Traffic deaths', x= "Year") +
  guides(linetype = F) +
  theme_bw()
```

The number of deaths is count data, so we use Poisson observation model. First we use log-linear model for the Poisson intensity, which corresponds to assuming constant proportional change in the rate.
```{r, comment=NA}
fit_lin <- stan_glm(deaths ~ year, deaths, family="poisson", refresh=1000, iter=1000, chains=4, seed=583829)
```
```{r, comment=NA}
summary(fit_lin)
```

n_eff's and Rhat's are ok. Let's look at the posterior predictive distribution (median and 5% and 95% intervals).
```{r}
x_predict=seq(1993,2023)
N_predict=length(x_predict)
y_predict <- posterior_predict(fit_lin, newdata=data.frame(year=x_predict))
mu <- apply(t(y_predict), 1, quantile, c(0.05, 0.5, 0.95)) %>%
  t() %>% data.frame(x = x_predict, .) %>% gather(pct, y, -x)
pfit <- ggplot() +
  geom_point(aes(year, deaths), data = deaths, size = 1) +
  geom_line(aes(x, y, linetype = pct), data = mu, color = 'red') +
  scale_linetype_manual(values = c(2,1,2)) +
  labs(x = 'Year', y = 'Traffic deaths') +
  guides(linetype = F) +
  theme_bw()
(pfit)
```

Next we test non-linear spline model with stan_gamm4
```{r}
fit_gam <- stan_gamm4(deaths ~ s(year), data=deaths, family="poisson", refresh=1000, iter=1000, chain=4, seed=583829)
summary(fit_gam)
```

n_eff is clearly smaller than for the linear model, but Rhat's are ok.

Let's look at the posterior predictive distribution.
```{r}
x_predict=seq(1993,2023)
N_predict=length(x_predict)
y_predict <- posterior_predict(fit_gam, newdata=data.frame(year=x_predict))
mu <- apply(t(y_predict), 1, quantile, c(0.05, 0.5, 0.95)) %>%
  t() %>% data.frame(x = x_predict, .) %>% gather(pct, y, -x)
pfit <- ggplot() +
  geom_point(aes(year, deaths), data = deaths, size = 1) +
  geom_line(aes(x, y, linetype = pct), data = mu, color = 'red') +
  scale_linetype_manual(values = c(2,1,2)) +
  labs(x = 'Year', y = 'Traffic deaths') +
  guides(linetype = F) +
  theme_bw()
(pfit)
```

The predictive median is clearly nonlinear. The predictive mean for future years stays at the same level as the most recent observations, but uncertainty increases quickly.

Finally we test Gaussian process centered on linear model. This is not yet available in rstanarm, and has been written directly in Stan language:
```{r, comment=NA}
writeLines(readLines("poisson_gp.stan"))
```

```{r}
N<-nrow(deaths)
Ey<-mean(deaths$deaths)
d_data <- list(N=N, x=deaths$year, y=deaths$deaths, Ey=Ey, N_predict=N_predict, x_predict=x_predict, alpha0=2, beta0=4)
fit_gp <- stan(file='poisson_gp.stan', data=d_data, refresh=1000, iter=1000,
                     chains=4, seed=583829, init=0, control=list(adapt_delta=0.99))
```
```{r}
check_treedepth(fit_gp)
check_energy(fit_gp)
check_div(fit_gp)
```

```{r}
gp_params <- rstan::extract(fit_gp)
mu <- apply(t(gp_params$y_predict), 1, quantile, c(0.05, 0.5, 0.95)) %>%
  t() %>% data.frame(x = x_predict, .) %>% gather(pct, y, -x)
pfit <- ggplot() +
   geom_point(aes(year, deaths), data = deaths, size = 1) +
  geom_line(aes(x, y, linetype = pct), data = mu, color = 'red') +
  scale_linetype_manual(values = c(2,1,2)) +
  labs(x = 'Year', y = 'Traffic deaths') +
  guides(linetype = F) +
  theme_bw()
(pfit)
```

Finally PSIS-LOO estimates. 

```{r}
(loo_lin<-loo(fit_lin))
(loo_gam<-loo(fit_gam))
compare(loo_lin,loo_gam)
(loo_gp<-loo(gp_params$log_lik))
compare(loo_lin,loo_gp)
```

There are no practical differences in predictive performance, which is partially due to small number of observations. Based on the posterior predictive distributions there are clear differences in the future predictions.

<br />

### Appendix: Session information

```{r}
sessionInfo()
```

<br />


### Appendix: Licenses

* Code &copy; 2017, Aki Vehtari, licensed under BSD-3.
* Text &copy; 2017, Aki Vehtari, licensed under CC-BY-NC 4.0.