###########################################################################
#                                                                         #
#   RAPID DROP IN THE REPRODUCTION NUMBER DURING THE EBOLA OUTBREAK       #
#   IN THE DEMOCRATIC REPUBLIC OF CONGO                                   #
#                                                                         #
#   REFERENCE: ALTHAUS CL. PEERJ. 2015;3:e1418.                           #
#   DOI: 10.7717/peerj.1418                                               #
#                                                                         #
#   Copyright (C) 2015 by Christian L. Althaus                            #
#   (christian.althaus@alumni.ethz.ch)                                    #
#                                                                         #
#   This script is free software; you can redistribute it and/or modify   #
#   it under the terms of the GNU General Public License as published by  #
#   the Free Software Foundation; either version 2 of the License, or     #
#   (at your option) any later version.                                   #
#                                                                         #
#   This script is distributed in the hope that it will be useful,        #
#   but WITHOUT ANY WARRANTY; without even the implied warranty of        #
#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         #
#   GNU General Public License for more details.                          #
#                                                                         #
###########################################################################

# init
rm(list = ls())

# Necessary libraries
library(chron)
library(R0)

# Initialize random number generator
set.seed(679236)

# Read the data
ebola <- read.csv("Ebola_outbreak_DRC_data.csv")
ebola$Date <- chron(as.character(ebola$Date), format=c(dates = "day mon year"))
ebola$Date <- as.Date(ebola$Date)
ebola <- ebola[1:71,]

# Estimating case reproduction number
mGT <- generation.time("gamma", c(15.3,9.3)) # Generation time (serial interval) as reported by WHO Ebola Response Team (2014, NEJM)
TD <- est.R0.TD(ebola$Cases, mGT, t = ebola$Date, begin = 1, end = as.numeric(max(ebola$Date)-min(ebola$Date)), nsim = 1e4, correct = FALSE)

# Plot the case reproduction number during the outbreak
l <- length(ebola$Date)-1
plot(TD)
plot(NA,ty="n",frame=FALSE,axes=FALSE,xlim=c(as.Date("2014-07-26"),as.Date("2014-10-18")),ylim=c(0,8),xlab=NA,ylab=bquote("Case reproduction number " ~ italic("R")),cex.sub = 0.75)
polygon(x = c(ebola$Date[1:l], rev(ebola$Date[1:l])), y = c(TD$conf.int[1:l,1], rev(TD$conf.int[1:l,2])), col = rgb(0, 0, 1, alpha=0.2), border = NA)
lines(ebola$Date,c(TD$R,NA), col = "blue")
axis(1, seq(as.Date("2014-07-26"),as.Date("2014-10-18"),14), c("26 Jul 2014",NA,"23 Aug 2014",NA,"20 Sep 2014",NA,"18 Oct 2014"))
axis(2, 0:8)
control <- as.Date("2014-08-16")
lines(c(0,control),c(1,1))
lines(c(control,control),c(0,1))
points(control,1,pch=19)
