Create a gist now

Instantly share code, notes, and snippets.

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

What would you like to do?
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[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)
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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment