###################################################
### chunk number 1: yellowfever
###################################################
x = read.table("yellowfever1.dat",header=TRUE)
a = (x$age1+x$age2)/2  ## average age of category
plot(prot~age1,data=x,type="s",  ## "stair-step" plot
     xlim=c(0,70),ylim=c(0,100),
     main="Yellow fever protection in Amazonia",
     sub="Muench 1934",
     xlab="age",ylab="protection")
points(prot~a,data=x,pch=16)     ## midpoints
unprot = (100-x$prot)/100        ## unprotected
fit1 = lm(log(unprot)~a-1,data=x,subset=prot<100) ## fit
curve(100*(1-exp(coef(fit1)*x)),add=TRUE)  ## add the curve


###################################################
### chunk number 2: getndata
###################################################
ndata = read.csv("niamey_weekly.csv")


###################################################
### chunk number 3: plotndata
###################################################
plot(tot_cases~absweek,data=ndata,log="y",lwd=2,type="l")
matlines(ndata$absweek,ndata[c("cases_1","cases_2","cases_3")],
         col=2:4)
legend("topleft",
       c("tot","cases_1","cases_2","cases_3"),
       col=1:4,lty=c(1,1:3),lwd=c(2,rep(1,3)))


###################################################
### chunk number 4: fit1
###################################################
fit1 = lm(log(tot_cases)~absweek,data=ndata,subset=absweek<=17)
lmc1 = coef(fit1)
curve(exp(lmc1[1]+lmc1[2]*x),add=TRUE) 
abline(v=17,lty=2) ## show the cutoff


###################################################
### chunk number 5: 
###################################################
casedata = ndata[c("tot_cases","cases_1","cases_2","cases_3")]
cumcases = apply(casedata,2,cumsum)  ## cumulative sum of each column
finalsize = cumcases[nrow(cumcases),] ## last row
cumpropcases = scale(cumcases,center=FALSE,scale=finalsize) 
logitcases = qlogis(cumpropcases)


###################################################
### chunk number 6: 
###################################################
matplot(ndata$absweek,cumcases)
abline(h=finalsize,col=1:4)


###################################################
### chunk number 7: 
###################################################
matplot(ndata$absweek,logitcases,
        xlab="Time",ylab="logit(cum prop cases)")
logitcase2 = data.frame(absweek=ndata$absweek,logitcases)
(logitfit1 = lm(tot_cases~absweek,data=logitcase2,subset=absweek<=27))
abline(logitfit1)


###################################################
### chunk number 8: 
###################################################
bombay = read.csv("bombayplague.csv")
plot(mort~week,data=bombay)
curve(890/cosh(0.2*(x-1)-3.4)^2,add=TRUE) ## from Kermack and McKendrick