R ABSEIR - R0 and EAR0 Issue.Rmd

---
title: "R ABSEIR - R0 and EA-R0 issue"
output: html_document
---

## Issue

Modified code below 

```{r Setup, echo=TRUE}
# load the ABSEIR library
library(ABSEIR)

# read in the data set
data(Kikwit1995)

# Create a model to relate observe data to epidemic process
data_model = DataModel(Kikwit1995$Count,
                       type = "identity",
                       compartment="I_star",
                       cumulative=FALSE)

# Create a model to describe relationship of any covariates to epidemic intensity
# In this analysis we know on which date interventions began, so we may include a linear term
# starting on that date. 
intervention_term = cumsum(Kikwit1995$Date >  as.Date("05-09-1995", "%m-%d-%Y"))
exposure_model_2 = ExposureModel(cbind(1,intervention_term),
                                 nTpt = nrow(Kikwit1995),
                                 nLoc = 1,
                                 betaPriorPrecision = 0.5,
                                 betaPriorMean = 0)

# There's no reinfection in this case, so we just use a "SEIR" model. 
reinfection_model = ReinfectionModel("SEIR")

distance_model = DistanceModel(list(matrix(0)))

# Set initial population sizes
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 = 100,
                                                p_ir_ess = 100)

# Set algorithm configuration
sampling_control = SamplingControl(seed = 123123, 
                                   n_cores = 8, # changed since I have fewer cores
                                   algorithm="Beaumont2009",
                                   list(init_batch_size = 10000,
                                        batch_size = 2000,
                                        epochs = 1e6,
                                        max_batches = 100,
                                        shrinkage = 0.99,
                                        multivariate_perturbation=FALSE ,
                                        keep_compartments=TRUE
                                   )
)

# Reasonable intensity
runtime2 = system.time(result2 <- SpatialSEIRModel(data_model,
                                                   exposure_model_2,
                                                   reinfection_model,
                                                   distance_model,
                                                   transition_priors,
                                                   initial_value_container,
                                                   sampling_control,
                                                   samples = 100,
                                                   verbose = FALSE))

simulations = epidemic.simulations(result2, replicates = 50)

simulations2 <- ComputeR0(simulations, cores = 7)

plotR0 = function(simulations, main)
{
  #allSimulatedR0 = sapply(simulations2$simulationResults, function(x){x$R0t})
  allSimulatedEA_R0 = sapply(simulations$simulationResults, function(x){x$R_EA})
  
  plot(apply(allSimulatedEA_R0, 1, mean), type = "l", ylim = c(0, 3), lwd =2,
       ylab = "Reproductive Number", main = main)
  
  lines(apply(allSimulatedEA_R0, 1, quantile, probs = 0.025), type = "l", ylim = c(0, 3), 
       lwd =2,lty=2)
  lines(apply(allSimulatedEA_R0, 1, quantile, probs = 0.975), type = "l", ylim = c(0, 3), 
       lwd =2,lty=2)
}

# error pops up when running this line
plotR0(simulations2, "Model 2: R0(t) and EA-R(t)")


```

Traditional R0 isn't estimated automatically anymore, as it isn't defined unambiguously for all the models implemented by ABSEIR. We can calculate it for this model though:

```{r}
params <- as.data.frame(result2$param.samples)
mean((exp(params$Beta_SE_1)/params$gamma_IR))
quantile((exp(params$Beta_SE_1)/params$gamma_IR), probs = c(0.05,0.95))

```