###################################################
### chunk number 1: input
###################################################
##############################
# input cases
##############################
niamey_cases1<-c(4,1,12,17,37,55,51,331,323,370,145,224,162,26,6)
niamey_cases2<-c(1,2,3,6,9,24,40,55,88,158,173,155,141,36,20,2)
niamey_cases3<-c(2,7,4,14,21,49,73,151,279,245,185,129,49,19,2)


###################################################
### chunk number 2: plot1
###################################################
## basic plot
plot(niamey_cases1,type="l",xlab="biweek",ylab="measles cases")
## add lines, one at a time, different line types
lines(niamey_cases2, lty=2)
lines(niamey_cases3, lty=3)
legend("topleft",
       legend=c("Center 1","Center 2","Center 3"),lty=1:3)
## IF case sequences were the same length, we could use
## matplot(...,type="l")


###################################################
### chunk number 3: likfuns
###################################################
likelihood <- function(S0,beta,I){
    ## cat(S0,beta,"\n") ## debugging statement (commented out)
    n <- length(I)
    S <- floor(S0-cumsum(I[-n]))
    p <- 1-exp(-beta*(I[-n]))	
    L <- dbinom(I[-1],S,p,log=TRUE)	
    Lik<-sum(-L,na.rm=TRUE)
}
## OR -- define vectors
##likelihood <- function(par,vecs){
##     S0<-par[1]
##     beta<-par[2]
##     I<-vecs$I
##     n <- length(I)
##     S<-floor(S0-cumsum(I[-n]))
##     p<-1-exp(-beta*(I[-n]))	
##     L<-dbinom(I[-1],S,p,log=TRUE)	
##     Lik<-sum(-L,na.rm=TRUE)
## }


###################################################
### chunk number 4: calcsurf
###################################################
tot1 <- sum(niamey_cases1)
S0<-floor(seq(tot1+1,3*tot1,length=40))
beta<-seq(1.5/S0[1],5/S0[1],length=40)
lik<-matrix(NA,nrow=40,ncol=40)
for(i in 1:40){
    for(j in 1:40){
        lik[i,j]<-likelihood(S0=S0[i],beta=beta[j],
                             I=niamey_cases1)
    }
}


###################################################
### chunk number 5: calcsurf2
###################################################
library(emdbook)
lik <- apply2d(likelihood,S0,beta,I=niamey_cases1)


###################################################
### chunk number 6: calcsurf3
###################################################
lik <- curve3d(likelihood(x,y,I=niamey_cases1),from=c(tot1+1,1.5/S0[1]),
             to=c(3*tot1,5/S0[1]),sys3d="contour",n=c(40,40),
             levels=200*(2:11))$z


###################################################
### chunk number 7: calcmin
###################################################
mle.ind<-which(lik==min(lik),arr.ind=TRUE)
init<-list(S0=S0[mle.ind[1]],beta=beta[mle.ind[2]])


###################################################
### chunk number 8: plotsurf
###################################################
symbols(rep(S0,40),rep(beta,each=40),
        circles=50*as.vector(lik)/max(lik),
        xlab=~S[0],ylab=~beta,
        inches=FALSE)
## ~ is a shorthand to get R to interpret labels as math expressions
points(S0[mle.ind[1]],beta[mle.ind[2]],col=2,pch=16)
contour(S0,beta,lik,levels=(200*(2:11)),add=TRUE)


###################################################
### chunk number 9: fit1 eval=FALSE
###################################################
## mle.fit<-optim(init,likelihood,
##                method="L-BFGS-B", 
##                lower=c(sum(niamey_cases1),1e-6),upper=c(Inf,0.1),
##                vecs=list(I=niamey_cases1),control=list(ndeps=c(1,1e-5)),
##                hessian=TRUE)


###################################################
### chunk number 10: fit1B
###################################################
library(bbmle)
mle.fit2 <- mle2(start=init,
                 likelihood,
                 method="L-BFGS-B", 
                 lower=c(tot1,1e-6),upper=c(Inf,0.1),
                 data=list(I=niamey_cases1),
                 control=list(ndeps=c(1,1e-5)))


###################################################
### chunk number 11: 
###################################################
newstart <- list(S0=2.153976e+03,beta=9.697629e-04)
mle.fit2 <- mle2(start=newstart,
                 likelihood,
                 method="L-BFGS-B", 
                 lower=c(tot1,1e-6),upper=c(Inf,0.1),
                 data=list(I=niamey_cases1),
                 control=list(ndeps=c(1,1e-5)))


###################################################
### chunk number 12: coefs
###################################################
## mle.fit$par
## OR
coef(mle.fit2)


###################################################
### chunk number 13: addlikcont
###################################################
contour(S0,beta,lik,levels=-logLik(mle.fit2)+qchisq(0.95,2)/2, 
        col="red", add=TRUE)


###################################################
### chunk number 14: vcov
###################################################
(covmat=vcov(mle.fit2))


###################################################
### chunk number 15: cormat
###################################################
(cor.mat<-cov2cor(covmat))


###################################################
### chunk number 16: plot2
###################################################
symbols(rep(S0,40),rep(beta,each=40),
        circles=50*as.vector(lik)/max(lik),
        xlab=~S[0],ylab=~beta,
        inches=FALSE)
## ~ is a shorthand to get R to interpret labels as math expressions
points(S0[mle.ind[1]],beta[mle.ind[2]],col=2,pch=16)
contour(S0,beta,lik,levels=(200*(2:11)),add=TRUE)
contour(S0,beta,lik,levels=-logLik(mle.fit2)+qchisq(0.95,2)/2, 
        col="red", add=TRUE)
library(ellipse)
lines(ellipse(vcov(mle.fit2),centre=coef(mle.fit2)),col=4)


###################################################
### chunk number 17: plot3
###################################################
c2 = curve3d(likelihood(x,y,I=niamey_cases1),from=c(tot1+1,1.5/S0[1]),
             to=c(2500,0.00125),sys3d="contour",n=c(40,40))
contour(c2$x,c2$y,c2$z,levels=-logLik(mle.fit2)+qchisq(0.95,2)/2, 
        col="red", add=TRUE)
lines(ellipse(vcov(mle.fit2),centre=coef(mle.fit2)),col=4)


###################################################
### chunk number 18: 
###################################################
p1 = profile(mle.fit2)


###################################################
### chunk number 19: plotprof
###################################################
plot(p1)


###################################################
### chunk number 20: 
###################################################
confint(mle.fit2)


###################################################
### chunk number 21: 
###################################################
confint(mle.fit2,method="quad")


###################################################
### chunk number 22: calcR0
###################################################
(R0.mle<-prod(coef(mle.fit2)))


###################################################
### chunk number 23: calcR0mat
###################################################
R0.fun<-function(S0,beta){
	S0*beta
	}
mmat<-outer(S0,beta,R0.fun)


###################################################
### chunk number 24: plotR0
###################################################
symbols(rep(S0,40),rep(beta,each=40),
        circles=50*as.vector(lik)/max(lik),xlab=~S[0],ylab=~beta,inches=FALSE)
contour(S0,beta,mmat,add=TRUE,col="blue",labcex=2)
## add half-integer values, with a dotted line
contour(S0,beta,mmat,add=TRUE,col="blue",lty=2,drawlabels=FALSE,
        levels=seq(1.5,12.5))
lines(ellipse(vcov(mle.fit2),centre=coef(mle.fit2)),col=5)


###################################################
### chunk number 25: 
###################################################
likelihood_R0 <- function(S0,R0,I){
    ## cat(S0,beta,"\n") ## debugging statement (commented out)
    n <- length(I)
    beta <- R0/S0
    S <- floor(S0-cumsum(I[-n]))
    p <- 1-exp(-beta*(I[-n]))	
    L <- dbinom(I[-1],S,p,log=TRUE)	
    Lik<-sum(-L,na.rm=TRUE)
}
R0start <- list(S0=coef(mle.fit2)["S0"],
                R0=coef(mle.fit2)["S0"]/coef(mle.fit2)["beta"])
mle.fit3 = mle2(start=R0start,
                 likelihood_R0,
                 method="L-BFGS-B",
                 lower=c(tot1,0.5),upper=c(10*tot1,10),
                 data=list(I=niamey_cases1),
                 control=list(ndeps=c(1,0.001)))
(c1 = confint(mle.fit3))


###################################################
### chunk number 26: plotR0prof
###################################################
p3 = plot(profile(mle.fit3,which="R0"))


###################################################
### chunk number 27: 
###################################################
curve3d(likelihood_R0(x,y,I=niamey_cases2),from=c(900,1),
        to=c(2000,3),sys3d="contour",n=c(40,40),
        levels=60:100)


###################################################
### chunk number 28: 
###################################################
mle.fit.s2 = mle2(start=list(S0=1100,R0=2),
                 likelihood_R0,
                 method="L-BFGS-B",
                 lower=c(tot1,0.5),upper=c(10*tot1,10),
                 data=list(I=niamey_cases2),
                 control=list(ndeps=c(1,0.001)))
c2 = confint(mle.fit.s2)
mle.fit.s3 = mle2(start=list(S0=1450,R0=2.2),
  likelihood_R0,
  method="L-BFGS-B",
  lower=c(tot1,0.5),upper=c(10*tot1,10),
  data=list(I=niamey_cases3),
  control=list(ndeps=c(1,0.001)))
c3 = confint(mle.fit.s3)


###################################################
### chunk number 29: 
###################################################
library(plotrix)
R0vals = c(coef(mle.fit3)[2],coef(mle.fit.s2)[2],coef(mle.fit.s3)[2])
cint = cbind(c1[2,],c2[2,],c3[2,])
suppressWarnings(plotCI(1:3,R0vals,li=cint[1,],ui=cint[2,],axes=FALSE,
       ylab=~R[0],xlab="Site"))
box()
axis(side=1,at=1:3)
axis(side=2)


###################################################
### chunk number 30:  eval=FALSE
###################################################
## likfun2<-function(R,vecs){
## ###############################
## #Profile likelihood for R0
## ###############################
##     vecs<<-vecs
##     likfun3<-function(S0){
##         beta<-R/S0
##         likelihood(c(S0,beta),vecs=vecs)
##     }
## 
##     o<-optimise(likfun3,c(sum(vecs$I),5*sum(vecs$I)))
##     r<-c(o$objective,o$minimum,R/o$minimum)
##     names(r)<-c("NLL","S0","beta")
##     r
## }


###################################################
### chunk number 31:  eval=FALSE
###################################################
## R0.vec<-seq(1,5,by=.01)
## lmat<-t(sapply(R0.vec,likfun2,vecs=list(I=niamey_cases1)))


###################################################
### chunk number 32:  eval=FALSE
###################################################
## par(mfrow=c(1,1))
## plot(R0.vec,lmat[,"NLL"],xlab="R0",ylab="NLL")
## critval<- mle.fit$value + qchisq(0.95,1)/2
## abline(h=critval)
## R0.mle<-prod(mle.fit$par)
## critfun<-function(m)likfun2(m,vecs=list(I=niamey_cases2))[1]-critval
## proflower<-uniroot(critfun,c(1,R0.mle))$root
## profupper<-uniroot(critfun,c(R0.mle,5))$root  
## abline(v=c(proflower,profupper))