# Bootstrap estimates for Roberts et al.

library(MASS)
library(boot)

dat = read.csv("iraq.csv",as.is=T)
dat$cmr.pre = 12000*dat$d.pre/dat$pm.pre    # pre-invasion CMR
dat$cmr.post = 12000*dat$d.post/dat$pm.post # post-invasion CMR
dat$cex.size = sqrt(dat$pm.post/max(dat$pm.post)) # cex scaling for labeling
attach(dat)

win.graph()
plot(cmr.pre,cmr.post,ylim=c(-1,30),xlim=c(-1,30),
   xlab="Pre-invasion death rate",ylab="Post-invasion death rate",
   las=1,main="Pre- and post-invasion death rates by cluster",
   type="n",panel.last=abline(0,1,lty=3))
text(cmr.pre,cmr.post,cluster,cex=cex.size)
mtext("Excluding Falluja",line=.5,cex=.8)
text(19.2,20.8,"Higher post-invasion",srt=44,cex=.8,col="cornflowerblue")
text(20,19.25,"Lower post-invasion",srt=44,cex=.8,col="cornflowerblue")
rug(jitter(cmr.pre),side=3,ticksize=.02)
rug(jitter(cmr.post),side=4,ticksize=.02)
points(c(5,31),c(31,7.9),pch=20,col="red",cex=.7)

detach(dat)

# bootstrap estimates
# we'll do it once each with and w/o Falluja
# 17.8 months, 24.4 million total population
# The estimate w/o Falluja is only for 32/33ths of the country

set.seed(1) # make this replicable
dat2 = dat[-21,]  # same thing as dat, but excluding falluja

excessdeaths.nof = function(f,w) 17.8*24400*(32/33)*((sum(f$d.post*w)/sum(f$pm.post*w))-(sum(f$d.pre*w)/sum(f$pm.pre*w)))
excessdeaths.wf = function(f,w) 17.8*24400*((sum(f$d.post*w)/sum(f$pm.post*w))-(sum(f$d.pre*w)/sum(f$pm.pre*w)))

boot.wf = boot(dat,excessdeaths.wf,R=5000,stype="w")
boot.nof = boot(dat2,excessdeaths.nof,R=5000,stype="w")
boot.ci.wf = boot.ci(boot.wf)
boot.ci.nof = boot.ci(boot.nof)

win.graph()
par(mfrow=c(2,1))

truehist(boot.nof$t,h=10,xlab="thousands",ylab="",yaxt="n",
   main="Bootstrap estimates of excess mortality excluding Falluja")
lines(density(boot.nof$t))
abline(v=mean(boot.nof$t),lty=3,col="red")
abline(v=boot.ci.nof$bca[c(4,5)],lty=3,col="red")

truehist(boot.wf$t,h=20,xlab="thousands",ylab="",yaxt="n",
   main="Bootstrap estimates of excess mortality including Falluja",xlim=c(0,1050))
lines(density(boot.wf$t))
abline(v=mean(boot.wf$t),lty=3,col="red")
abline(v=boot.ci.wf$bca[c(4,5)],lty=3,col="red")