## @knitr set-opts
##glop <- options(keep.source=TRUE,width=80,continue=" ",prompt=" ")
glop <- options()
pdf.options(useDingbats=FALSE)
set.seed(5384959)
require(ggplot2)
require(plyr)


## @knitr load-bbs-data-alt

baseURL <- "http://www.math.mcmaster.ca/bolker/eeid/data"
flu <- read.csv(url(paste(baseURL,"boarding_school_flu.csv",sep="/")))



## @knitr load-bbs-data

flu <- read.csv("boarding_school_flu.csv")



## @knitr bbs-ggplot

require(ggplot2)
ggplot(data=flu,mapping=aes(x=day,y=flu))+geom_point(size=2)+geom_line()



## @knitr bbs-plot


require(ggplot2)
ggplot(data=flu,mapping=aes(x=day,y=flu))+geom_point(size=2)+geom_line()




## @knitr unnamed-chunk-1

require(deSolve)



## @knitr open-sir-model-defn

sibr.model <- function (t, x, params) {
  ## first extract the state variables
  S <- x[1]
  I <- x[2]
  B <- x[3]
  ## now extract the parameters
  beta <- params["beta"]
  gamma <- params["gamma"]
  delta <- params["delta"]
  N <- 763
  ## now code the model equations
  dS.dt <- -beta*S*I/N
  dI.dt <- beta*S*I/N-gamma*I
  dB.dt <- gamma*I-delta*B
  ## combine results into a single vector
  dxdt <- c(dS.dt,dI.dt,dB.dt) 
  ## return result as a list!
  list(dxdt)                              
}



## @knitr set-closed-params

params <- c(beta=2,gamma=1/3,delta=1/3)



## @knitr set-times-ics

times <- seq(from=0,to=15,by=1/4) ## returns a sequence
xstart <- c(S=762,I=1,B=0)       ## initial conditions



## @knitr unnamed-chunk-2

out <- ode(
           func=sibr.model,
           y=xstart,
           times=times,
           parms=params
           )
class(out)
head(out)



## @knitr solve-closed-sibr

out <- as.data.frame(
                     ode(
                         func=sibr.model,
                         y=xstart,
                         times=times,
                         parms=params
                         )
                     )



## @knitr epi-curve-ggplot
require(ggplot2)
require(reshape)
out <- subset(out,select=c(I,B,time))
ggplot(data=melt(out,id.var="time"),
       mapping=aes(x=time,y=value,group=variable,color=variable))+
  geom_line()


## @knitr epi-curve-plot

require(ggplot2)
require(reshape)
out <- subset(out,select=c(I,B,time))
ggplot(data=melt(out,id.var="time"),
       mapping=aes(x=time,y=value,group=variable,color=variable))+
  geom_line()



## @knitr full-sim

times <- seq(from=0,to=14,by=1)
params <- c(beta=2,gamma=1/3,delta=1/3,N=763,p=0.3,k=1000)
out <- as.data.frame(
                     ode(
                         func=sibr.model,
                         y=xstart,
                         times=times,
                         parms=params
                         )
                     )
within(
       out,
       C <- rnbinom(n=length(B),mu=params["p"]*B,size=params["k"])
       ) -> out
ggplot(data=out,mapping=aes(x=time,y=C))+geom_point()+geom_line()



## @knitr sir-model-likelihood

sibr.nll <- function (beta, gamma, delta, I.0, p, k) {
  times <- c(flu$day[1]-1,flu$day)
  ode.params <- c(beta=beta,gamma=gamma,delta=delta)
  xstart <- c(S=763-I.0,I=I.0,B=0)
  out <- ode(
             func=sibr.model,
             y=xstart,
             times=times,
             parms=ode.params
             )
  ## 'out' is a matrix
  ll <- dnbinom(x=flu$flu,size=params["k"],mu=p*out[-1,"B"],log=TRUE)
  -sum(ll)
}



## @knitr beta-loglik

nll <- function (par) {
  sibr.nll(beta=par[1],gamma=1/3,delta=1/3,
             I.0=1,p=0.5,k=100)
}
betacurve <- data.frame(beta=seq(1/3,10,length=100))
within(betacurve,nll <- sapply(beta,nll)) -> betacurve
ggplot(data=betacurve,mapping=aes(x=beta,y=nll))+geom_line()



## @knitr optim-beta

fit <- optim(fn=nll,par=2,method="Brent",lower=1.5,upper=3)
fit



## @knitr beta-lr-ci

crit.lr <- pchisq(q=0.05,df=1,lower.tail=FALSE)
betacurve <- data.frame(beta=seq(2,3,length=100))
betacurve <- within(betacurve,nll <- sapply(beta,nll))
ggplot(data=betacurve,mapping=aes(x=beta,y=nll))+geom_line()+
  ylim(fit$value+c(0,10))+
  geom_vline(xintercept=fit$par,color='red')+
  geom_hline(yintercept=fit$value+crit.lr,color='blue')



## @knitr beta-gamma-fit

nll <- function (par) {
  sibr.nll(beta=par[1],gamma=par[2],delta=1/3,
             I.0=1,p=0.5,k=100)
}
fit <- optim(fn=nll,par=c(2.6,1/3),method="Nelder-Mead")
fit



## @knitr mle2-beta-gamma

require(bbmle)
fit <- mle2(sibr.nll,start=list(beta=2.5,gamma=0.7),
            method="Nelder-Mead",
            fixed=list(delta=1/3,k=100,p=0.5,I.0=1))
coef(fit)
pfit <- profile(fit)
plot(pfit)
confint(pfit)



## @knitr mle2-beta-gamma-p-k

require(bbmle)
fit <- mle2(sibr.nll,
            start=list(beta=1.5,gamma=1/3,p=0.5),
            method="L-BFGS-B",
            lower=c(0,0,0),
            upper=c(Inf,Inf,1),
            fixed=list(delta=1/3,k=100,I.0=1))
coef(fit)
pfit <- profile(fit)
plot(pfit)
confint(pfit)



## @knitr mle-sim

mle <- coef(fit)

times <- seq(from=0,to=14,by=1)
xstart <- c(S=763-mle["I.0"],I=mle["I.0"],B=0)
out <- as.data.frame(
                     ode(
                         func=sibr.model,
                         y=xstart,
                         times=times,
                         parms=mle
                         )
                     )
within(
       subset(out,time>0),
       C <- rnbinom(n=length(B),mu=mle["p"]*B,size=mle["k"])
       ) -> out

ggplot(data=out,mapping=aes(x=time,y=C))+geom_line()+
  geom_point(data=flu,mapping=aes(x=day,y=flu))




## @knitr get-plague-data
plague <- read.csv(url(paste(baseURL,"bombay_plague.csv",sep="/")))


## @knitr restore-opts
options(glop)