Instantly share code, notes, and snippets.

# tslumley/README.md Last active Oct 13, 2017

Simple Bayesian model for road death trends

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.

 library(rjags) safety<-read.csv("safety2.csv") model2.string <-" model { mu0 ~ dgamma(5,5/10) delta~dnorm(-tau2/2,2/tau2) e~ dnorm(-sigma2/2,2/sigma2) mu<-mu0*exp(delta) 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[],d=safety[])) 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[][,17]>b[][,13]) summary(b) table(apply(b[][,1:17],1,which.min)) round(table(apply(b[][,1:17],1,which.min))/1e4,2)
Year Cars occupant deaths Light passenger travel (billion km) Deaths / billion km Relative to 2001 Relative to 2013 2001 358 28.89272819 12.39066099 0 1.215470745 2002 316 29.90935491 10.56525629 -0.147321011 0.889085355 2003 368 30.89352281 11.91188206 -0.038640306 1.129864277 2004 355 31.69724243 11.19971243 -0.096116628 1.002527166 2005 326 31.88484245 10.22429389 -0.174838703 0.828120714 2006 299 31.65550207 9.445435403 -0.237697213 0.688859523 2007 321 32.00855355 10.02856938 -0.190634835 0.793124846 2008 267 31.30079459 8.530134891 -0.311567406 0.525202273 2009 295 31.37139729 9.403470216 -0.241084053 0.681356078 2010 278 31.31393591 8.877836398 -0.283505827 0.587371879 2011 211 30.83284008 6.84335272 -0.447700754 0.223602822 2012 217 30.80635266 7.044001684 -0.43150719 0.259479189 2013 174 31.11148872 5.592789261 -0.548628659 0 2014 197 31.73623193 6.207416194 -0.499024612 0.109896316 2015 232 32.70117862 7.094545512 -0.427428003 0.26851651 2016 240 34.01011008 7.056725176 -0.430480328 0.261754171 2017 206 26.11350115 7.888639627 -0.363339887 0.410501855