Measles Example
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.
Step 0: Setup
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## Loading required package: parallel## Loading required package: compilerlibrary(surveillance)## Loading required package: sp## Loading required package: xtable## This is surveillance 1.19.1. 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.")
}Step 1: Data Model and Initial Values
# Make cumulative
weser.data_model = DataModel(Y=apply(cases, 2,cumsum),
type = "identity", # Assume data is correct
compartment = "I_star", # Data related to new infections
cumulative = TRUE # Not reported on cumulative scale
)## Warning in if (class(Y) != "matrix") {: the condition has length > 1 and only
## the first element will be usedI0 = (apply(cases[1:3,], 2, max) > 0)*2
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")Step 2: Exposure Process
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:
- An intercept term
- The population density
- The proportion of students with a vaccine card
- The proportion of students receiving at least one vaccination
- The proportion of students who were doubly vaccinated
- A three degree of freedom trigonometric temporal basis, to capture seasonality
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."Step 3: Contact Structure
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))
gravityModel <- (pop.matrix * t(pop.matrix))/neighbourhood^2
diag(gravityModel) <- 1
# Rescale
maxRowSum <- max(apply(gravityModel,1,sum))
gravityModel <- gravityModel/maxRowSum
weser.distance_model <- DistanceModel(list(gravityModel),
priorAlpha = 1,
priorBeta = 1)
# 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)*1),
priorAlpha = 1,
priorBeta = 1
)Step 4: Specify Latent and Infectious Period Transition Models
# 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
)Step 5: Configure ABC Sampler
weser.sampling_control <- SamplingControl(seed=123124,
n_cores = 14,
algorithm = "Beaumont2009",
params = list(
batch_size = 2000,
init_batch_size = 1000000,
epochs = 1e6,
shrinkage = 0.99,
max_batches = 200,
multivariate_perturbation=FALSE
))Step 6: Run Our Two Models
# Cache the model objects so we don't need to re-run the analysis for things like format tweaks
if (!file.exists( "weserModel1.rda")){
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 = 50,
verbose = FALSE)
}
)
save("weser.model1", "t1", file = "weserModel1.rda")
} else {
load("weserModel1.rda")
}
if (!file.exists( "weserModel2.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,
samples = 50,
verbose = 2)
}
)
save("weser.model2", "t2", file = "weserModel2.rda")
} else {
load("weserModel2.rda")
}## Initializing Model Components
## ...Building data model
## ...Building distance model
## ...Building exposure model
## ...Building initial value container
## ...Building reinfection model
## ...Building sampling control model
## ...Building transition priors
## ...Preparing model object
## Running main simulation
## Starting sampler
## Number of iterations requested: 1000000
## Data Model Summary
## ------------------
## Number of locations: 17
## Number of time points: 96
## Data Model Compartment: I_star
## Data Model Type: identity
##
## Exposure Model Summary
## ---------------------- dim(X): (1632, 8)
## beta prior means: (0,
## 0,
## 0,
## 0,
## 0,
## 0,
## 0,
## 0)
## beta prior precision: (1,
## 1,
## 1,
## 1,
## 1,
## 1,
## 1,
## 1)
##
## Transition Priors Summary
## -------------------------
## transition mode: exponential
## E to I transition prior params: (10, 6.21335)
## I to R transition prior params: (10, 6.21335)
## Reinfection Model Summary
## -------------------------
## reinfection mode: SEIR
##
## Distance Model Summary
## ----------------------
## Number of locations: 17
## Number of distance structures: 1
## Number of time varying distance structures: 96
## Number of time varying distance lags: 0
## Initial Compartment Values - Summary
## ------------------------------------
## S0, E0, I0, R0
## 75986, 0, 0, 0
## 51445, 0, 0, 0
## 158340, 0, 0, 0
## 165517, 0, 0, 0
## 84586, 0, 0, 0
## 114524, 0, 0, 0
## 189648, 2, 2, 0
## 153283, 0, 0, 0
## 307734, 0, 0, 0
## 101657, 0, 0, 0
## 132975, 0, 0, 0
## 164536, 2, 2, 0
## 124564, 0, 0, 0
## 358041, 0, 0, 0
## 130471, 0, 0, 0
## 94242, 0, 0, 0
## 57672, 0, 0, 0
##
## Sampling Control Summary
## ------------------------
## algorithm: 2
## simulation_width: 1
## random_seed: 123124
## CPU_cores: 14
## init_batch_size: 1000000
## batch_size: 2000
## epochs: 1000000
## max_batches: 200
## multivariatePerturbation: 0
## m: 1
## accept_fraction: 1
## shrinkage: 0.99
## lpow: 2
## target_eps: 0
## Note: not all parameters are used for all algorithms.
##
## Generating starting parameters from prior
## Iteration 0 [2458.9, 3186.43] eps: 3187.43
## batch 0, 0/50 accepted
## batch 1, 0/50 accepted
## batch 2, 0/50 accepted
## batch 3, 1/50 accepted
## batch 4, 3/50 accepted
## batch 5, 3/50 accepted
## batch 6, 4/50 accepted
## batch 7, 7/50 accepted
## batch 8, 9/50 accepted
## batch 9, 11/50 accepted
## batch 10, 12/50 accepted
## batch 11, 15/50 accepted
## batch 12, 15/50 accepted
## batch 13, 15/50 accepted
## batch 14, 15/50 accepted
## batch 15, 15/50 accepted
## batch 16, 18/50 accepted
## batch 17, 18/50 accepted
## batch 18, 19/50 accepted
## batch 19, 20/50 accepted
## batch 20, 22/50 accepted
## batch 21, 23/50 accepted
## batch 22, 26/50 accepted
## batch 23, 26/50 accepted
## batch 24, 27/50 accepted
## batch 25, 27/50 accepted
## batch 26, 29/50 accepted
## batch 27, 31/50 accepted
## batch 28, 31/50 accepted
## batch 29, 33/50 accepted
## batch 30, 33/50 accepted
## batch 31, 33/50 accepted
## batch 32, 34/50 accepted
## batch 33, 34/50 accepted
## batch 34, 36/50 accepted
## batch 35, 38/50 accepted
## batch 36, 40/50 accepted
## batch 37, 41/50 accepted
## batch 38, 41/50 accepted
## batch 39, 41/50 accepted
## batch 40, 42/50 accepted
## batch 41, 43/50 accepted
## batch 42, 43/50 accepted
## batch 43, 44/50 accepted
## batch 44, 47/50 accepted
## batch 45, 47/50 accepted
## batch 46, 48/50 accepted
## batch 47, 49/50 accepted
## Iteration 1 [2283.95, 3147.19] eps: 3155.56
## batch 0, 0/50 accepted
## batch 1, 0/50 accepted
## batch 2, 1/50 accepted
## batch 3, 1/50 accepted
## batch 4, 1/50 accepted
## batch 5, 1/50 accepted
## batch 6, 1/50 accepted
## batch 7, 3/50 accepted
## batch 8, 3/50 accepted
## batch 9, 3/50 accepted
## batch 10, 3/50 accepted
## batch 11, 3/50 accepted
## batch 12, 4/50 accepted
## batch 13, 6/50 accepted
## batch 14, 6/50 accepted
## batch 15, 6/50 accepted
## batch 16, 6/50 accepted
## batch 17, 6/50 accepted
## batch 18, 6/50 accepted
## batch 19, 8/50 accepted
## batch 20, 8/50 accepted
## batch 21, 8/50 accepted
## batch 22, 9/50 accepted
## batch 23, 10/50 accepted
## batch 24, 10/50 accepted
## batch 25, 10/50 accepted
## batch 26, 11/50 accepted
## batch 27, 12/50 accepted
## batch 28, 13/50 accepted
## batch 29, 13/50 accepted
## batch 30, 13/50 accepted
## batch 31, 13/50 accepted
## batch 32, 13/50 accepted
## batch 33, 14/50 accepted
## batch 34, 14/50 accepted
## batch 35, 14/50 accepted
## batch 36, 14/50 accepted
## batch 37, 16/50 accepted
## batch 38, 16/50 accepted
## batch 39, 16/50 accepted
## batch 40, 17/50 accepted
## batch 41, 17/50 accepted
## batch 42, 17/50 accepted
## batch 43, 18/50 accepted
## batch 44, 18/50 accepted
## batch 45, 18/50 accepted
## batch 46, 18/50 accepted
## batch 47, 18/50 accepted
## batch 48, 18/50 accepted
## batch 49, 18/50 accepted
## batch 50, 18/50 accepted
## batch 51, 19/50 accepted
## batch 52, 19/50 accepted
## batch 53, 20/50 accepted
## batch 54, 20/50 accepted
## batch 55, 20/50 accepted
## batch 56, 20/50 accepted
## batch 57, 21/50 accepted
## batch 58, 21/50 accepted
## batch 59, 22/50 accepted
## batch 60, 22/50 accepted
## batch 61, 23/50 accepted
## batch 62, 24/50 accepted
## batch 63, 24/50 accepted
## batch 64, 24/50 accepted
## batch 65, 26/50 accepted
## batch 66, 26/50 accepted
## batch 67, 26/50 accepted
## batch 68, 27/50 accepted
## batch 69, 27/50 accepted
## batch 70, 29/50 accepted
## batch 71, 29/50 accepted
## batch 72, 29/50 accepted
## batch 73, 30/50 accepted
## batch 74, 31/50 accepted
## batch 75, 31/50 accepted
## batch 76, 31/50 accepted
## batch 77, 31/50 accepted
## batch 78, 31/50 accepted
## batch 79, 31/50 accepted
## batch 80, 31/50 accepted
## batch 81, 31/50 accepted
## batch 82, 32/50 accepted
## batch 83, 32/50 accepted
## batch 84, 32/50 accepted
## batch 85, 32/50 accepted
## batch 86, 32/50 accepted
## batch 87, 32/50 accepted
## batch 88, 32/50 accepted
## batch 89, 32/50 accepted
## batch 90, 32/50 accepted
## batch 91, 32/50 accepted
## batch 92, 33/50 accepted
## batch 93, 34/50 accepted
## batch 94, 34/50 accepted
## batch 95, 34/50 accepted
## batch 96, 34/50 accepted
## batch 97, 35/50 accepted
## batch 98, 35/50 accepted
## batch 99, 36/50 accepted
## batch 100, 36/50 accepted
## batch 101, 37/50 accepted
## batch 102, 37/50 accepted
## batch 103, 37/50 accepted
## batch 104, 37/50 accepted
## batch 105, 37/50 accepted
## batch 106, 37/50 accepted
## batch 107, 38/50 accepted
## batch 108, 39/50 accepted
## batch 109, 39/50 accepted
## batch 110, 41/50 accepted
## batch 111, 42/50 accepted
## batch 112, 42/50 accepted
## batch 113, 42/50 accepted
## batch 114, 42/50 accepted
## batch 115, 43/50 accepted
## batch 116, 43/50 accepted
## batch 117, 43/50 accepted
## batch 118, 43/50 accepted
## batch 119, 44/50 accepted
## batch 120, 44/50 accepted
## batch 121, 44/50 accepted
## batch 122, 45/50 accepted
## batch 123, 45/50 accepted
## batch 124, 45/50 accepted
## batch 125, 45/50 accepted
## batch 126, 46/50 accepted
## batch 127, 47/50 accepted
## batch 128, 49/50 accepted
## Iteration 2 [2435.95, 3123.01] eps: 3124
## batch 0, 0/50 accepted
## batch 1, 0/50 accepted
## batch 2, 0/50 accepted
## batch 3, 0/50 accepted
## batch 4, 1/50 accepted
## batch 5, 1/50 accepted
## batch 6, 1/50 accepted
## batch 7, 2/50 accepted
## batch 8, 2/50 accepted
## batch 9, 2/50 accepted
## batch 10, 2/50 accepted
## batch 11, 2/50 accepted
## batch 12, 2/50 accepted
## batch 13, 2/50 accepted
## batch 14, 2/50 accepted
## batch 15, 3/50 accepted
## batch 16, 4/50 accepted
## batch 17, 4/50 accepted
## batch 18, 4/50 accepted
## batch 19, 4/50 accepted
## batch 20, 4/50 accepted
## batch 21, 5/50 accepted
## batch 22, 6/50 accepted
## batch 23, 6/50 accepted
## batch 24, 6/50 accepted
## batch 25, 6/50 accepted
## batch 26, 6/50 accepted
## batch 27, 6/50 accepted
## batch 28, 6/50 accepted
## batch 29, 7/50 accepted
## batch 30, 8/50 accepted
## batch 31, 8/50 accepted
## batch 32, 9/50 accepted
## batch 33, 9/50 accepted
## batch 34, 10/50 accepted
## batch 35, 10/50 accepted
## batch 36, 11/50 accepted
## batch 37, 11/50 accepted
## batch 38, 11/50 accepted
## batch 39, 13/50 accepted
## batch 40, 13/50 accepted
## batch 41, 14/50 accepted
## batch 42, 16/50 accepted
## batch 43, 16/50 accepted
## batch 44, 17/50 accepted
## batch 45, 17/50 accepted
## batch 46, 17/50 accepted
## batch 47, 17/50 accepted
## batch 48, 18/50 accepted
## batch 49, 18/50 accepted
## batch 50, 18/50 accepted
## batch 51, 18/50 accepted
## batch 52, 18/50 accepted
## batch 53, 18/50 accepted
## batch 54, 18/50 accepted
## batch 55, 18/50 accepted
## batch 56, 18/50 accepted
## batch 57, 18/50 accepted
## batch 58, 19/50 accepted
## batch 59, 19/50 accepted
## batch 60, 19/50 accepted
## batch 61, 21/50 accepted
## batch 62, 23/50 accepted
## batch 63, 23/50 accepted
## batch 64, 23/50 accepted
## batch 65, 23/50 accepted
## batch 66, 24/50 accepted
## batch 67, 24/50 accepted
## batch 68, 24/50 accepted
## batch 69, 24/50 accepted
## batch 70, 24/50 accepted
## batch 71, 24/50 accepted
## batch 72, 25/50 accepted
## batch 73, 25/50 accepted
## batch 74, 25/50 accepted
## batch 75, 25/50 accepted
## batch 76, 26/50 accepted
## batch 77, 26/50 accepted
## batch 78, 27/50 accepted
## batch 79, 28/50 accepted
## batch 80, 29/50 accepted
## batch 81, 29/50 accepted
## batch 82, 29/50 accepted
## batch 83, 29/50 accepted
## batch 84, 29/50 accepted
## batch 85, 29/50 accepted
## batch 86, 29/50 accepted
## batch 87, 29/50 accepted
## batch 88, 29/50 accepted
## batch 89, 29/50 accepted
## batch 90, 29/50 accepted
## batch 91, 30/50 accepted
## batch 92, 30/50 accepted
## batch 93, 30/50 accepted
## batch 94, 31/50 accepted
## batch 95, 31/50 accepted
## batch 96, 32/50 accepted
## batch 97, 33/50 accepted
## batch 98, 33/50 accepted
## batch 99, 33/50 accepted
## batch 100, 33/50 accepted
## batch 101, 33/50 accepted
## batch 102, 33/50 accepted
## batch 103, 33/50 accepted
## batch 104, 33/50 accepted
## batch 105, 33/50 accepted
## batch 106, 33/50 accepted
## batch 107, 33/50 accepted
## batch 108, 34/50 accepted
## batch 109, 34/50 accepted
## batch 110, 34/50 accepted
## batch 111, 34/50 accepted
## batch 112, 34/50 accepted
## batch 113, 34/50 accepted
## batch 114, 35/50 accepted
## batch 115, 35/50 accepted
## batch 116, 36/50 accepted
## batch 117, 37/50 accepted
## batch 118, 37/50 accepted
## batch 119, 37/50 accepted
## batch 120, 37/50 accepted
## batch 121, 37/50 accepted
## batch 122, 37/50 accepted
## batch 123, 38/50 accepted
## batch 124, 39/50 accepted
## batch 125, 39/50 accepted
## batch 126, 39/50 accepted
## batch 127, 39/50 accepted
## batch 128, 39/50 accepted
## batch 129, 40/50 accepted
## batch 130, 40/50 accepted
## batch 131, 40/50 accepted
## batch 132, 40/50 accepted
## batch 133, 42/50 accepted
## batch 134, 42/50 accepted
## batch 135, 42/50 accepted
## batch 136, 42/50 accepted
## batch 137, 42/50 accepted
## batch 138, 42/50 accepted
## batch 139, 44/50 accepted
## batch 140, 44/50 accepted
## batch 141, 44/50 accepted
## batch 142, 44/50 accepted
## batch 143, 44/50 accepted
## batch 144, 44/50 accepted
## batch 145, 44/50 accepted
## batch 146, 45/50 accepted
## batch 147, 45/50 accepted
## batch 148, 46/50 accepted
## batch 149, 46/50 accepted
## batch 150, 46/50 accepted
## batch 151, 46/50 accepted
## batch 152, 47/50 accepted
## batch 153, 47/50 accepted
## batch 154, 47/50 accepted
## batch 155, 47/50 accepted
## batch 156, 47/50 accepted
## batch 157, 47/50 accepted
## batch 158, 47/50 accepted
## batch 159, 47/50 accepted
## batch 160, 47/50 accepted
## batch 161, 47/50 accepted
## batch 162, 47/50 accepted
## batch 163, 48/50 accepted
## batch 164, 48/50 accepted
## batch 165, 48/50 accepted
## batch 166, 48/50 accepted
## batch 167, 48/50 accepted
## batch 168, 48/50 accepted
## batch 169, 49/50 accepted
## batch 170, 49/50 accepted
## batch 171, 49/50 accepted
## batch 172, 49/50 accepted
## Iteration 3 [2273.82, 3080.63] eps: 3092.76
## batch 0, 0/50 accepted
## batch 1, 1/50 accepted
## batch 2, 3/50 accepted
## batch 3, 3/50 accepted
## batch 4, 3/50 accepted
## batch 5, 4/50 accepted
## batch 6, 4/50 accepted
## batch 7, 4/50 accepted
## batch 8, 4/50 accepted
## batch 9, 5/50 accepted
## batch 10, 5/50 accepted
## batch 11, 6/50 accepted
## batch 12, 6/50 accepted
## batch 13, 6/50 accepted
## batch 14, 6/50 accepted
## batch 15, 7/50 accepted
## batch 16, 7/50 accepted
## batch 17, 7/50 accepted
## batch 18, 7/50 accepted
## batch 19, 7/50 accepted
## batch 20, 8/50 accepted
## batch 21, 9/50 accepted
## batch 22, 9/50 accepted
## batch 23, 10/50 accepted
## batch 24, 11/50 accepted
## batch 25, 11/50 accepted
## batch 26, 11/50 accepted
## batch 27, 11/50 accepted
## batch 28, 11/50 accepted
## batch 29, 11/50 accepted
## batch 30, 11/50 accepted
## batch 31, 12/50 accepted
## batch 32, 12/50 accepted
## batch 33, 13/50 accepted
## batch 34, 13/50 accepted
## batch 35, 15/50 accepted
## batch 36, 15/50 accepted
## batch 37, 15/50 accepted
## batch 38, 15/50 accepted
## batch 39, 15/50 accepted
## batch 40, 15/50 accepted
## batch 41, 15/50 accepted
## batch 42, 15/50 accepted
## batch 43, 15/50 accepted
## batch 44, 15/50 accepted
## batch 45, 15/50 accepted
## batch 46, 16/50 accepted
## batch 47, 16/50 accepted
## batch 48, 18/50 accepted
## batch 49, 18/50 accepted
## batch 50, 18/50 accepted
## batch 51, 19/50 accepted
## batch 52, 20/50 accepted
## batch 53, 20/50 accepted
## batch 54, 21/50 accepted
## batch 55, 21/50 accepted
## batch 56, 21/50 accepted
## batch 57, 21/50 accepted
## batch 58, 21/50 accepted
## batch 59, 21/50 accepted
## batch 60, 21/50 accepted
## batch 61, 21/50 accepted
## batch 62, 21/50 accepted
## batch 63, 21/50 accepted
## batch 64, 22/50 accepted
## batch 65, 22/50 accepted
## batch 66, 22/50 accepted
## batch 67, 22/50 accepted
## batch 68, 22/50 accepted
## batch 69, 22/50 accepted
## batch 70, 22/50 accepted
## batch 71, 22/50 accepted
## batch 72, 22/50 accepted
## batch 73, 22/50 accepted
## batch 74, 23/50 accepted
## batch 75, 24/50 accepted
## batch 76, 24/50 accepted
## batch 77, 24/50 accepted
## batch 78, 24/50 accepted
## batch 79, 25/50 accepted
## batch 80, 26/50 accepted
## batch 81, 26/50 accepted
## batch 82, 26/50 accepted
## batch 83, 26/50 accepted
## batch 84, 27/50 accepted
## batch 85, 27/50 accepted
## batch 86, 27/50 accepted
## batch 87, 27/50 accepted
## batch 88, 27/50 accepted
## batch 89, 29/50 accepted
## batch 90, 29/50 accepted
## batch 91, 29/50 accepted
## batch 92, 29/50 accepted
## batch 93, 29/50 accepted
## batch 94, 29/50 accepted
## batch 95, 30/50 accepted
## batch 96, 30/50 accepted
## batch 97, 30/50 accepted
## batch 98, 31/50 accepted
## batch 99, 31/50 accepted
## batch 100, 31/50 accepted
## batch 101, 31/50 accepted
## batch 102, 33/50 accepted
## batch 103, 33/50 accepted
## batch 104, 33/50 accepted
## batch 105, 33/50 accepted
## batch 106, 33/50 accepted
## batch 107, 33/50 accepted
## batch 108, 33/50 accepted
## batch 109, 34/50 accepted
## batch 110, 35/50 accepted
## batch 111, 35/50 accepted
## batch 112, 35/50 accepted
## batch 113, 35/50 accepted
## batch 114, 35/50 accepted
## batch 115, 36/50 accepted
## batch 116, 36/50 accepted
## batch 117, 36/50 accepted
## batch 118, 36/50 accepted
## batch 119, 36/50 accepted
## batch 120, 36/50 accepted
## batch 121, 36/50 accepted
## batch 122, 36/50 accepted
## batch 123, 37/50 accepted
## batch 124, 37/50 accepted
## batch 125, 37/50 accepted
## batch 126, 37/50 accepted
## batch 127, 38/50 accepted
## batch 128, 38/50 accepted
## batch 129, 38/50 accepted
## batch 130, 38/50 accepted
## batch 131, 38/50 accepted
## batch 132, 38/50 accepted
## batch 133, 39/50 accepted
## batch 134, 39/50 accepted
## batch 135, 40/50 accepted
## batch 136, 40/50 accepted
## batch 137, 40/50 accepted
## batch 138, 41/50 accepted
## batch 139, 41/50 accepted
## batch 140, 41/50 accepted
## batch 141, 42/50 accepted
## batch 142, 42/50 accepted
## batch 143, 42/50 accepted
## batch 144, 43/50 accepted
## batch 145, 43/50 accepted
## batch 146, 43/50 accepted
## batch 147, 43/50 accepted
## batch 148, 43/50 accepted
## batch 149, 44/50 accepted
## batch 150, 44/50 accepted
## batch 151, 44/50 accepted
## batch 152, 44/50 accepted
## batch 153, 45/50 accepted
## batch 154, 45/50 accepted
## batch 155, 46/50 accepted
## batch 156, 46/50 accepted
## batch 157, 46/50 accepted
## batch 158, 47/50 accepted
## batch 159, 47/50 accepted
## batch 160, 47/50 accepted
## batch 161, 48/50 accepted
## batch 162, 49/50 accepted
## batch 163, 49/50 accepted
## batch 164, 49/50 accepted
## batch 165, 49/50 accepted
## Iteration 4 [2466.9, 3061.57] eps: 3061.83
## batch 0, 0/50 accepted
## batch 1, 0/50 accepted
## batch 2, 0/50 accepted
## batch 3, 0/50 accepted
## batch 4, 0/50 accepted
## batch 5, 1/50 accepted
## batch 6, 2/50 accepted
## batch 7, 2/50 accepted
## batch 8, 2/50 accepted
## batch 9, 2/50 accepted
## batch 10, 3/50 accepted
## batch 11, 3/50 accepted
## batch 12, 3/50 accepted
## batch 13, 3/50 accepted
## batch 14, 3/50 accepted
## batch 15, 3/50 accepted
## batch 16, 3/50 accepted
## batch 17, 3/50 accepted
## batch 18, 3/50 accepted
## batch 19, 3/50 accepted
## batch 20, 3/50 accepted
## batch 21, 3/50 accepted
## batch 22, 3/50 accepted
## batch 23, 3/50 accepted
## batch 24, 3/50 accepted
## batch 25, 3/50 accepted
## batch 26, 3/50 accepted
## batch 27, 4/50 accepted
## batch 28, 4/50 accepted
## batch 29, 5/50 accepted
## batch 30, 5/50 accepted
## batch 31, 5/50 accepted
## batch 32, 5/50 accepted
## batch 33, 5/50 accepted
## batch 34, 5/50 accepted
## batch 35, 5/50 accepted
## batch 36, 6/50 accepted
## batch 37, 6/50 accepted
## batch 38, 7/50 accepted
## batch 39, 7/50 accepted
## batch 40, 7/50 accepted
## batch 41, 7/50 accepted
## batch 42, 7/50 accepted
## batch 43, 7/50 accepted
## batch 44, 7/50 accepted
## batch 45, 8/50 accepted
## batch 46, 8/50 accepted
## batch 47, 8/50 accepted
## batch 48, 8/50 accepted
## batch 49, 8/50 accepted
## batch 50, 8/50 accepted
## batch 51, 8/50 accepted
## batch 52, 8/50 accepted
## batch 53, 10/50 accepted
## batch 54, 11/50 accepted
## batch 55, 11/50 accepted
## batch 56, 11/50 accepted
## batch 57, 12/50 accepted
## batch 58, 12/50 accepted
## batch 59, 13/50 accepted
## batch 60, 13/50 accepted
## batch 61, 13/50 accepted
## batch 62, 13/50 accepted
## batch 63, 13/50 accepted
## batch 64, 13/50 accepted
## batch 65, 13/50 accepted
## batch 66, 14/50 accepted
## batch 67, 14/50 accepted
## batch 68, 14/50 accepted
## batch 69, 14/50 accepted
## batch 70, 14/50 accepted
## batch 71, 14/50 accepted
## batch 72, 14/50 accepted
## batch 73, 14/50 accepted
## batch 74, 14/50 accepted
## batch 75, 14/50 accepted
## batch 76, 14/50 accepted
## batch 77, 14/50 accepted
## batch 78, 14/50 accepted
## batch 79, 14/50 accepted
## batch 80, 15/50 accepted
## batch 81, 15/50 accepted
## batch 82, 16/50 accepted
## batch 83, 16/50 accepted
## batch 84, 16/50 accepted
## batch 85, 16/50 accepted
## batch 86, 17/50 accepted
## batch 87, 17/50 accepted
## batch 88, 17/50 accepted
## batch 89, 17/50 accepted
## batch 90, 17/50 accepted
## batch 91, 17/50 accepted
## batch 92, 17/50 accepted
## batch 93, 17/50 accepted
## batch 94, 17/50 accepted
## batch 95, 17/50 accepted
## batch 96, 17/50 accepted
## batch 97, 17/50 accepted
## batch 98, 17/50 accepted
## batch 99, 17/50 accepted
## batch 100, 17/50 accepted
## batch 101, 17/50 accepted
## batch 102, 17/50 accepted
## batch 103, 17/50 accepted
## batch 104, 18/50 accepted
## batch 105, 18/50 accepted
## batch 106, 18/50 accepted
## batch 107, 19/50 accepted
## batch 108, 19/50 accepted
## batch 109, 19/50 accepted
## batch 110, 19/50 accepted
## batch 111, 19/50 accepted
## batch 112, 19/50 accepted
## batch 113, 19/50 accepted
## batch 114, 19/50 accepted
## batch 115, 19/50 accepted
## batch 116, 20/50 accepted
## batch 117, 20/50 accepted
## batch 118, 20/50 accepted
## batch 119, 20/50 accepted
## batch 120, 20/50 accepted
## batch 121, 20/50 accepted
## batch 122, 21/50 accepted
## batch 123, 21/50 accepted
## batch 124, 21/50 accepted
## batch 125, 21/50 accepted
## batch 126, 22/50 accepted
## batch 127, 22/50 accepted
## batch 128, 22/50 accepted
## batch 129, 22/50 accepted
## batch 130, 22/50 accepted
## batch 131, 22/50 accepted
## batch 132, 22/50 accepted
## batch 133, 22/50 accepted
## batch 134, 22/50 accepted
## batch 135, 22/50 accepted
## batch 136, 23/50 accepted
## batch 137, 23/50 accepted
## batch 138, 23/50 accepted
## batch 139, 23/50 accepted
## batch 140, 23/50 accepted
## batch 141, 23/50 accepted
## batch 142, 23/50 accepted
## batch 143, 23/50 accepted
## batch 144, 23/50 accepted
## batch 145, 23/50 accepted
## batch 146, 23/50 accepted
## batch 147, 24/50 accepted
## batch 148, 24/50 accepted
## batch 149, 24/50 accepted
## batch 150, 24/50 accepted
## batch 151, 24/50 accepted
## batch 152, 24/50 accepted
## batch 153, 25/50 accepted
## batch 154, 26/50 accepted
## batch 155, 26/50 accepted
## batch 156, 26/50 accepted
## batch 157, 27/50 accepted
## batch 158, 27/50 accepted
## batch 159, 27/50 accepted
## batch 160, 27/50 accepted
## batch 161, 27/50 accepted
## batch 162, 28/50 accepted
## batch 163, 28/50 accepted
## batch 164, 28/50 accepted
## batch 165, 28/50 accepted
## batch 166, 28/50 accepted
## batch 167, 28/50 accepted
## batch 168, 28/50 accepted
## batch 169, 28/50 accepted
## batch 170, 28/50 accepted
## batch 171, 29/50 accepted
## batch 172, 29/50 accepted
## batch 173, 29/50 accepted
## batch 174, 29/50 accepted
## batch 175, 29/50 accepted
## batch 176, 29/50 accepted
## batch 177, 29/50 accepted
## batch 178, 29/50 accepted
## batch 179, 29/50 accepted
## batch 180, 29/50 accepted
## batch 181, 29/50 accepted
## batch 182, 29/50 accepted
## batch 183, 29/50 accepted
## batch 184, 30/50 accepted
## batch 185, 30/50 accepted
## batch 186, 30/50 accepted
## batch 187, 30/50 accepted
## batch 188, 31/50 accepted
## batch 189, 31/50 accepted
## batch 190, 31/50 accepted
## batch 191, 32/50 accepted
## batch 192, 32/50 accepted
## batch 193, 32/50 accepted
## batch 194, 33/50 accepted
## batch 195, 33/50 accepted
## batch 196, 35/50 accepted
## batch 197, 35/50 accepted
## batch 198, 35/50 accepted
## batch 199, 35/50 accepted
##
## Maximum batches exceeded: 36/50 acceptances in 200 batches of max 200
## Returning last params
## Simulation completeprint(t1)## user system elapsed
## 1580.882 9.890 122.312print(t2)## user system elapsed
## 1932.335 11.781 149.215Parameter Summaries
summary(weser.model1)## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Data Model Parameters: 0
## 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 -0.292 1.468 -2.798 2.225
## Beta_SE_2 -0.621 2.786 -7.024 4.365
## Beta_SE_3 -0.694 1.722 -3.969 2.634
## Beta_SE_4 -0.529 1.762 -3.919 2.205
## Beta_SE_5 0.414 1.386 -2.346 3.066
## Beta_SE_6 0.637 1.308 -1.945 3.730
## Beta_SE_7 0.132 0.922 -1.757 1.823
## Beta_SE_8 0.640 2.427 -3.816 5.423
## rho_1 0.137 0.087 0.017 0.345
## gamma_EI 0.589 0.721 0.057 2.644
## gamma_IR 1.214 0.843 0.109 2.788
## S0_1 75986.000 0.000 75986.000 75986.000
## S0_2 51445.000 0.000 51445.000 51445.000
## S0_3 158340.000 0.000 158340.000 158340.000
## S0_4 165517.000 0.000 165517.000 165517.000
## S0_5 84586.000 0.000 84586.000 84586.000
## S0_6 114524.000 0.000 114524.000 114524.000
## S0_7 189648.000 0.000 189648.000 189648.000
## S0_8 153283.000 0.000 153283.000 153283.000
## S0_9 307734.000 0.000 307734.000 307734.000
## S0_10 101657.000 0.000 101657.000 101657.000
## S0_11 132975.000 0.000 132975.000 132975.000
## S0_12 164536.000 0.000 164536.000 164536.000
## S0_13 124564.000 0.000 124564.000 124564.000
## S0_14 358041.000 0.000 358041.000 358041.000
## S0_15 130471.000 0.000 130471.000 130471.000
## S0_16 94242.000 0.000 94242.000 94242.000
## S0_17 57672.000 0.000 57672.000 57672.000
## E0_1 0.000 0.000 0.000 0.000
## E0_2 0.000 0.000 0.000 0.000
## E0_3 0.000 0.000 0.000 0.000
## E0_4 0.000 0.000 0.000 0.000
## E0_5 0.000 0.000 0.000 0.000
## E0_6 0.000 0.000 0.000 0.000
## E0_7 2.000 0.000 2.000 2.000
## E0_8 0.000 0.000 0.000 0.000
## E0_9 0.000 0.000 0.000 0.000
## E0_10 0.000 0.000 0.000 0.000
## E0_11 0.000 0.000 0.000 0.000
## E0_12 2.000 0.000 2.000 2.000
## E0_13 0.000 0.000 0.000 0.000
## E0_14 0.000 0.000 0.000 0.000
## E0_15 0.000 0.000 0.000 0.000
## E0_16 0.000 0.000 0.000 0.000
## E0_17 0.000 0.000 0.000 0.000
## I0_1 0.000 0.000 0.000 0.000
## I0_2 0.000 0.000 0.000 0.000
## I0_3 0.000 0.000 0.000 0.000
## I0_4 0.000 0.000 0.000 0.000
## I0_5 0.000 0.000 0.000 0.000
## I0_6 0.000 0.000 0.000 0.000
## I0_7 2.000 0.000 2.000 2.000
## I0_8 0.000 0.000 0.000 0.000
## I0_9 0.000 0.000 0.000 0.000
## I0_10 0.000 0.000 0.000 0.000
## I0_11 0.000 0.000 0.000 0.000
## I0_12 2.000 0.000 2.000 2.000
## I0_13 0.000 0.000 0.000 0.000
## I0_14 0.000 0.000 0.000 0.000
## I0_15 0.000 0.000 0.000 0.000
## I0_16 0.000 0.000 0.000 0.000
## I0_17 0.000 0.000 0.000 0.000
## R0_1 0.000 0.000 0.000 0.000
## R0_2 0.000 0.000 0.000 0.000
## R0_3 0.000 0.000 0.000 0.000
## R0_4 0.000 0.000 0.000 0.000
## R0_5 0.000 0.000 0.000 0.000
## R0_6 0.000 0.000 0.000 0.000
## R0_7 0.000 0.000 0.000 0.000
## R0_8 0.000 0.000 0.000 0.000
## R0_9 0.000 0.000 0.000 0.000
## R0_10 0.000 0.000 0.000 0.000
## R0_11 0.000 0.000 0.000 0.000
## R0_12 0.000 0.000 0.000 0.000
## R0_13 0.000 0.000 0.000 0.000
## R0_14 0.000 0.000 0.000 0.000
## R0_15 0.000 0.000 0.000 0.000
## R0_16 0.000 0.000 0.000 0.000
## R0_17 0.000 0.000 0.000 0.000summary(weser.model2)## Summary: SEIR Model
##
## Locations: 17
## Time Points: 96
## Data Model Parameters: 0
## 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.167 2.376 -6.749 2.162
## Beta_SE_2 0.280 7.582 -14.423 11.979
## Beta_SE_3 0.481 2.718 -4.113 6.083
## Beta_SE_4 -2.283 1.728 -5.001 1.070
## Beta_SE_5 3.244 2.029 -0.511 6.775
## Beta_SE_6 1.204 1.151 -1.215 3.425
## Beta_SE_7 -0.603 1.748 -4.344 2.432
## Beta_SE_8 -0.743 2.742 -5.501 4.499
## rho_1 0.069 0.030 0.012 0.130
## gamma_EI 1.587 1.328 0.061 4.454
## gamma_IR 1.878 1.300 0.105 4.337
## S0_1 75986.000 0.000 75986.000 75986.000
## S0_2 51445.000 0.000 51445.000 51445.000
## S0_3 158340.000 0.000 158340.000 158340.000
## S0_4 165517.000 0.000 165517.000 165517.000
## S0_5 84586.000 0.000 84586.000 84586.000
## S0_6 114524.000 0.000 114524.000 114524.000
## S0_7 189648.000 0.000 189648.000 189648.000
## S0_8 153283.000 0.000 153283.000 153283.000
## S0_9 307734.000 0.000 307734.000 307734.000
## S0_10 101657.000 0.000 101657.000 101657.000
## S0_11 132975.000 0.000 132975.000 132975.000
## S0_12 164536.000 0.000 164536.000 164536.000
## S0_13 124564.000 0.000 124564.000 124564.000
## S0_14 358041.000 0.000 358041.000 358041.000
## S0_15 130471.000 0.000 130471.000 130471.000
## S0_16 94242.000 0.000 94242.000 94242.000
## S0_17 57672.000 0.000 57672.000 57672.000
## E0_1 0.000 0.000 0.000 0.000
## E0_2 0.000 0.000 0.000 0.000
## E0_3 0.000 0.000 0.000 0.000
## E0_4 0.000 0.000 0.000 0.000
## E0_5 0.000 0.000 0.000 0.000
## E0_6 0.000 0.000 0.000 0.000
## E0_7 2.000 0.000 2.000 2.000
## E0_8 0.000 0.000 0.000 0.000
## E0_9 0.000 0.000 0.000 0.000
## E0_10 0.000 0.000 0.000 0.000
## E0_11 0.000 0.000 0.000 0.000
## E0_12 2.000 0.000 2.000 2.000
## E0_13 0.000 0.000 0.000 0.000
## E0_14 0.000 0.000 0.000 0.000
## E0_15 0.000 0.000 0.000 0.000
## E0_16 0.000 0.000 0.000 0.000
## E0_17 0.000 0.000 0.000 0.000
## I0_1 0.000 0.000 0.000 0.000
## I0_2 0.000 0.000 0.000 0.000
## I0_3 0.000 0.000 0.000 0.000
## I0_4 0.000 0.000 0.000 0.000
## I0_5 0.000 0.000 0.000 0.000
## I0_6 0.000 0.000 0.000 0.000
## I0_7 2.000 0.000 2.000 2.000
## I0_8 0.000 0.000 0.000 0.000
## I0_9 0.000 0.000 0.000 0.000
## I0_10 0.000 0.000 0.000 0.000
## I0_11 0.000 0.000 0.000 0.000
## I0_12 2.000 0.000 2.000 2.000
## I0_13 0.000 0.000 0.000 0.000
## I0_14 0.000 0.000 0.000 0.000
## I0_15 0.000 0.000 0.000 0.000
## I0_16 0.000 0.000 0.000 0.000
## I0_17 0.000 0.000 0.000 0.000
## R0_1 0.000 0.000 0.000 0.000
## R0_2 0.000 0.000 0.000 0.000
## R0_3 0.000 0.000 0.000 0.000
## R0_4 0.000 0.000 0.000 0.000
## R0_5 0.000 0.000 0.000 0.000
## R0_6 0.000 0.000 0.000 0.000
## R0_7 0.000 0.000 0.000 0.000
## R0_8 0.000 0.000 0.000 0.000
## R0_9 0.000 0.000 0.000 0.000
## R0_10 0.000 0.000 0.000 0.000
## R0_11 0.000 0.000 0.000 0.000
## R0_12 0.000 0.000 0.000 0.000
## R0_13 0.000 0.000 0.000 0.000
## R0_14 0.000 0.000 0.000 0.000
## R0_15 0.000 0.000 0.000 0.000
## R0_16 0.000 0.000 0.000 0.000
## R0_17 0.000 0.000 0.000 0.000Graphical Summary Function
This function will plot cases each of the spatial locations for a model, add a black line for the mean posterior predictive number of cases, and add blue dashed lines for the 99% credible interval for the number of cases.
makePlots <- function(modelObj, nm){
sims <- epidemic.simulations(modelObj, replicates = 50)
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(Isub[,x])),
main = paste(model, ": Observed and Posterior Predictive Simulation.\n location ",
colnames(neighbourhood)[x], sep = ""))
lines(Ism[,x], col = rgb(0,0,0,0.8), lwd = 2)
lines(Islb[,x], col = rgb(0,0,0.5,0.8), lwd = 1, lty = 2)
lines(Isub[,x], col = rgb(0,0,0.5,0.8), lwd = 1, lty = 2)
#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)
}
}Graphical Illustration of Model 1:
makePlots(weser.model1, "Dist")
Graphical Illustration of Model 2:
makePlots(weser.model2, "CAR")
Bayes Factors Comparing the Models (row v. column):
comps <- compareModels(modelList = list(weser.model1, weser.model2), n_samples = 100)## [1] "Assuming equal prior probabilities."rownames(comps) <- colnames(comps) <- c("Distance", "CAR")
print(comps)## Distance CAR
## Distance 1.000000 0.8947368
## CAR 1.117647 1.0000000It looks like the preferred model is the CAR model.