Simple Bayesian model for road death trends

README.md

This is a model for NZ road deaths per billion km travelled, using data complied by Sam Warburton (@economissive on Twitter)

The model has an underlying smoothed trend which is a random walk with Gaussian increments, and a year-specific deviation that is also Gaussian. That is, it's a state space model with a random walk in the latent variable and a measurement model that's a Poisson-logNormal mixture to get overdispersed counts.

The variances of the random walk steps and the year-specific deviations have diffuse hyperpriors and the means are set up so that the mean of the exponential of the increment is 1 -- ie, so mu is genuinely the trend in means.

roads.R

library(rjags)
safety<-read.csv("safety2.csv")

model2.string <-"
  model {
  mu0 ~ dgamma(5,5/10)
  delta[1]~dnorm(-tau2/2,2/tau2)
  e[1]~ dnorm(-sigma2/2,2/sigma2)
  mu[1]<-mu0*exp(delta[1])

  for (i in 2:17){
  	delta[i]~dnorm(-tau2/2, 2/tau2)
     e[i]~ dnorm(-sigma2/2, 2/sigma2)
     mu[i]<-mu[i-1]*exp(delta[i])
     x[i] ~ dpois(mu[i]*d[i]*exp(e[i]))
  }
  sigma2~dgamma(1,1)
  tau2~dgamma(1,1)
}
"
model2.spec<-textConnection(model2.string)
jj<-jags.model(model2.spec,data=list(x=safety[[2]],d=safety[[3]]))
update(jj, 1000)
b<-coda.samples(jj,#
             c('mu','sigma2','tau2'),#
             10000)

matplot(2001:2017,summary(b)$quantiles[1:17,],type="l",lty=c(3,2,1,2,3),col=1,ylab="Deaths/billion km",xlab="Year")
qq<-summary(b)$quantiles[1:17,]
polygon(c(2001:2017,2017:2001),c(qq[,2],rev(qq[,4])),col=adjustcolor("blue",alpha.f=0.3),border=NA)
polygon(c(2001:2017,2017:2001),c(qq[,1],rev(qq[,5])),col=adjustcolor("blue",alpha.f=0.1),border=NA)
mean(b[[1]][,17]>b[[1]][,13])
summary(b)
table(apply(b[[1]][,1:17],1,which.min))
round(table(apply(b[[1]][,1:17],1,which.min))/1e4,2)

safety2.csv