KikwitPredictionExample.R

library(ABSEIR)
# read in the data set
data(Kikwit1995)

lastTpt <- 80

# Create a model to relate observe data to epidemic process
count <- Kikwit1995$Count
count[(lastTpt+1):length(count)] <- NA
data_model = DataModel(count,
                       type = "identity",
                       compartment="I_star",
                       cumulative=FALSE)


intervention_term = cumsum(Kikwit1995$Date >  as.Date("05-09-1995", "%m-%d-%Y"))
exposure_model = ExposureModel(cbind(1,intervention_term),
                                 nTpt = nrow(Kikwit1995),
                                 nLoc = 1,
                                 betaPriorPrecision = 0.5,
                                 betaPriorMean = 0)
reinfection_model = ReinfectionModel("SEIR")
distance_model = DistanceModel(list(matrix(0)))
initial_value_container = InitialValueContainer(S0=5.36e6,
                                                E0=2,
                                                I0=2,
                                                R0=0)
transition_priors = ExponentialTransitionPriors(p_ei = 1-exp(-1/5), 
                                                p_ir= 1-exp(-1/7),
                                                p_ei_ess = 10,
                                                p_ir_ess = 10)
sampling_control = SamplingControl(seed = 123123, 
                                   n_cores = 8,
                                   algorithm="Beaumont2009",
                                   list(batch_size = 2000,
                                        init_batch_size = 1000,
                                        epochs = 1e6,
                                        max_batches = 100,
                                        shrinkage = 0.99,
                                        keep_compartments=TRUE,
                                        multivariate_perturbation=FALSE
                                   )
)
runtime = system.time(result <- SpatialSEIRModel(data_model,
                                                 exposure_model,
                                                   reinfection_model,
                                                   distance_model,
                                                   transition_priors,
                                                   initial_value_container,
                                                   sampling_control,
                                                   samples = 100,
                                                   verbose = 2))
simulations <- epidemic.simulations(result, replicates = 50)

plotPosteriorPredictive = function(simulations, rawData, main, lastTime)
{
  allSimulatedI_star = sapply(simulations$simulationResults, function(x){x$I_star})
  
  lowerQuantile = apply(allSimulatedI_star, 1, quantile, probs = c(0.025))
  posteriorMean = apply(allSimulatedI_star, 1, mean)
  upperQuantile = apply(allSimulatedI_star, 1, quantile, probs = c(0.975))
  
  
  plot(rawData, ylim = c(0, max(rawData)*2),
       xlab = "Epidemic Day", ylab = "New Cases", main = main,
       col = ifelse(1:length(rawData) <= lastTime, "black", "red"))
  lines(upperQuantile, lty = 2, col = "blue")
  lines(lowerQuantile, lty = 2, col = "blue")
  lines(posteriorMean, lty = 1, col = "blue")
  
  legend(x = 100, y = 12, legend = c("Mean", "95% CI", "Observed", "Future"), lty = c(1,2,0,0), 
         pch = c(NA,NA,1,1), col = c("blue", "blue", "black","red"), cex = 1)
}

plotPosteriorPredictive(result, Kikwit1995$Count, "Model 1: Posterior Distribution", lastTpt)
plotPosteriorPredictive(simulations, Kikwit1995$Count, "Model 1: Posterior Predictive Distribution", lastTpt)