{{ message }}

Instantly share code, notes, and snippets.

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.

This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
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