In grantbrown/ABSEIR: Posterior Approximation to Fully Bayesian Spatial SEIR Models Using Multithreaded C++

Introduction

This document presents an analysis of the 2014 and 2015 Ebola epidemic in West Africa, centered on Guinea, Liberia and Sierra Leone. More information is available on Wikipedia.

Note: This tutorial is incomplete, containing example code without sufficient explanation. Documentation efforts are underway.

library(ABSEIR)
library(splines)

data(WestAfrica2015)
WestAfrica2015 = WestAfrica2015[rev(1:nrow(WestAfrica2015)),]
timeIdx = as.Date(WestAfrica2015$Date, format = "%m/%d/%Y")
timeIdx = as.numeric(timeIdx - min(timeIdx)) + 1
modelIdx = 1:length(timeIdx)


I_star = matrix(NA, nrow = max(timeIdx), ncol = 3)
I_star[timeIdx[modelIdx], 1] = WestAfrica2015$Cases_Guinea[modelIdx]
I_star[timeIdx[modelIdx], 2] = WestAfrica2015$Cases_Liberia[modelIdx]
I_star[timeIdx[modelIdx], 3] = WestAfrica2015$Cases_SierraLeone[modelIdx]

# Fill in initial zeros
I_star[1:(min(which(!is.na(I_star[,2])))-1),2] <- 0
I_star[1:(min(which(!is.na(I_star[,3])))-1),3] <- 0


# Linearly interpolate missing data:
I_star <- apply(I_star, 2, function(i){round(approx(1:length(i), i, n = length(i))$y)})

# Make sure these make sense as cumulative counts (sometimes decreases were observed)
currentVals = I_star[nrow(I_star),]
for (i in (nrow(I_star)-1):1)
{
  badIdx = I_star[i,] > currentVals 
  badIdx = ifelse(is.na(badIdx), FALSE, badIdx)
  I_star[i, ] = ifelse(badIdx, currentVals, I_star[i,])
  currentVals = ifelse(is.na(I_star[i,]), currentVals, I_star[i,])
}

# Set up starting val
I0 <- I_star[1,]
# Thin the data to weekly:
I_star <- I_star[seq(2,nrow(I_star),7),]

data_model = DataModel(Y = I_star,
                       type = "identity",
                       compartment = "I_star",
                       cumulative = TRUE)

intercepts = diag(3)[rep(1:ncol(I_star), each = nrow(I_star)),]
#intercepts = 1
timeBasis = bs(1:nrow(I_star), degree = 3)[rep(1:nrow(I_star), ncol(I_star)),]
X = cbind(intercepts, timeBasis)
exposure_model = ExposureModel(X, nTpt = nrow(I_star),
                               nLoc = ncol(I_star),
                               betaPriorPrecision = 0.5,
                               betaPriorMean = c(rep(-1, ncol(intercepts)),
                                                 rep(0, ncol(timeBasis))))

reinfection_model = ReinfectionModel("SEIR")
DM1 = matrix(c(0,1,0,
               1,0,1,
               0,1,0), nrow = 3, byrow = TRUE)

distance_model = DistanceModel(list(DM1), priorAlpha = 1, priorBeta = 25)
N = c(10057975, 4128572, 6190280)
E0 = apply(I_star[1:4,], 2, sum, na.rm = TRUE)
initial_value_container = InitialValueContainer(S0=N - I0 - E0,
                                                E0 = E0,
                                                I0 = I0,
                                                R0 = rep(0, ncol(I_star)))

transition_priors1 = ExponentialTransitionPriors(p_ei = 1-exp(-1/5), 
                                     p_ir= 1-exp(-1/7),
                                     p_ei_ess = 100,
                                     p_ir_ess = 100)




sampling_control = SamplingControl(seed = 123124, 
                                   n_cores = 14,
                                   algorithm="Beaumont2009",
                                   list(batch_size = 2500,
                                           epochs = 1e6,
                                           #max_batches = 2000,
                                           max_batches = 100,
                                           shrinkage = 0.99,
                                           multivariate_perturbation=FALSE
                                         )
                                   )

## Exponential Results:
system.time(result <- SpatialSEIRModel(data_model,
                          exposure_model,
                          reinfection_model,
                          distance_model,
                          transition_priors1,
                          initial_value_container,
                          sampling_control,
                          samples = 50,
                          verbose = 0))

sims = epidemic.simulations(result, replicates = 25)
sim_I_star = lapply(sims$simulationResults, function(x){apply(x$I_star, 2, cumsum)})
ism = array(Reduce(c, sim_I_star), dim = c(nrow(sim_I_star[[1]]),
                                           ncol(sim_I_star[[2]]),
                                           length(sim_I_star)))
ismLB = apply(ism,1:2,quantile,probs = c(0.05))
ismUB = apply(ism,1:2,quantile,probs = c(0.95))
ismMean = apply(ism,1:2,mean)

casePlot = function(idx, main){
  plot(I_star[,idx], ylim = c(0, max(ismUB)), 
       ylab = "Total Cases", xlab = "Time",
       main = main)
  lines(ismMean[,idx], col = "blue", lty = 2, lwd = 2)
  lines(ismLB[,idx], col = "blue", lty = 3, lwd = 1)
  lines(ismUB[,idx], col = "blue", lty = 3, lwd = 1)
}

casePlot(1, "Guinea - Exponential ")
casePlot(2, "Liberia - Exponential ")
casePlot(3, "Sierra Leone - Exponential ")

R0.result <- ComputeR0(sims, cores = 8)


sim_EAR0 = lapply(R0.result$simulationResults, function(x){x$R_EA})
eaR0ar = array(Reduce(c, sim_EAR0), dim = c(nrow(sim_I_star[[1]]),
                                           ncol(sim_I_star[[2]]),
                                           length(sim_I_star)))

EAR0_mean = apply(eaR0ar,1:2,mean)
EAR0_LB = apply(eaR0ar,1:2,quantile, probs = 0.025)
EAR0_UB = apply(eaR0ar,1:2,quantile, probs = 0.975)
plotRepNum <- function(idx, lbl){
  plot(EAR0_mean[,idx], type = "l", ylim = c(0, max(EAR0_UB[,idx])), col = "blue", lwd = 2,
       ylab = "Empirically Adjusted R0",
       xlab = "Time Index", main = lbl)
  lines(EAR0_UB[,idx], col = "blue", lty = 2)
  lines(EAR0_LB[,idx], col = "blue", lty = 2)
  abline(h = seq(0,100,0.5), lty = 2, col = "lightgrey")
  abline(h = 1, lty = 2)
}

plotRepNum(1, "Guinea")
plotRepNum(2, "Liberia")
plotRepNum(3, "Sierra Leone")

grantbrown/ABSEIR documentation built on Oct. 14, 2021, 2:32 p.m.