########################################################################### # # # RAPID DROP IN THE REPRODUCTION NUMBER DURING THE EBOLA OUTBREAK # # IN THE DEMOCRATIC REPUBLIC OF CONGO # # # # REFERENCE: ALTHAUS CL. PEERJ PREPRINTS. 2015;3:E1279. # # DOI: 10.7287/peerj.preprints.1041v1 # # # # 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(deSolve) library(chron) library(bbmle) library(mvtnorm) # Initialize random number generator set.seed(529342) # 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$Cases[1] <- 0 # Definition of the intervention function exponential <- function(t,beta0,beta1,k,tau1) { if(t < tau1) { beta <- beta0 } else { beta <- beta1 + (beta0-beta1)*exp(-k*(t-tau1)) } return(beta) } Exponential <- Vectorize(exponential, vectorize.args = 't') # Definition of the SEIR model SEIR <- function(t, x, parms) { with(as.list(c(parms,x)),{ beta <- exponential(t,beta0,beta1,k,tau1) # Density-dependent transmission N <- S + E + I + R + D dS <- - beta/N*S*I dE <- beta/N*S*I - sigma*E dI <- sigma*E - gamma*I dR <- (1-f)*gamma*I dD <- f*gamma*I dC <- sigma*E der <- c(dS,dE,dI,dR,dD,dC) list(der) }) } # Negative log-likelihood nll <- function(beta0,beta1,k,f,tau0,tau1,sigma,gamma) { pars <- c(beta0=beta0,beta1=beta1,k=k,f=f,tau0=tau0,tau1=tau1,sigma=sigma,gamma=gamma) pars <- trans(pars) times <- c(data$times+pars["tau0"],max(data$times+pars["tau0"])+1) simulation <- as.data.frame(ode(init,times,SEIR,parms=pars)) ll <- sum(dpois(data$cases,diff(simulation$C),log=TRUE)) return(-ll) } # Parameter transformation trans <- function(pars) { pars["beta0"] <- exp(pars["beta0"]) pars["beta1"] <- exp(pars["beta1"]) pars["k"] <- exp(pars["k"]) pars["f"] <- plogis(pars["f"]) pars["tau0"] <- exp(pars["tau0"]) if(is.na(pars["tau1"])) { pars["tau1"] <- pars["tau0"] } return(pars) } # Prepare the data and set the initial values data <- na.omit(ebola[c("Date","Cases")]) names(data) <- c("times","cases") begin <- chron("26 Jul 2014", format=c(dates = "day mon year")) data$times <- data$times - data$times[1] N <- 1e6 init <- c(S = N - 1, E = 0, I = 1, R = 0, D = 0, C = 0) # Fit the model to the data fixed <- c(sigma = 1/9.312799, gamma = 1/7.411374, tau0 = -Inf, beta1 = -Inf, f = qlogis(49/69)) free <- c(beta0 = log(1), k = log(1), tau1 = 14) fit <- mle2(nll,start=as.list(free),fixed=as.list(fixed),method="Nelder-Mead",control=list(maxit=1e3)) fit # Parameter sampling (for CIs and prediction interval) n_sim <- 1e4 m <- coef(fit,exclude.fixed=TRUE) sigma <- vcov(fit) sim_coef <- data.frame(rmvnorm(n_sim, mean = m, sigma = sigma)) # Timepoints pars <- trans(coef(fit)) end <- chron("18 Oct 2014", format=c(dates = "day mon year")) weeks <- seq.dates(from = chron("26 Jul 2014", format = c(dates = "day mon year")), to = chron("18 Oct 2014", format = c(dates = "day mon year")), by = 14) timepoints <- seq(0,end-begin,1) # Plot the daily incidence and cumulative number of new cases par(mfrow=c(1,2)) # Simulations mod_sim_Inc <- numeric() obs_sim_Inc <- numeric() mod_sim_Cum <- numeric() obs_sim_Cum <- numeric() for(i in 1:n_sim) { pars <- trans(c(unlist(sim_coef[i,]),fixed)) simulation <- as.data.frame(ode(init,timepoints,SEIR,parms=pars)) mod_sim_Inc <- cbind(mod_sim_Inc, c(0,diff(simulation$C))) mod_sim_Cum <- cbind(mod_sim_Cum, simulation$C) obs_sim_dC <- rpois(length(timepoints) - 1, lambda = diff(simulation$C)) obs_sim_Inc <- cbind(obs_sim_Inc, c(0, obs_sim_dC)) obs_sim_Cum <- cbind(obs_sim_Cum, cumsum(c(0, obs_sim_dC))) } # Calculate the point-wise quantiles to construct an approximate 95% CI and a 95% prediction interval (cases). mod_sim <- mod_sim_Inc obs_sim <- obs_sim_Inc x_grid <- timepoints mod_q <- apply(mod_sim, MAR = 1, FUN = quantile, probs = c(0.025, 0.975)) obs_q <- apply(obs_sim, MAR = 1, FUN = quantile, probs = c(0.025, 0.975)) simulation <- as.data.frame(ode(init,x_grid,SEIR,parms=trans(coef(fit)))) plot(x_grid, mod_q[2,], ylim = c(0,7), ty="l", lty=2, lwd = 1, col = "red", xlab=NA, ylab="Daily incidence of new cases",frame=FALSE,axes=FALSE) axis(1, weeks - begin, weeks) axis(2) polygon(x = c(x_grid, rev(x_grid)), y = c(obs_q[1,], rev(obs_q[2,])), col = rgb(1, 0, 0, alpha=0.2), border = NA) lines(x_grid, mod_q[1,], lty = 2, lwd = 1, col = "red") lines(x_grid, c(0,diff(simulation$C)), lwd = 1, col = "red") points(data$times[1:85],data$cases[1:85],col="red") mod_sim <- mod_sim_Cum obs_sim <- obs_sim_Cum x_grid <- timepoints mod_q <- apply(mod_sim, MAR = 1, FUN = quantile, probs = c(0.025, 0.975)) obs_q <- apply(obs_sim, MAR = 1, FUN = quantile, probs = c(0.025, 0.975)) simulation <- as.data.frame(ode(init,x_grid,SEIR,parms=trans(coef(fit)))) plot(x_grid, 1+mod_q[2,], ylim = c(0,100), ty="l", lty=2, lwd = 1, col = "red", xlab=NA, ylab="Cumulative number of cases",frame=FALSE,axes=FALSE) axis(1, weeks - begin, weeks) axis(2) polygon(x = c(x_grid, rev(x_grid)), y = c(1+obs_q[1,], rev(1+obs_q[2,])), col = rgb(1, 0, 0, alpha=0.2), border = NA) lines(x_grid, 1+mod_q[1,], lty = 2, lwd = 1, col = "red") lines(x_grid, 1+simulation$C, lwd = 1, col = "red") points(data$times[1:85],1+cumsum(data$cases)[1:85],col="red") # Calculate the effective reproduction number vs. time (including approximate pointwise CI) beta <- function(t, pars) Exponential(t, pars[["beta0"]],pars[["beta1"]],pars[["k"]],pars[["tau1"]]) pars <- trans(coef(fit)) timepoints <- seq(0,end-begin,0.1) Re <- beta(timepoints, pars) / pars[["gamma"]] Re_sim <- numeric() for(i in 1:n_sim) { pars_i <- trans(c(unlist(sim_coef[i,]), fixed)) Re_i <- beta(timepoints, pars_i) / pars_i[["gamma"]] Re_sim <- cbind(Re_sim, Re_i) } Re_q <- apply(Re_sim, MAR = 1, FUN = quantile, probs = c(0.025, 0.975)) # Plot the effective reproduction number par(mfrow=c(1,1)) control <- (min(which(Re<1)) - 1)/10 control_l <- (min(which(Re_q[1,]<1)) - 1)/10 control_u <- (min(which(Re_q[2,]<1)) - 1)/10 plot(timepoints,Re,type="n",col="blue",ylim=c(0,8),xlab=NA,ylab=bquote("Net reproduction number " ~ italic("R"["t"])),axes=FALSE,frame=FALSE) polygon(x = c(timepoints, rev(timepoints)), y = c(Re_q[1,], rev(Re_q[2,])), col = rgb(0, 0, 1, alpha=0.2), border = NA) lines(timepoints, Re, col = "blue") lines(timepoints, Re_q[1,], col = 'blue', lty = 2) lines(timepoints, Re_q[2,], col = 'blue', lty = 2) axis(1, weeks - begin, weeks) axis(2, 0:8) lines(c(0,control),c(1,1)) lines(c(control,control),c(0,1)) points(control,1,pch=19)