################### # WinBUGS code for 1918 flu epidemic model # See Carlin and Louis (2008, Sec 8.3) ################### # S - Susceptibles # E - Asymptomatic infectives # I - Symptomatic infectives # R - Removed = Recovered + Deceased # D - Deceased model { S[1] <- popsz - E1 p <- theta * epsilon E[1] <- E1 + NI[2] R[1] <- R.1 I[1] <- I1 NE[1] <- NE1 NR[1] <- NR1 # day 1, i=2 EI[1] <- round(epsilon*E[1]) r[2] <- exp(beta0+ beta1*(1-centerday)+ beta2*(1-centerday)*(1-centerday)+ beta3*(1-centerday)*(1-centerday)*(1-centerday)) lambda[2]<-r[2]*S[1]*(EI[1]+rho*delta*I[1]) /(S[1]+E[1]+delta*I[1]+gamma*R[1]) NE[2] ~ dpois(lambda[2]) NI[2] ~ dbin(theta,E[1]) NI.mean[2] <- E[1]*theta NI.sample[2] ~ dbin(theta,E[1]) NR[2] <- NR2 ND.mean[2] <- round(NR[2]*(1-gamma)) S[2] <- S[1] - NE[2] E[2] <- E[1] + NE[2] - NI[2] I[2] <- I[1] + NI[2] - NR[2] R[2] <- R[1] + NR[2] RN[2] <- epsilon*[2]/theta+rho*delta*r[2]/phi # day 2 to day 91 for (i in 3:92) { EI[i-1] <- round(epsilon*E[i-1]) r[i] <- exp(beta0+ beta1*(i-1-centerday)+ beta2*((i-1-centerday)*(i-1-centerday))+ beta3*((i-1-centerday)*(i-1-centerday)*(i-1-centerday))) lambda[i]<-r[i]*S[i-1]*(EI[i-1]+rho*I[i-1]) /(S[i-1]+E[i-1]+delta*I[i-1]+gamma*R[i-1]) NE[i] ~ dpois(lambda[i]) NI[i] ~ dbin(theta, E[i-1]) NI.sample[i] ~ dbin(theta,E[i-1]) NI.mean[i] <- E[i-1]*theta ND[i] ~ dbin(eta,I[i-1]) ND.sample[i] ~ dbin(eta, I[i-1]) ND.mean[i] <- I[i-1] * eta NR[i] <- round(ND[i]/(1-gamma)) S[i] <- S[i-1] - NE[i] E[i] <- E[i-1] + NE[i] - NI[i] I[i] <- I[i-1] + NI[i] - NR[i] R[i] <- R[i-1] + NR[i] RN[i] <- epsilon*r[i]/theta*S[i]/(S[i]+E[i]+delta*I[i]+gamma*R[i]) +r[i]/phi*rho*delta*S[i]/(S[i]+E[i]+delta*I[i]+gamma*R[i]) # The variables ``NI.sample" and ``ND.sample" are the simulated # incidence and mortality. } beta0 ~ dnorm(0,.0001) beta1 ~ dnorm(0,.0001) beta2 ~ dnorm(0,.0001) beta3 ~ dnorm(0,.0001) eta ~ dbeta(1,1) gamma ~ dbeta(1,1) phi <- eta/(1-gamma) dura <- 1/phi } # initial values list(eta=.0026, NE=c(NA,1,1,1,1,1,1,1,1,1,1,1,1,2,3,5,10, 8,15,9,18,27,18,18,22,29,30,40,50,65,79,98,106,124,164,202,267, 343,294,328,344,386,400,400,400,400,400,400,400,400,400,400, 300,300,300,200,200,200,200,200,200,200,100,100,100,100,100, 100,100,100,100,100,100,100,100,100,100,100,100,100,100,100, 100,100,100,100,100,100,100,100,100,100,100,100,100,100,100,50,50 ),gamma=.98, beta0=-0.16, beta1=-0.016, beta2=0, beta3=0) # data list(NI=c(NA,2,0,1,3,1,2,0,8,1,7,18,3,5,12,9,2,17,4,21,4,11,37,14,13,37,17, 43,21,18,39,40,88,70,41,87,62,107,481,225,294,255, 337,260,360,327,228,443,186,264,164,205,628,182,131,172,100, 161,112,90,80,75,135,57,50,83,55,73,34,60,41,29,22,14,24,13,14, 10,19,5,5,9,10,17,10,10,12,7,10,3,2,4,5,10,5,6,3,6,1),popsz = 33776,ND=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0, 0,0,2,0,0,0,0,0,3,0,3,5,5,1,10,3,13,8,4,5,10,5,7,5,3,7,3,2,5,1, 2,2,5,2,4,0,0,2,1,1,1,0,1,2,1,0,1,1,0,1,0,1,1,2,0,0,0,1,0,0, 1,0,0,0,0,0,0),R.1 =0, NR2 =0, I1=0, NE1=0, NR1=0, E1=2, theta=0.5 epsilon=0.5, rho=0.41, delta=0.001)