Note: This tutorial assumes that you have successfully installed ABSEIR, and are at least passingly familiar with compartmental models. Some introductory information is available in this vignette.
Data for this example were obtained from the measlesWeserEms data set from the surveillance package.
# load the ABSEIR library
library(ABSEIR)
## Loading required package: Rcpp
##
## Attaching package: 'ABSEIR'
##
## The following object is masked _by_ '.GlobalEnv':
##
## WestAfrica2015
library(surveillance)
## Loading required package: sp
## Loading required package: xtable
## Loading required package: polyCub
## This is surveillance 1.11.0. For overview type 'help(surveillance)'.
# Read in data
data(measlesWeserEms)
# Identify cases
cases<-measlesWeserEms@observed
epidemic.start = min(which(apply(cases, 1, max) > 0))
cases = cases[(epidemic.start-1):nrow(cases),]
# Obtain distance matrix
neighbourhood<-measlesWeserEms@neighbourhood
# Week
week <- measlesWeserEms@epoch[(epidemic.start-1):length(measlesWeserEms@epoch)]
# Vaccination and population data set
vaccine.data <- measlesWeserEms@map@data
vaccine.data$adminID <- rownames(vaccine.data)
# Population size
N <- vaccine.data$POPULATION
# Check that spatial unit ordering makes sense
if (!all(vaccine.data$adminID == colnames(neighbourhood))){
stop("Error, make sure spatial unit ordering is consistent.")
}
weser.data_model = DataModel(Y=cases,
type = "identity", # Assume data is correct
compartment = "I_star", # Data related to new infections
cumulative = FALSE # Not reported on cumulative scale
)
I0 = apply(cases[1:3,], 2, max) > 0
E0 = I0
R0 = 0*I0
S0 = N-E0-I0-R0
weser.initial_values = InitialValueContainer(S0 = S0,
E0 = E0,
I0 = I0,
R0 = R0)
# No reinfection
weser.reinfection_model = ReinfectionModel("SEIR")
To specify the exposure process, we need to create the exposure design matrix. This matrix of temporal and spatial covariates must have a particular structure. Each of the \(J\) spatial locations receives a \(T\times p\) design matrix, where \(p\) is the number of covariates. These matrices are then 'stacked' row-wise in the same spatial ordering as presented in the columns and rows of the neighborhood matrix (defined above).
Here, we include the following covariates:
For more information on these covariates, see the documentation for the 'measlesWeserEMS' data.
As you can see, while ABSEIR doesn't directly fit SEIRV models (which incorporate a vaccination compartment), we can include vaccination data into the exposure probability structure, effectively giving vaccinated persons a different probability of contracting the pathogen of interest.
n.locations <- ncol(neighbourhood)
n.timepoints <- length(week)
time_invariant_covariates <- data.frame(intercept = rep(1, nrow(vaccine.data)),
popDensity = vaccine.data$POPULATION/vaccine.data$AREA,
proportionVaccineCard = vaccine.data$vaccdoc.2004,
proportionVaccine1 = vaccine.data$vacc1.2004,
proportionVaccine2 = vaccine.data$vacc2.2004)
time_varying_covariates <- data.frame(sin_component = sin(week/52*2*pi),
cos_component = cos(week/52*2*pi),
trig_interact = sin(week/52*2*pi)*cos(week/52*2*pi))
exposure.design.matrix <- as.matrix(
cbind(
time_invariant_covariates[rep(1:n.locations, each = n.timepoints),],
time_varying_covariates[rep(1:n.timepoints, n.locations),]
)
)
## Build the exposure model
weser.exposure_model <- ExposureModel(X = exposure.design.matrix,
nTpt = n.timepoints,
nLoc = n.locations,
betaPriorPrecision = rep(1, ncol(exposure.design.matrix)),
betaPriorMean = rep(0, ncol(exposure.design.matrix)))
## [1] "Assuming equally spaced count data."
We now need to specify the contact structures. Here, we'll build both a gravity model (which depends on population size and distance), and a simple CAR model (which just depends on shared borders between regions).
# Build a gravity model, in which the contact process between a pair of
# spatial locations is proportional to the product of their populations divided
# by the squared 'distance'
pop.matrix <- matrix(N, nrow = length(N), ncol = length(N))
unscaledGravityModel <- (pop.matrix * t(pop.matrix))/neighbourhood^2
diag(unscaledGravityModel) <- 1
diag(gravityModel) <- 0
weser.distance_model <- DistanceModel(list(gravityModel),
priorAlpha = 1,
priorBeta = 1,
scaleMode = "rowscale")
# Build a simpler contact model, in which the contact probability between
# a pair of spatial locations is only nonzero when the distance is equal to 1
# (CAR specification)
weser.CAR_model <- DistanceModel(list((neighbourhood == 1)),
priorAlpha = 1,
priorBeta = 1
)
# 9-12 day latent period
# Infectious
weser.transition_priors = ExponentialTransitionPriors(p_ei = 0.8, # Guess at E to I transition probability (per week)
p_ir = 0.8, # Guess at I to R transition probability (per week)
p_ei_ess = 10, # confidence
p_ir_ess = 10 # confidence
)
# Initial sampling step: simulate 100k epidemics from prior.
weser.sampling_control.init <- SamplingControl(seed=123123,
n_cores = 16,
algorithm = "Beaumont2009",
params = list(
batch_size = 500000,
epochs = 1,
shrinkage = 0.98,
max_batches = 1,
multivariate_perturbation=FALSE
))
# SMC Sampling configuration
weser.sampling_control <- SamplingControl(seed=123126,
n_cores = 14,
algorithm = "Beaumont2009",
params = list(
batch_size = 5000,
epochs = 1e6,
shrinkage = 0.99,
max_batches = 1000,
multivariate_perturbation=TRUE
))
t1 <- system.time(
{
weser.model1 <- SpatialSEIRModel(data_model = weser.data_model,
exposure_model = weser.exposure_model,
reinfection_model = weser.reinfection_model,
distance_model = weser.distance_model,
transition_priors = weser.transition_priors,
initial_value_container = weser.initial_values,
sampling_control = weser.sampling_control,
samples = 100,
verbose = FALSE)
weser.model1 <- update(weser.model1, weser.sampling_control,
verbose = FALSE)
save(weser.model1, file = "weserModel1.rda")
}
)
t2 <- system.time({
weser.model2 <- SpatialSEIRModel(data_model = weser.data_model,
exposure_model = weser.exposure_model,
reinfection_model = weser.reinfection_model,
distance_model = weser.CAR_model,
transition_priors = weser.transition_priors,
initial_value_container = weser.initial_values,
sampling_control = weser.sampling_control.init,
samples = 100,
verbose = FALSE)
weser.model2 <- update(weser.model2,
sampling_control = weser.sampling_control,
verbose=FALSE)
save(weser.model2, file = "weserModel2.rda")
}
)
summary(weser.model1)
## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Exposure Process Parameters: 8
## Reinfection Model Parameters: 0
## Spatial Parameters: 1
## Transition Parameters: 2
##
##
## Parameter Estimates:
## Mean SD 95% LB 95% UB
## Beta_SE_1 -2.031 1.852 -5.187 1.258
## Beta_SE_2 -0.135 2.388 -4.553 4.319
## Beta_SE_3 1.053 1.466 -2.169 3.555
## Beta_SE_4 -0.156 1.597 -3.383 2.518
## Beta_SE_5 -0.020 1.083 -2.025 1.929
## Beta_SE_6 1.384 0.814 -0.125 3.009
## Beta_SE_7 -0.494 0.856 -2.198 1.066
## Beta_SE_8 -1.264 1.020 -3.410 0.753
## rho_1 0.176 0.061 0.068 0.299
## gamma_EI 1.449 0.289 0.885 1.957
## gamma_IR 2.079 0.431 1.160 2.764
summary(weser.model2)
## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Exposure Process Parameters: 8
## Reinfection Model Parameters: 0
## Spatial Parameters: 1
## Transition Parameters: 2
##
##
## Parameter Estimates:
## Mean SD 95% LB 95% UB
## Beta_SE_1 -1.351 2.259 -5.952 1.924
## Beta_SE_2 0.491 2.082 -3.755 4.002
## Beta_SE_3 -0.209 1.697 -3.768 2.930
## Beta_SE_4 -1.313 1.981 -4.974 2.706
## Beta_SE_5 2.255 1.363 -0.440 4.701
## Beta_SE_6 1.358 0.490 0.482 2.280
## Beta_SE_7 -0.276 0.426 -0.969 0.638
## Beta_SE_8 -1.211 0.721 -2.704 0.048
## rho_1 0.032 0.012 0.013 0.055
## gamma_EI 1.995 0.466 1.197 3.053
## gamma_IR 1.385 0.285 0.898 1.954
makePlots <- function(modelObj, nm){
sims <- epidemic.simulations(modelObj, replicates = 10)
Is <- lapply(sims$simulationResults, function(x){x$I_star})
Is <- array(Reduce(c, Is), dim = c(nrow(Is[[1]]),
ncol(Is[[2]]),
length(Is)))
Ism <- apply(Is, 1:2, mean)
Islb <- apply(Is, 1:2, quantile, probs = c(0.025))
Isub <- apply(Is, 1:2, quantile, probs = c(0.975))
plotLocation <- function(x, model){
plot(cases[,x], ylim = c(0, max(Ism[,x])),
main = paste(model, ": Observed and Posterior Predictive Simulation.\n location ",
colnames(neighbourhood)[x], sep = ""))
apply(Is, 3, function(i){
lines(i[,x], col = rgb(0,0,0,0.1))
})
points(cases[,x], pch = 16, col = "blue")
}
for (i in 1:ncol(neighbourhood)){
plotLocation(i, nm)
}
}
makePlots(weser.model1, "Dist")
makePlots(weser.model2, "CAR")
compareModels(modelList = list(weser.model1, weser.model2), n_samples = 100)
## [1] "Assuming equal prior probabilities."
## [,1] [,2]
## [1,] 1.000 0.8888889
## [2,] 1.125 1.0000000