Package 'SimInf'

Description

Approximate Bayesian computation

Usage

abc(
  model,
  priors = NULL,
  n_particles = NULL,
  n_init = NULL,
  distance = NULL,
  tolerance = NULL,
  data = NULL,
  verbose = getOption("verbose", FALSE),
  post_gen = NULL,
  init_model = NULL
)

## S4 method for signature 'SimInf_model'
abc(
  model,
  priors = NULL,
  n_particles = NULL,
  n_init = NULL,
  distance = NULL,
  tolerance = NULL,
  data = NULL,
  verbose = getOption("verbose", FALSE),
  post_gen = NULL,
  init_model = NULL
)

Arguments

Value

A SimInf_abc object.

References

T. Toni, D. Welch, N. Strelkowa, A. Ipsen, and M. P. H. Stumpf. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface 6, 187–202, 2009. doi:10.1098/rsif.2008.0172

U. Simola, J. Cisewski-Kehe, M. U. Gutmann, J. Corander. Adaptive Approximate Bayesian Computation Tolerance Selection. Bayesian Analysis, 16(2), 397–423, 2021. doi: 10.1214/20-BA1211

See Also

continue_abc.

Examples

## Not run: 
## Let us consider an SIR model in a closed population with N = 100
## individuals of whom one is initially infectious and the rest are
## susceptible. First, generate one realisation (with a specified
## seed) from the model with known parameters \code{beta = 0.16} and
## \code{gamma = 0.077}. Then, use \code{abc} to infer the (known)
## parameters from the simulated data.
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the SIR model and plot the number of infectious.
set.seed(22)
infectious <- trajectory(run(model), "I")$I
plot(infectious, type = "s")

## The distance function to accept or reject a proposal. Each node
## in the simulated trajectory (contained in the 'result' object)
## represents one proposal.
distance <- function(result, ...) {
    ## Extract the time-series of infectious in each node as a
    ## data.frame.
    sim <- trajectory(result, "I")

    ## Split the 'sim' data.frame by node and calculate the sum of the
    ## squared distance at each time-point for each node.
    dist <- tapply(sim$I, sim$node, function(sim_infectious) {
        sum((infectious - sim_infectious)^2)
    })

    ## Return the distance for each node. Each proposal will be
    ## accepted or rejected depending on if the distance is less than
    ## the tolerance for the current generation.
    dist
}

## Fit the model parameters using ABC-SMC and adaptive tolerance
## selection. The priors for the parameters are specified using a
## formula notation. Here we use a uniform distribtion for each
## parameter with lower bound = 0 and upper bound = 1. Note that we
## use a low number particles here to keep the run-time of the example
## short. In practice you would want to use many more to ensure better
## approximations.
fit <- abc(model = model,
           priors = c(beta ~ uniform(0, 1), gamma ~ uniform(0, 1)),
           n_particles = 100,
           n_init = 1000,
           distance = distance,
           verbose = TRUE)

## Print a brief summary.
fit

## Display the ABC posterior distribution.
plot(fit)

## End(Not run)

Description

Coerce to data frame

Usage

## S3 method for class 'SimInf_abc'
as.data.frame(x, ...)

Arguments

Description

Coerce events to a data frame

Usage

## S3 method for class 'SimInf_events'
as.data.frame(x, ...)

Arguments

Description

Coerce to data frame

Usage

## S3 method for class 'SimInf_indiv_events'
as.data.frame(x, ...)

Arguments

Description

Coerce to data frame

Usage

## S3 method for class 'SimInf_pmcmc'
as.data.frame(x, ...)

Arguments

Description

Produce box-and-whisker plot(s) of the number of individuals in each model compartment.

Usage

## S4 method for signature 'SimInf_model'
boxplot(x, compartments = NULL, index = NULL, ...)

Arguments

Examples

## Create an 'SIR' model with 10 nodes and initialise
## it with 99 susceptible individuals and one infected
## individual. Let the model run over 100 days.
model <- SIR(u0 = data.frame(S = rep(99, 10),
                             I = rep(1, 10),
                             R = rep(0, 10)),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the model and save the result.
result <- run(model)

## Create a boxplot that includes all compartments in all nodes.
boxplot(result)

## Create a boxplot that includes the S and I compartments in
## nodes 1 and 2.
boxplot(result, ~S+I, 1:2)

Description

Extract the C code from a SimInf_model object

Usage

C_code(model)

Arguments

Value

Character vector with C code for the model.

Examples

## Use the model parser to create a 'SimInf_model' object that
## expresses an SIR model, where 'b' is the transmission rate and
## 'g' is the recovery rate.
model <- mparse(transitions = c("S -> b*S*I/(S+I+R) -> I", "I -> g*I -> R"),
                compartments = c("S", "I", "R"),
                gdata = c(b = 0.16, g = 0.077),
                u0 = data.frame(S = 99, I = 1, R = 0),
                tspan = 1:10)

## View the C code.
C_code(model)

Description

Run more generations of ABC SMC

Usage

continue_abc(
  object,
  tolerance = NULL,
  data = NULL,
  verbose = getOption("verbose", FALSE),
  post_gen = NULL
)

## S4 method for signature 'SimInf_abc'
continue_abc(
  object,
  tolerance = NULL,
  data = NULL,
  verbose = getOption("verbose", FALSE),
  post_gen = NULL
)

Arguments

Value

A SimInf_abc object.

Description

Run more iterations of PMCMC

Usage

continue_pmcmc(
  object,
  obs_process,
  n_iterations,
  init_model = NULL,
  post_particle = NULL,
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'SimInf_pmcmc'
continue_pmcmc(
  object,
  obs_process,
  n_iterations,
  init_model = NULL,
  post_particle = NULL,
  verbose = getOption("verbose", FALSE)
)

Arguments

Description

Calculate the euclidian distances beween coordinates for all coordinates within the cutoff.

Usage

distance_matrix(x, y, cutoff, min_dist = NULL, na_fail = TRUE)

Arguments

Value

dgCMatrix

Examples

## Generate a grid 10 x 10 and place one node in each cell
## separated by 100m.
nodes <- expand.grid(x = (0:9) * 100, y = (0:9) * 100)
plot(y ~ x, nodes)

## Define the cutoff to only include neighbors within 300m.
d <- distance_matrix(x = nodes$x, y = nodes$y, cutoff = 300)

## View the first 10 rows and columns in the distance matrix
d[1:10, 1:10]

Description

A utility function to facilitate preparing edge properties for ldata in a model.

Usage

edge_properties_to_matrix(edges, n_nodes)

Arguments

Details

The edge properties will be converted to a matrix where each row in edges will become a sequence of (index, value_1, value_2, ..., value_n) where 'index' is the zero-based index of the from node. The reason for a zero-based index is to facilitate it's usage in C code. The sequence will be added to the 'to' column in the matrix. There will always be at least one stop value=-1 in each column. All other values in the matrix will be set to NaN. See ‘Examples’.

Value

a numeric matrix with the number of rows equal to max(table(edges$to)) * (ncol(edges) - 1) + 1 and the number of columns equal to n_nodes.

Examples

## Let us consider the following edge properties.
edges <- data.frame(
    from  = c(  2,    3,     4,  1,   4,    5,   1,   3,   1,   3),
    to    = c(  1,    1,     1,  2,   3,    3,   4,   4,   5,   5),
    rate  = c(0.2, 0.01,  0.79,  1, 0.2, 0.05, 0.2, 0.8, 0.2, 0.8),
    count = c(  5,    5,     5, 50,  10,   10,   5,   5,   5,   5))

## Converting the edge properties into a matrix
edge_properties_to_matrix(edges, 6)

## Gives the following output. The first column contains first the
## properties for the edge from = 2 --> to = 1, where the first
## row is the zero-based index of from, i.e., 1. The second row
## contains the rate=0.2 and the third row count=5. On the fourth
## row starts the next sequence with the values in the second row
## in the edges data.frame. The stop value in the first column is
## on row 10. As can be seen in column 6, there are no edge
## properties for node=6.
##        [,1] [,2]  [,3] [,4] [,5] [,6]
##  [1,]  1.00    0  3.00  0.0  0.0   -1
##  [2,]  0.20    1  0.20  0.2  0.2  NaN
##  [3,]  5.00   50 10.00  5.0  5.0  NaN
##  [4,]  2.00   -1  4.00  2.0  2.0  NaN
##  [5,]  0.01  NaN  0.05  0.8  0.8  NaN
##  [6,]  5.00  NaN 10.00  5.0  5.0  NaN
##  [7,]  3.00  NaN -1.00 -1.0 -1.0  NaN
##  [8,]  0.79  NaN   NaN  NaN  NaN  NaN
##  [9,]  5.00  NaN   NaN  NaN  NaN  NaN
## [10,] -1.00  NaN   NaN  NaN  NaN  NaN

Description

Extract the scheduled events from a SimInf_model object.

Usage

events(object, ...)

## S4 method for signature 'SimInf_model'
events(object, ...)

Arguments

Value

SimInf_events object.

Examples

## Create an SIR model that includes scheduled events.
model <- SIR(u0     = u0_SIR(),
             tspan  = 1:(4 * 365),
             events = events_SIR(),
             beta   = 0.16,
             gamma  = 0.077)

## Extract the scheduled events from the model and display summary
summary(events(model))

## Extract the scheduled events from the model and plot them
plot(events(model))

Description

Example data to initialize scheduled events for a population of 1600 nodes and demonstrate the SEIR model.

Usage

events_SEIR()

Details

Example data to initialize scheduled events (see SimInf_events) for a population of 1600 nodes and demonstrate the SEIR model. The dataset contains 466692 events for 1600 nodes distributed over 4 * 365 days. The events are divided into three types: ‘Exit’ events remove individuals from the population (n = 182535), ‘Enter’ events add individuals to the population (n = 182685), and ‘External transfer’ events move individuals between nodes in the population (n = 101472). The vignette contains a detailed description of how scheduled events operate on a model.

Value

A data.frame

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SEIR' model with 1600 nodes and initialize
## it to run over 4*365 days. Add one infected individual
## to the first node.
u0 <- u0_SEIR()
u0$I[1] <- 1
tspan <- seq(from = 1, to = 4*365, by = 1)
model <- SEIR(u0      = u0,
              tspan   = tspan,
              events  = events_SEIR(),
              beta    = 0.16,
              epsilon = 0.25,
              gamma   = 0.01)

## Display the number of individuals affected by each event type
## per day.
plot(events(model))

## Run the model to generate a single stochastic trajectory.
result <- run(model)
plot(result)

## Summarize the trajectory. The summary includes the number of
## events by event type.
summary(result)

Description

Example data to initialize scheduled events for a population of 1600 nodes and demonstrate the SIR model.

Usage

events_SIR()

Details

Example data to initialize scheduled events (see SimInf_events) for a population of 1600 nodes and demonstrate the SIR model. The dataset contains 466692 events for 1600 nodes distributed over 4 * 365 days. The events are divided into three types: ‘Exit’ events remove individuals from the population (n = 182535), ‘Enter’ events add individuals to the population (n = 182685), and ‘External transfer’ events move individuals between nodes in the population (n = 101472). The vignette contains a detailed description of how scheduled events operate on a model.

Value

A data.frame

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIR' model with 1600 nodes and initialize
## it to run over 4*365 days. Add one infected individual
## to the first node.
u0 <- u0_SIR()
u0$I[1] <- 1
tspan <- seq(from = 1, to = 4*365, by = 1)
model <- SIR(u0     = u0,
             tspan  = tspan,
             events = events_SIR(),
             beta   = 0.16,
             gamma  = 0.01)

## Display the number of individuals affected by each event type
## per day.
plot(events(model))

## Run the model to generate a single stochastic trajectory.
result <- run(model)
plot(result)

## Summarize the trajectory. The summary includes the number of
## events by event type.
summary(result)

Description

Example data to initialize scheduled events for a population of 1600 nodes and demonstrate the SIS model.

Usage

events_SIS()

Details

Example data to initialize scheduled events (see SimInf_events) for a population of 1600 nodes and demonstrate the SIS model. The dataset contains 466692 events for 1600 nodes distributed over 4 * 365 days. The events are divided into three types: ‘Exit’ events remove individuals from the population (n = 182535), ‘Enter’ events add individuals to the population (n = 182685), and ‘External transfer’ events move individuals between nodes in the population (n = 101472). The vignette contains a detailed description of how scheduled events operate on a model.

Value

A data.frame

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIS' model with 1600 nodes and initialize
## it to run over 4*365 days. Add one infected individual
## to the first node.
u0 <- u0_SIS()
u0$I[1] <- 1
tspan <- seq(from = 1, to = 4*365, by = 1)
model <- SIS(u0     = u0,
             tspan  = tspan,
             events = events_SIS(),
             beta   = 0.16,
             gamma  = 0.01)

## Display the number of individuals affected by each event type
## per day.
plot(events(model))

## Run the model to generate a single stochastic trajectory.
result <- run(model)
plot(result)

## Summarize the trajectory. The summary includes the number of
## events by event type.
summary(result)

Description

Example data to initialize scheduled events for a population of 1600 nodes and demonstrate the SISe model.

Usage

events_SISe()

Details

Example data to initialize scheduled events (see SimInf_events) for a population of 1600 nodes and demonstrate the SISe model. The dataset contains 466692 events for 1600 nodes distributed over 4 * 365 days. The events are divided into three types: ‘Exit’ events remove individuals from the population (n = 182535), ‘Enter’ events add individuals to the population (n = 182685), and ‘External transfer’ events move individuals between nodes in the population (n = 101472). The vignette contains a detailed description of how scheduled events operate on a model.

Value

A data.frame

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SISe' model with 1600 nodes and initialize
## it to run over 4*365 days. Add one infected individual
## to the first node.
u0 <- u0_SISe()
u0$I[1] <- 1
tspan <- seq(from = 1, to = 4*365, by = 1)
model <- SISe(u0 = u0, tspan = tspan, events = events_SISe(),
              phi = 0, upsilon = 1.8e-2, gamma = 0.1, alpha = 1,
              beta_t1 = 1.0e-1, beta_t2 = 1.0e-1, beta_t3 = 1.25e-1,
              beta_t4 = 1.25e-1, end_t1 = 91, end_t2 = 182,
              end_t3 = 273, end_t4 = 365, epsilon = 0)

## Display the number of individuals affected by each event type
## per day.
plot(events(model))

## Run the model to generate a single stochastic trajectory.
result <- run(model)

## Summarize the trajectory. The summary includes the number of
## events by event type.
summary(result)

Description

Example data to initialize scheduled events for a population of 1600 nodes and demonstrate the SISe3 model.

Usage

data(events_SISe3)

Format

A data.frame

Details

Example data to initialize scheduled events (see SimInf_events) for a population of 1600 nodes and demonstrate the SISe3 model. The dataset contains 783773 events for 1600 nodes distributed over 4 * 365 days. The events are divided into three types: ‘Exit’ events remove individuals from the population (n = 182535), ‘Enter’ events add individuals to the population (n = 182685), ‘Internal transfer’ events move individuals between compartmens within one node e.g. ageing (n = 317081), and ‘External transfer’ events move individuals between nodes in the population (n = 101472). The vignette contains a detailed description of how scheduled events operate on a model.

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SISe3' model with 1600 nodes and initialize
## it to run over 4*365 days. Add one infected individual
## to the first node.
data("u0_SISe3", package = "SimInf")
data("events_SISe3", package = "SimInf")
u0_SISe3$I_1[1] <- 1
tspan <- seq(from = 1, to = 4*365, by = 1)
model <- SISe3(u0 = u0_SISe3, tspan = tspan, events = events_SISe3,
               phi = rep(0, nrow(u0_SISe3)), upsilon_1 = 1.8e-2,
               upsilon_2 = 1.8e-2, upsilon_3 = 1.8e-2,
               gamma_1 = 0.1, gamma_2 = 0.1, gamma_3 = 0.1,
               alpha = 1, beta_t1 = 1.0e-1, beta_t2 = 1.0e-1,
               beta_t3 = 1.25e-1, beta_t4 = 1.25e-1, end_t1 = 91,
               end_t2 = 182, end_t3 = 273, end_t4 = 365, epsilon = 0)

## Display the number of individuals affected by each event type
## per day.
plot(events(model))

## Run the model to generate a single stochastic trajectory.
result <- run(model)

## Summarize the trajectory. The summary includes the number of
## events by event type.
summary(result)

Description

The global data is a numeric vector that is common to all nodes. The global data vector is passed as an argument to the transition rate functions and the post time step function.

Usage

gdata(model)

## S4 method for signature 'SimInf_model'
gdata(model)

Arguments

Value

a numeric vector

Examples

## Create an SIR model
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:5, beta = 0.16, gamma = 0.077)

## Set 'beta' to a new value
gdata(model, "beta") <- 2

## Extract the global data vector that is common to all nodes
gdata(model)

Description

The global data is a numeric vector that is common to all nodes. The global data vector is passed as an argument to the transition rate functions and the post time step function.

Usage

gdata(model, parameter) <- value

## S4 replacement method for signature 'SimInf_model'
gdata(model, parameter) <- value

Arguments

Value

a SimInf_model object

Examples

## Create an SIR model
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:5, beta = 0.16, gamma = 0.077)

## Set 'beta' to a new value
gdata(model, "beta") <- 2

## Extract the global data vector that is common to all nodes
gdata(model)

Description

Lookup individuals, in which node they are located and their age at a specified time-point.

Usage

get_individuals(x, time = NULL)

## S4 method for signature 'SimInf_indiv_events'
get_individuals(x, time = NULL)

Arguments

Value

a data.frame with the columns id, node, and age.

Description

The number of nodes with inward external transfer events to each node.

Usage

indegree(model)

Arguments

Value

vector with in-degree for each node.

Examples

## Create an 'SIR' model with 1600 nodes and initialize
## it with example data.
model <- SIR(u0 = u0_SIR(), tspan = 1:1460, events = events_SIR(),
             beta   = 0.16, gamma  = 0.077)

## Display indegree for each node in the model.
plot(indegree(model))

Description

In many countries, individual-based livestock data are collected to enable contact tracing during disease outbreaks. However, the livestock databases are not always structured in such a way that relevant information for disease spread simulations is easily retrieved. The aim of this function is to facilitate cleaning livestock event data and prepare it for usage in SimInf.

Usage

individual_events(events)

Arguments

Details

The argument events in individual_events must be a data.frame with the following columns:

• id: an integer or character identifier of the individual.

• event: four event types are supported: exit, enter, internal transfer, and external transfer. When assigning the events, they can either be coded as a numerical value or a character string: exit; 0 or 'exit', enter; 1 or 'enter', internal transfer; 2 or 'intTrans', and external transfer; 3 or 'extTrans'.

• time: an integer, character, or date (of class Date) for when the event occured. If it's a character it must be able to coerce to Date.

• node: an integer or character identifier of the source node.

• dest: an integer or character identifier of the destination node for movement events, and 'dest' will be replaced with NA for non-movement events.

Value

SimInf_indiv_events

See Also

node_events.

Description

The local data is a numeric vector that is specific to a node. The local data vector is passed as an argument to the transition rate functions and the post time step function.

Usage

ldata(model, node)

## S4 method for signature 'SimInf_model'
ldata(model, node)

Arguments

Value

a numeric vector

Examples

## Create an 'SISe' model with 1600 nodes.
model <- SISe(u0 = u0_SISe(), tspan = 1:100, events = events_SISe(),
              phi = 0, upsilon = 1.8e-2, gamma = 0.1, alpha = 1,
              beta_t1 = 1.0e-1, beta_t2 = 1.0e-1, beta_t3 = 1.25e-1,
              beta_t4 = 1.25e-1, end_t1 = c(91, 101), end_t2 = c(182, 185),
              end_t3 = c(273, 275), end_t4 = c(365, 360), epsilon = 0)

## Display local data from the first two nodes.
ldata(model, node = 1)
ldata(model, node = 2)

Description

Length of the MCMC chain

Usage

## S4 method for signature 'SimInf_pmcmc'
length(x)

Arguments

Description

Extract the estimated log likelihood from a SimInf_pfilter object.

Usage

## S4 method for signature 'SimInf_pfilter'
logLik(object)

Arguments

Value

the estimated log likelihood.

Description

Describe your model in a logical way in R. mparse creates a SimInf_model object with your model definition that is ready to run.

Usage

mparse(
  transitions = NULL,
  compartments = NULL,
  ldata = NULL,
  gdata = NULL,
  u0 = NULL,
  v0 = NULL,
  tspan = NULL,
  events = NULL,
  E = NULL,
  N = NULL,
  pts_fun = NULL,
  use_enum = FALSE
)

Arguments

Value

a SimInf_model object

Examples

## Not run: 
## Use the model parser to create a 'SimInf_model' object that
## expresses the SIR model, where 'beta' is the transmission rate
## and 'gamma' is the recovery rate.
model  <- mparse(transitions = c("S -> beta*S*I/N -> I",
                                 "I -> gamma*I -> R",
                                 "N <- S+I+R"),
                 compartments = c("S", "I", "R"),
                 gdata = c(beta = 0.16, gamma = 0.077),
                 u0 = data.frame(S = 100, I = 1, R = 0),
                 tspan = 1:100)

## Run and plot the result
set.seed(22)
result <- run(model)
plot(result)

## End(Not run)

Description

Determine the number of compartments in a model

Usage

n_compartments(model)

## S4 method for signature 'SimInf_model'
n_compartments(model)

Arguments

Value

the number of compartments in the model.

Examples

## Create an 'SIR' model with 100 nodes, with 99 susceptible,
## 1 infected and 0 recovered in each node.
u0 <- data.frame(S = rep(99, 100), I = rep(1, 100), R = rep(0, 100))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)

## Display the number of compartments in the model.
n_compartments(model)

Description

Determine the number of generations

Usage

n_generations(object)

## S4 method for signature 'SimInf_abc'
n_generations(object)

Arguments

Value

an integer with the number of generations.

Description

Determine the number of nodes in a model

Usage

n_nodes(model)

## S4 method for signature 'SimInf_model'
n_nodes(model)

## S4 method for signature 'SimInf_pfilter'
n_nodes(model)

## S4 method for signature 'SimInf_pmcmc'
n_nodes(model)

Arguments

Value

the number of nodes in the model.

Examples

## Create an 'SIR' model with 100 nodes, with 99 susceptible,
## 1 infected and 0 recovered in each node.
u0 <- data.frame(S = rep(99, 100), I = rep(1, 100), R = rep(0, 100))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)

## Display the number of nodes in the model.
n_nodes(model)

Description

Determine the number of replicates in a model

Usage

n_replicates(model)

## S4 method for signature 'SimInf_model'
n_replicates(model)

Arguments

Value

the number of replicates in the model.

Examples

## Create an 'SIR' model with 100 nodes, with 99 susceptible,
## 1 infected and 0 recovered in each node.
u0 <- data.frame(S = rep(99, 100), I = rep(1, 100), R = rep(0, 100))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)

## Display the number of replicates in the model.
n_replicates(model)

Description

In many countries, individual-based livestock data are collected to enable contact tracing during disease outbreaks. However, the livestock databases are not always structured in such a way that relevant information for disease spread simulations is easily retrieved. The aim of this function is to facilitate cleaning livestock event data and prepare it for usage in SimInf.

Usage

node_events(x, time = NULL, target = NULL, age = NULL)

## S4 method for signature 'SimInf_indiv_events'
node_events(x, time = NULL, target = NULL, age = NULL)

Arguments

Details

The individual-based events will be aggregated on node-level. The select value is determined by the event type and age category. If there is only one age category, i.e., age=NULL, then select=1 for the enter events, and select=2 for all other events. If there are two age categories, then select=1 for the enter events in the first age category, and select=2 for the enter events in the second age category. Similarly, select=3 for all other events in the first age category, and select=4 for all other events in the first second category. With three age categories, it works similarly with select=1,2,3 for the enter events in each age category, respectively. And select=4,5,6 for all other events.

Value

a data.frame with the columns event, time, node, dest, n, proportion, select, and shift.

See Also

individual_events.

Description

Example data to initialize a population of 1600 nodes and demonstrate various models.

Usage

data(nodes)

Format

A data.frame

Examples

## Not run: 
## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIR' model with 1600 nodes and initialize
## it to run over 4*365 days. Add one infected individual
## to the first node.
u0 <- u0_SIR()
u0$I[1] <- 1
tspan <- seq(from = 1, to = 4*365, by = 1)
model <- SIR(u0     = u0,
             tspan  = tspan,
             events = events_SIR(),
             beta   = 0.16,
             gamma  = 0.077)

## Run the model to generate a single stochastic trajectory.
result <- run(model)

## Determine nodes with one or more infected individuals in the
## trajectory. Extract the 'I' compartment and check for any
## infected individuals in each node.
infected <- colSums(trajectory(result, ~ I, format = "matrix")) > 0

## Display infected nodes in 'blue' and non-infected nodes in 'yellow'.
data("nodes", package = "SimInf")
col <- ifelse(infected, "blue", "yellow")
plot(y ~ x, nodes, col = col, pch = 20, cex = 2)

## End(Not run)

Description

The number nodes that are connected with external transfer events from each node.

Usage

outdegree(model)

Arguments

Value

vector with out-degree for each node.

Examples

## Create an 'SIR' model with 1600 nodes and initialize
## it with example data.
model <- SIR(u0 = u0_SIR(), tspan = 1:1460, events = events_SIR(),
             beta   = 0.16, gamma  = 0.077)

## Display outdegree for each node in the model.
plot(outdegree(model))

Description

Describe your model in a logical way in R, then mparse creates a SimInf_model object with your model definition that can be installed as an add-on R package.

Usage

package_skeleton(
  model,
  name = NULL,
  path = ".",
  author = NULL,
  email = NULL,
  maintainer = NULL,
  license = "GPL-3"
)

Arguments

Value

invisible NULL.

References

Read the Writing R Extensions manual for more details.

Once you have created a source package you need to install it: see the R Installation and Administration manual, INSTALL and install.packages.

Description

A matrix of scatterplots with the number of individuals in each compartment is produced. The ijth scatterplot contains x[,i] plotted against x[,j].

Usage

## S4 method for signature 'SimInf_model'
pairs(x, compartments = NULL, index = NULL, ...)

Arguments

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIR' model with 10 nodes and initialise
## it with 99 susceptible individuals and one infected
## individual. Let the model run over 100 days.
model <- SIR(u0 = data.frame(S = rep(99, 10),
                             I = rep(1, 10),
                             R = rep(0, 10)),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the model and save the result.
result <- run(model)

## Create a scatter plot that includes all compartments in all
## nodes.
pairs(result)

## Create a scatter plot that includes the S and I compartments in
## nodes 1 and 2.
pairs(result, ~S+I, 1:2)

Description

The bootstrap filtering algorithm. Systematic resampling is performed at each observation.

Usage

pfilter(model, obs_process, data, n_particles, init_model = NULL)

## S4 method for signature 'SimInf_model'
pfilter(model, obs_process, data, n_particles, init_model = NULL)

Arguments

Value

A SimInf_pfilter object.

References

N. J. Gordon, D. J. Salmond, and A. F. M. Smith. Novel Approach to Nonlinear/Non-Gaussian Bayesian State Estimation. Radar and Signal Processing, IEE Proceedings F, 140(2) 107–113, 1993. doi:10.1049/ip-f-2.1993.0015

Examples

## Not run: 
## Let us consider an SIR model in a closed population with N = 100
## individuals of whom one is initially infectious and the rest are
## susceptible. First, generate one realisation (with a specified
## seed) from the model with known parameters 'beta = 0.16' and
## 'gamma = 0.077'. Then, use 'pfilter' to apply the bootstrap
## particle algorithm on the simulated data.
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = seq(1, 100, by = 3),
             beta = 0.16,
             gamma = 0.077)

## Run the SIR model to generate simulated observed data for the
## number of infected individuals.
set.seed(22)
infected <- trajectory(run(model), "I")[, c("time", "I")]
colnames(infected) <- c("time", "Iobs")

## Use a Poison observation process for the infected individuals, such
## that 'Iobs ~ poison(I + 1e-6)'. A small constant '1e-6' is added to
## prevent numerical errors, since the simulated counts 'I' could be
## zero, which would result in the Poisson rate parameter being zero,
## which violates the conditions of the Poisson distribution. Use 1000
## particles.
pf <- pfilter(model,
              obs_process = Iobs ~ poisson(I + 1e-6),
              data = infected,
              npart = 1000)

## Print a brief summary.
pf

## Compare the number infected 'I' in the filtered trajectory with the
## infected 'Iobs' in the observed data.
plot(pf, ~I)
lines(Iobs ~ time, infected, col = "blue", lwd = 2, type = "s")

## End(Not run)

Description

Display the ABC posterior distribution

Usage

## S4 method for signature 'SimInf_abc'
plot(x, y, ...)

Arguments

Description

Display the distribution of scheduled events over time

Usage

## S4 method for signature 'SimInf_events'
plot(x, frame.plot = FALSE, ...)

Arguments

Description

Display the distribution of individual events over time

Usage

## S4 method for signature 'SimInf_indiv_events'
plot(x, frame.plot = FALSE, ...)

Arguments

Description

Plot either the median and the quantile range of the counts in all nodes, or plot the counts in specified nodes.

Usage

## S4 method for signature 'SimInf_model'
plot(
  x,
  y,
  level = 1,
  index = NULL,
  range = 0.5,
  type = "s",
  lwd = 2,
  frame.plot = FALSE,
  legend = TRUE,
  ...
)

Arguments

Examples

## Not run: 
## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIR' model with 100 nodes and initialise
## it with 990 susceptible individuals and 10 infected
## individuals in each node. Run the model over 100 days.
model <- SIR(u0 = data.frame(S = rep(990, 100),
                             I = rep(10, 100),
                             R = rep(0, 100)),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the model and save the result.
result <- run(model)

## Plot the median and interquartile range of the number
## of susceptible, infected and recovered individuals.
plot(result)

## Plot the median and the middle 95\
## number of susceptible, infected and recovered individuals.
plot(result, range = 0.95)

## Plot the median and interquartile range of the  number
## of infected individuals.
plot(result, "I")

## Use the formula notation instead to plot the median and
## interquartile range of the number of infected individuals.
plot(result, ~I)

## Plot the number of susceptible, infected
## and recovered individuals in the first
## three nodes.
plot(result, index = 1:3, range = FALSE)

## Use plot type line instead.
plot(result, index = 1:3, range = FALSE, type = "l")

## Plot the number of infected individuals in the first node.
plot(result, "I", index = 1, range = FALSE)

## Plot the proportion of infected individuals (cases)
## in the population.
plot(result, I ~ S + I + R)

## Plot the proportion of nodes with infected individuals.
plot(result, I ~ S + I + R, level = 2)

## Plot the median and interquartile range of the proportion
## of infected individuals in each node
plot(result, I ~ S + I + R, level = 3)

## Plot the proportion of infected individuals in the first
## three nodes.
plot(result, I ~ S + I + R, level = 3, index = 1:3, range = FALSE)

## End(Not run)

Description

Diagnostic plot of a particle filter object

Usage

## S4 method for signature 'SimInf_pfilter'
plot(x, y, ...)

Arguments

Description

Display the (approximate) posterior distributions obtained from fitting a particle Markov chain Monte Carlo algorithm, or the corresponding trace plots.

Usage

## S4 method for signature 'SimInf_pmcmc'
plot(x, y, start = 1, end = NULL, thin = 1, ...)

Arguments

Description

Particle Markov chain Monte Carlo (PMCMC) algorithm

Usage

pmcmc(
  model,
  obs_process,
  data,
  priors,
  n_particles,
  n_iterations,
  theta = NULL,
  covmat = NULL,
  adaptmix = 0.05,
  adaptive = 100,
  init_model = NULL,
  post_particle = NULL,
  verbose = getOption("verbose", FALSE)
)

## S4 method for signature 'SimInf_model'
pmcmc(
  model,
  obs_process,
  data,
  priors,
  n_particles,
  n_iterations,
  theta = NULL,
  covmat = NULL,
  adaptmix = 0.05,
  adaptive = 100,
  init_model = NULL,
  post_particle = NULL,
  verbose = getOption("verbose", FALSE)
)

Arguments

References

C. Andrieu, A. Doucet and R. Holenstein. Particle Markov chain Monte Carlo methods. Journal of the Royal Statistical Society, Series B 72, 269–342, 2010. doi:10.1111/j.1467-9868.2009.00736.x

G. O. Roberts and J. S. Rosenthal. Examples of adaptive MCMC. Journal of computational and graphical statistics, 18(2), 349–367, 2009. doi:10.1198/jcgs.2009.06134

See Also

continue_pmcmc.

Description

Calculate the proportion of individuals with disease in the population, or the proportion of nodes with at least one diseased individual, or the proportion of individuals with disease in each node.

Usage

prevalence(model, formula, level = 1, index = NULL, ...)

Arguments

Description

Calculate the proportion of individuals with disease in the population, or the proportion of nodes with at least one diseased individual, or the proportion of individuals with disease in each node.

Usage

## S4 method for signature 'SimInf_model'
prevalence(model, formula, level, index, format = c("data.frame", "matrix"))

Arguments

Value

A data.frame if format = "data.frame", else a matrix.

Examples

## Create an 'SIR' model with 6 nodes and initialize
## it to run over 10 days.
u0 <- data.frame(S = 100:105, I = c(0, 1, 0, 2, 0, 3), R = rep(0, 6))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)

## Run the model to generate a single stochastic trajectory.
result <- run(model)

## Determine the proportion of infected individuals (cases)
## in the population at the time-points in 'tspan'.
prevalence(result, I ~ S + I + R)

## Identical result is obtained with the shorthand 'I~.'
prevalence(result, I ~ .)

## Determine the proportion of nodes with infected individuals at
## the time-points in 'tspan'.
prevalence(result, I ~ S + I + R, level = 2)

## Determine the proportion of infected individuals in each node
## at the time-points in 'tspan'.
prevalence(result, I ~ S + I + R, level = 3)

## Determine the proportion of infected individuals in each node
## at the time-points in 'tspan' when the number of recovered is
## zero.
prevalence(result, I ~ S + I + R | R == 0, level = 3)

Description

Extract prevalence from running a particle filter

Usage

## S4 method for signature 'SimInf_pfilter'
prevalence(model, formula, level, index, format = c("data.frame", "matrix"))

Arguments

Value

A data.frame if format = "data.frame", else a matrix.

Description

Extract prevalence from the filtered trajectories from a particle Markov chain Monte Carlo algorithm.

Usage

## S4 method for signature 'SimInf_pmcmc'
prevalence(model, formula, level, index, start = 1, end = NULL, thin = 1)

Arguments

Value

A data.frame where the first column is the iteration and the remaining columns are the result from calling prevalence,SimInf_model-method with the arguments formula, level and index for each iteration.

Description

Using a sparse result matrix can save a lot of memory if the model contains many nodes and time-points, but where only a few of the data points are of interest for post-processing.

Usage

punchcard(model) <- value

## S4 replacement method for signature 'SimInf_model'
punchcard(model) <- value

Arguments

Details

Using a sparse result matrix can save a lot of memory if the model contains many nodes and time-points, but where only a few of the data points are of interest for post-processing. To use this feature, a template has to be defined for which data points to record. This is done using a data.frame that specifies the time-points (column ‘time’) and nodes (column ‘node’) to record the state of the compartments, see ‘Examples’. The specified time-points, nodes and compartments must exist in the model, or an error is raised. Note that specifying a template only affects which data-points are recorded for post-processing, it does not affect how the solver simulates the trajectory.

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIR' model with 6 nodes and initialize it to run over 10 days.
u0 <- data.frame(S = 100:105, I = 1:6, R = rep(0, 6))
model <- SIR(u0 = u0, tspan = 1:10, beta = 0.16, gamma = 0.077)

## Run the model.
result <- run(model)

## Display the trajectory with data for every node at each
## time-point in tspan.
trajectory(result)

## Assume we are only interested in nodes '2' and '4' at the
## time-points '3' and '5'
df <- data.frame(time = c(3, 5, 3, 5),
                 node = c(2, 2, 4, 4),
                 S = c(TRUE, TRUE, TRUE, TRUE),
                 I = c(TRUE, TRUE, TRUE, TRUE),
                 R = c(TRUE, TRUE, TRUE, TRUE))
punchcard(model) <- df
result <- run(model)
trajectory(result)

## We can also specify to record only some of the compartments in
## each time-step.
df <- data.frame(time = c(3, 5, 3, 5),
                 node = c(2, 2, 4, 4),
                 S = c(FALSE, TRUE, TRUE, TRUE),
                 I = c(TRUE, FALSE, TRUE, FALSE),
                 R = c(TRUE, FALSE, TRUE, TRUE))
punchcard(model) <- df
result <- run(model)
trajectory(result)

## A shortcut to specify to record all of the compartments in
## each time-step is to only inlude node and time.
df <- data.frame(time = c(3, 5, 3, 5),
                 node = c(2, 2, 4, 4))
punchcard(model) <- df
result <- run(model)
trajectory(result)

## It is possible to use an empty 'data.frame' to specify
## that no data-points should be recorded for the trajectory.
punchcard(model) <- data.frame()
result <- run(model)
trajectory(result)

## Use 'NULL' to reset the model to record data for every node at
## each time-point in tspan.
punchcard(model) <- NULL
result <- run(model)
trajectory(result)

Description

Run the SimInf stochastic simulation algorithm

Usage

run(model, ...)

## S4 method for signature 'SimInf_model'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SEIR'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SIR'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SIS'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SISe'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SISe3'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SISe3_sp'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SISe_sp'
run(model, solver = c("ssm", "aem"), ...)

## S4 method for signature 'SimInf_abc'
run(model, ...)

Arguments

Value

SimInf_model object with result from simulation.

References

S. Widgren, P. Bauer, R. Eriksson and S. Engblom. SimInf: An R Package for Data-Driven Stochastic Disease Spread Simulations. Journal of Statistical Software, 91(12), 1–42. doi:10.18637/jss.v091.i12. An updated version of this paper is available as a vignette in the package.

P. Bauer, S. Engblom and S. Widgren. Fast Event-Based Epidemiological Simulations on National Scales. International Journal of High Performance Computing Applications, 30(4), 438–453, 2016. doi: 10.1177/1094342016635723

P. Bauer and S. Engblom. Sensitivity Estimation and Inverse Problems in Spatial Stochastic Models of Chemical Kinetics. In: A. Abdulle, S. Deparis, D. Kressner, F. Nobile and M. Picasso (eds.), Numerical Mathematics and Advanced Applications - ENUMATH 2013, pp. 519–527, Lecture Notes in Computational Science and Engineering, vol 103. Springer, Cham, 2015. doi:10.1007/978-3-319-10705-9_51

Examples

## For reproducibility, call the set.seed() function and specify
## the number of threads to use. To use all available threads,
## remove the set_num_threads() call.
set.seed(123)
set_num_threads(1)

## Create an 'SIR' model with 10 nodes and initialise
## it to run over 100 days.
model <- SIR(u0 = data.frame(S = rep(99, 10),
                             I = rep(1, 10),
                             R = rep(0, 10)),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the model and save the result.
result <- run(model)

## Plot the proportion of susceptible, infected and recovered
## individuals.
plot(result)

Description

Create an SEIR model to be used by the simulation framework.

Usage

SEIR(u0, tspan, events = NULL, beta = NULL, epsilon = NULL, gamma = NULL)

Arguments

Details

The SEIR model contains four compartments; number of susceptible (S), number of exposed (E) (those who have been infected but are not yet infectious), number of infectious (I), and number of recovered (R). Moreover, it has three state transitions,

$S \stackrel{\beta S I / N}{\longrightarrow} E$

$E \stackrel{\epsilon E}{\longrightarrow} I$

$I \stackrel{\gamma I}{\longrightarrow} R$

where $\beta$ is the transmission rate, $\epsilon$ is the incubation rate, $\gamma$ is the recovery rate, and $N=S+E+I+R$.

The argument u0 must be a data.frame with one row for each node with the following columns:

S

The number of sucsceptible in each node

E

The number of exposed in each node

I

The number of infected in each node

R

The number of recovered in each node

Value

A SimInf_model of class SEIR

Examples

## Create a SEIR model object.
model <- SEIR(u0 = data.frame(S = 99, E = 0, I = 1, R = 0),
              tspan = 1:100,
              beta = 0.16,
              epsilon = 0.25,
              gamma = 0.077)

## Run the SEIR model and plot the result.
set.seed(3)
result <- run(model)
plot(result)

Description

Class to handle the SEIR SimInf_model.

Description

Utility function to extract events@E from a SimInf_model object, see SimInf_events

Usage

select_matrix(model)

## S4 method for signature 'SimInf_model'
select_matrix(model)

Arguments

Value

dgCMatrix object.

Examples

## Create an SIR model
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:5, beta = 0.16, gamma = 0.077)

## Extract the select matrix from the model
select_matrix(model)

Description

Utility function to set events@E in a SimInf_model object, see SimInf_events

Usage

select_matrix(model) <- value

## S4 replacement method for signature 'SimInf_model'
select_matrix(model) <- value

Arguments

Examples

## Create an SIR model
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:5, beta = 0.16, gamma = 0.077)

## Set the select matrix
select_matrix(model) <- matrix(c(1, 0, 0, 1, 1, 1, 0, 0, 1), nrow = 3)

## Extract the select matrix from the model
select_matrix(model)

Description

Set the number of threads to be used in SimInf code that is parallelized with OpenMP (if available). The number of threads is initialized when SimInf is first loaded in the R session using optional envioronment variables (see ‘Details’). It is also possible to specify the number of threads by calling set_num_threads. If the environment variables that affect the number of threads change, then set_num_threads must be called again for it to take effect.

Usage

set_num_threads(threads = NULL)

Arguments

Details

The omp_get_num_procs() function is used to determine the number of processors that are available to the device at the time the routine is called. The number of threads is then limited by omp_get_thread_limit() and the current values of the environmental variables (if set)

• Sys.getenv("OMP_THREAD_LIMIT")

• Sys.getenv("OMP_NUM_THREADS")

• Sys.getenv("SIMINF_NUM_THREADS")

Additionally, the maximum number of threads can be controlled by the threads argument, given that its value is not above any of the limits described above.

Value

The previous value is returned (invisible).

Description

Utility function to extract the shift matrix events@N from a SimInf_model object, see SimInf_events

Usage

shift_matrix(model)

## S4 method for signature 'SimInf_model'
shift_matrix(model)

Arguments

Value

A mtrix.

Examples

## Create an SIR model
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:5, beta = 0.16, gamma = 0.077)

## Extract the shift matrix from the model
shift_matrix(model)

Description

Utility function to set events@N in a SimInf_model object, see SimInf_events

Usage

shift_matrix(model) <- value

## S4 replacement method for signature 'SimInf_model'
shift_matrix(model) <- value

Arguments

Value

SimInf_model object

Examples

## Create an SIR model
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:5, beta = 0.16, gamma = 0.077)

## Set the shift matrix
shift_matrix(model) <- matrix(c(2, 1, 0), nrow = 3)

## Extract the shift matrix from the model
shift_matrix(model)

Description

Print summary of a SimInf_abc object

Usage

## S4 method for signature 'SimInf_abc'
show(object)

Arguments

Value

invisible(object).

Description

Shows the number of scheduled events.

Usage

## S4 method for signature 'SimInf_events'
show(object)

Arguments

Value

None (invisible 'NULL').

Description

Print summary of a SimInf_indiv_events object

Usage

## S4 method for signature 'SimInf_indiv_events'
show(object)

Arguments

Value

invisible(object).

Description

Brief summary of SimInf_model

Usage

## S4 method for signature 'SimInf_model'
show(object)

Arguments

Value

None (invisible 'NULL').

Examples

## Create an 'SIR' model with 10 nodes and initialise
## it to run over 100 days.
model <- SIR(u0 = data.frame(S = rep(99, 10),
                             I = rep(1, 10),
                             R = rep(0, 10)),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Brief summary of the model
model

## Run the model and save the result
result <- run(model)

## Brief summary of the result. Note that 'U' and 'V' are
## non-empty after running the model.
result

Description

Brief summary of a SimInf_pfilter object

Usage

## S4 method for signature 'SimInf_pfilter'
show(object)

Arguments

Value

invisible(object).

Description

Brief summary of a SimInf_pmcmc object

Usage

## S4 method for signature 'SimInf_pmcmc'
show(object)

Arguments

Value

invisible(object).

Description

The SimInf package provides a flexible framework for data-driven spatio-temporal disease spread modeling, designed to efficiently handle population demographics and network data. The framework integrates infection dynamics in each subpopulation as continuous-time Markov chains (CTMC) using the Gillespie stochastic simulation algorithm (SSA) and incorporates available data such as births, deaths or movements as scheduled events. A scheduled event is used to modify the state of a subpopulation at a predefined time-point.

Details

The SimInf_model is central and provides the basis for the framework. A SimInf_model object supplies the state-change matrix, the dependency graph, the scheduled events, and the initial state of the system.

All predefined models in SimInf have a generating function, with the same name as the model, for example SIR.

A model can also be created from a model specification using the mparse method.

After a model is created, a simulation is started with a call to the run method and if execution is successful, it returns a modified SimInf_model object with a single stochastic solution trajectory attached to it.

SimInf provides several utility functions to inspect simulated data, for example, show, summary and plot. To facilitate custom analysis, it provides the trajectory,SimInf_model-method and prevalence methods.

One of our design goal was to make SimInf extendable and enable usage of the numerical solvers from other R extension packages in order to facilitate complex epidemiological research. To support this, SimInf has functionality to generate the required C and R code from a model specification, see package_skeleton

Author(s)

Maintainer: Stefan Widgren [email protected] (ORCID)

Authors:

• Robin Eriksson (ORCID)

• Stefan Engblom (ORCID)

• Pavol Bauer (ORCID)

Other contributors:

• Thomas Rosendal (ORCID) [contributor]

• Ivana Rodriguez Ewerlöf (ORCID) [contributor]

• Attractive Chaos (Author of 'kvec.h'.) [copyright holder]

References

S. Widgren, P. Bauer, R. Eriksson and S. Engblom. SimInf: An R Package for Data-Driven Stochastic Disease Spread Simulations. Journal of Statistical Software, 91(12), 1–42. doi:10.18637/jss.v091.i12. An updated version of this paper is available as a vignette in the package.

See Also

Useful links:

• https://github.com/stewid/SimInf

• http://stewid.github.io/SimInf/

• Report bugs at https://github.com/stewid/SimInf/issues

Description

Class "SimInf_abc"

Slots

model

The SimInf_model object to estimate parameters in.

priors

A data.frame containing the four columns parameter, distribution, p1 and p2. The column parameter gives the name of the parameter referred to in the model. The column distribution contains the name of the prior distribution. Valid distributions are 'gamma', 'normal' or 'uniform'. The column p1 is a numeric vector with the first hyperparameter for each prior: 'gamma') shape, 'lognormal') logmean, 'normal') mean, and 'uniform') lower bound. The column p2 is a numeric vector with the second hyperparameter for each prior: 'gamma') rate, 'lognormal') standard deviation on the log scale, 'normal') standard deviation, and 'uniform') upper bound.

target

Character vector (gdata or ldata) that determines if the ABC-SMC method estimates parameters in model@gdata or in model@ldata.

pars

Index to the parameters in target.

nprop

An integer vector with the number of simulated proposals in each generation.

fn

A function for calculating the summary statistics for the simulated trajectory and determine the distance for each particle, see abc for more details.

tolerance

A numeric matrix (number of summary statistics $\times$ number of generations) where each column contains the tolerances for a generation and each row contains a sequence of gradually decreasing tolerances.

x

A numeric array (number of particles $\times$ number of parameters $\times$ number of generations) with the parameter values for the accepted particles in each generation. Each row is one particle.

weight

A numeric matrix (number of particles $\times$ number of generations) with the weights for the particles x in the corresponding generation.

distance

A numeric array (number of particles $\times$ number of summary statistics $\times$ number of generations) with the distance for the particles x in each generation. Each row contains the distance for a particle and each column contains the distance for a summary statistic.

ess

A numeric vector with the effective sample size (ESS) in each generation. The effective sample size is computed as

$\left(\sum_{i=1}^N\!(w_{g}^{(i)})^2\right)^{-1},$

where $w_{g}^{(i)}$ is the normalized weight of particle $i$ in generation $g$.

init_model

An optional function that, if non-NULL, is applied before running each proposal. The function must accept one argument of type SimInf_model with the current model of the fitting process. This function can be useful to specify the initial state of u0 or v0 of the model before running a trajectory with proposed parameters.

See Also

abc and continue_abc.

Create a SimInf_events object

Description

The argument events must be a data.frame with the following columns:

event

Four event types are supported by the current solvers: exit, enter, internal transfer, and external transfer. When assigning the events, they can either be coded as a numerical value or a character string: exit; 0 or 'exit', enter; 1 or 'enter', internal transfer; 2 or 'intTrans', and external transfer; 3 or 'extTrans'. Internally in SimInf, the event type is coded as a numerical value.

time

When the event occurs i.e., the event is processed when time is reached in the simulation. Can be either an integer or a Date vector. A Date vector is coerced to a numeric vector as days, where t0 determines the offset to match the time of the events to the model tspan vector.

node

The node that the event operates on. Also the source node for an external transfer event. 1 <= node[i] <= Number of nodes.

dest

The destination node for an external transfer event i.e., individuals are moved from node to dest, where 1 <= dest[i] <= Number of nodes. Set event = 0 for the other event types. dest is an integer vector.

n

The number of individuals affected by the event. n[i] >= 0.

proportion

If n[i] equals zero, the number of individuals affected by event[i] is calculated by sampling the number of individuals from a binomial distribution using the proportion[i] and the number of individuals in the compartments. Numeric vector. 0 <= proportion[i] <= 1.

select

To process an event[i], the compartments affected by the event are specified with select[i] together with the matrix E, where select[i] determines which column in E to use. The specific individuals affected by the event are sampled from the compartments corresponding to the non-zero entries in the specified column in E[, select[i]], where select is an integer vector.

shift

Determines how individuals in internal transfer and external transfer events are shifted to enter another compartment. The sampled individuals are shifted according to column shift[i] in matrix N i.e., N[, shift[i]], where shift is an integer vector. See above for a description of N. Unsued for the other event types.

Usage

SimInf_events(E = NULL, N = NULL, events = NULL, t0 = NULL)

Arguments

Value

S4 class SimInf_events

Examples

## Let us illustrate how movement events can be used to transfer
## individuals from one node to another.  Use the built-in SIR
## model and start with 2 nodes where all individuals are in the
## first node (100 per compartment).
u0 <- data.frame(S = c(100, 0), I = c(100, 0), R = c(100, 0))

## Then create 300 movement events to transfer all individuals,
## one per day, from the first node to the second node. Use the
## fourth column in the select matrix where all compartments
## can be sampled with equal weight.
events <- data.frame(event      = rep("extTrans", 300),
                     time       = 1:300,
                     node       = 1,
                     dest       = 2,
                     n          = 1,
                     proportion = 0,
                     select     = 4,
                     shift      = 0)

## Create an SIR model without disease transmission to
## demonstrate the events.
model <- SIR(u0      = u0,
             tspan  = 1:300,
             events = events,
             beta   = 0,
             gamma  = 0)

## Run the model and plot the number of individuals in
## the second node.  As can be seen in the figure, all
## indivuduals have been moved to the second node when
## t = 300.
plot(run(model), index = 1:2, range = FALSE)

## Let us now double the weight to sample from the 'I'
## compartment and rerun the model.
model@events@E[2, 4] <- 2
plot(run(model), index = 1:2, range = FALSE)

## And much larger weight to sample from the I compartment.
model@events@E[2, 4] <- 10
plot(run(model), index = 1:2, range = FALSE)

## Increase the weight for the R compartment.
model@events@E[3, 4] <- 4
plot(run(model), index = 1:2, range = FALSE)

Description

Class to hold data for scheduled events to modify the discrete state of individuals in a node at a pre-defined time t.

Slots

E

Each row corresponds to one compartment in the model. The non-zero entries in a column indicates the compartments to include in an event. For the exit, internal transfer and external transfer events, a non-zero entry indicate the compartments to sample individuals from. For the enter event, all individuals enter first non-zero compartment. E is sparse matrix of class dgCMatrix.

N

Determines how individuals in internal transfer and external transfer events are shifted to enter another compartment. Each row corresponds to one compartment in the model. The values in a column are added to the current compartment of sampled individuals to specify the destination compartment, for example, a value of 1 in an entry means that sampled individuals in this compartment are moved to the next compartment. Which column to use for each event is specified by the shift vector (see below). N is an integer matrix.

event

Type of event: 0) exit, 1) enter, 2) internal transfer, and 3) external transfer. Other values are reserved for future event types and not supported by the current solvers. Integer vector.

time

Time of when the event occurs i.e., the event is processed when time is reached in the simulation. time is an integer vector.

node

The node that the event operates on. Also the source node for an external transfer event. Integer vector. 1 <= node[i] <= Number of nodes.

dest

The destination node for an external transfer event i.e., individuals are moved from node to dest, where 1 <= dest[i] <= Number of nodes. Set event = 0 for the other event types. dest is an integer vector.

n

The number of individuals affected by the event. Integer vector. n[i] >= 0.

proportion

If n[i] equals zero, the number of individuals affected by event[i] is calculated by sampling the number of individuals from a binomial distribution using the proportion[i] and the number of individuals in the compartments. Numeric vector. 0 <= proportion[i] <= 1.

select

To process event[i], the compartments affected by the event are specified with select[i] together with the matrix E, where select[i] determines which column in E to use. The specific individuals affected by the event are proportionally sampled from the compartments corresponding to the non-zero entries in the specified column in E[, select[i]], where select is an integer vector.

shift

Determines how individuals in internal transfer and external transfer events are shifted to enter another compartment. The sampled individuals are shifted according to column shift[i] in matrix N i.e., N[, shift[i]], where shift is an integer vector. See above for a description of N. Unsued for the other event types.

Description

Class "SimInf_indiv_events"

Slots

id

an integer or character identifier of the individual.

event

four event types are supported: exit, enter, internal transfer, and external transfer. When assigning the events, they can either be coded as a numerical value or a character string: exit; 0 or 'exit', enter; 1 or 'enter', internal transfer; 2 or 'intTrans', and external transfer; 3 or 'extTrans'.

time

an integer, character, or date (of class Date) for when the event occured. If it's a character it must be able to coerce to Date.

node

an integer or character identifier of the source node.

dest

an integer or character identifier of the destination node.

Description

Create a SimInf_model

Usage

SimInf_model(
  G,
  S,
  tspan,
  events = NULL,
  ldata = NULL,
  gdata = NULL,
  U = NULL,
  u0 = NULL,
  v0 = NULL,
  V = NULL,
  E = NULL,
  N = NULL,
  C_code = NULL
)

Arguments

Value

SimInf_model

Description

Class to handle data for the SimInf_model.

Slots

G

Dependency graph that indicates the transition rates that need to be updated after a given state transition has occured. A non-zero entry in element G[i, i] indicates that transition rate i needs to be recalculated if the state transition j occurs. Sparse matrix ($Nt \times Nt$) of object class dgCMatrix.

S

Each column corresponds to a state transition, and execution of state transition j amounts to adding the S[, j] column to the state vector u[, i] of node i where the transition occurred. Sparse matrix ($Nc \times Nt$) of object class dgCMatrix.

U

The result matrix with the number of individuals in each compartment in every node. U[, j] contains the number of individuals in each compartment at tspan[j]. U[1:Nc, j] contains the number of individuals in node 1 at tspan[j]. U[(Nc + 1):(2 * Nc), j] contains the number of individuals in node 2 at tspan[j] etc. Integer matrix ($N_n N_c \times$ length(tspan)).

U_sparse

If the model was configured to write the solution to a sparse matrix (dgCMatrix) the U_sparse contains the data and U is empty. The layout of the data in U_sparse is identical to U. Please note that U_sparse is numeric and U is integer.

V

The result matrix for the real-valued continuous state. V[, j] contains the real-valued state of the system at tspan[j]. Numeric matrix ($N_n$dim(ldata)[1] $\times$ length(tspan)).

V_sparse

If the model was configured to write the solution to a sparse matrix (dgCMatrix) the V_sparse contains the data and V is empty. The layout of the data in V_sparse is identical to V.

ldata

A matrix with local data for the nodes. The column ldata[, j] contains the local data vector for the node j. The local data vector is passed as an argument to the transition rate functions and the post time step function.

gdata

A numeric vector with global data that is common to all nodes. The global data vector is passed as an argument to the transition rate functions and the post time step function.

tspan

A vector of increasing time points where the state of each node is to be returned.

u0

The initial state vector ($N_c \times N_n$) with the number of individuals in each compartment in every node.

v0

The initial value for the real-valued continuous state. Numeric matrix (dim(ldata)[1] $\times N_n$).

events

Scheduled events SimInf_events

replicates

Number of replicates of the model.

C_code

Character vector with optional model C code. If non-empty, the C code is written to a temporary C-file when the run method is called. The temporary C-file is compiled and the resulting DLL is dynamically loaded. The DLL is unloaded and the temporary files are removed after running the model.

Description

Class "SimInf_pfilter"

Slots

model

A SimInf_model object with one filtered trajectory attached.

n_particles

An integer with the number of particles that was used at each timestep.

loglik

The estimated log likelihood.

ess

A numeric vector with the effective sample size (ESS). The effective sample size is computed as

$\left(\sum_{i=1}^N\!(w_{t}^{i})^2\right)^{-1},$

where $w_{t}^{i}$ is the normalized weight of particle $i$ at time $t$.

Description

Class "SimInf_pmcmc"

Slots

model

The SimInf_model object to estimate parameters in.

priors

A data.frame containing the four columns parameter, distribution, p1 and p2. The column parameter gives the name of the parameter referred to in the model. The column distribution contains the name of the prior distribution. Valid distributions are 'gamma', 'normal' or 'uniform'. The column p1 is a numeric vector with the first hyperparameter for each prior: 'gamma') shape, 'lognormal') logmean, 'normal') mean, and 'uniform') lower bound. The column p2 is a numeric vector with the second hyperparameter for each prior: 'gamma') rate, 'lognormal') standard deviation on the log scale, 'normal') standard deviation, and 'uniform') upper bound.

target

Character vector (gdata or ldata) that determines if the pmcmc method estimates parameters in model@gdata or in model@ldata.

pars

Index to the parameters in target.

n_particles

An integer with the number of particles (> 1) to use in the bootstrap particle filter.

data

A data.frame holding the time series data for the observation process.

chain

A matrix where each row contains logPost, logLik, logPrior, accept, and the parameters for each iteration.

covmat

A named numeric (npars x npars) matrix with covariances to use as initial proposal matrix.

adaptmix

Mixing proportion for adaptive proposal.

adaptive

Controls when to start adaptive update.

See Also

pmcmc and continue_pmcmc.

Description

Create an SIR model to be used by the simulation framework.

Usage

SIR(u0, tspan, events = NULL, beta = NULL, gamma = NULL)

Arguments

Details

The SIR model contains three compartments; number of susceptible (S), number of infectious (I), and number of recovered (R). Moreover, it has two state transitions,

$S \stackrel{\beta S I / N}{\longrightarrow} I$

$I \stackrel{\gamma I}{\longrightarrow} R$

where $\beta$ is the transmission rate, $\gamma$ is the recovery rate, and $N=S+I+R$.

The argument u0 must be a data.frame with one row for each node with the following columns:

S

The number of sucsceptible in each node

I

The number of infected in each node

R

The number of recovered in each node

Value

A SimInf_model of class SIR

Examples

## Create an SIR model object.
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the SIR model and plot the result.
set.seed(22)
result <- run(model)
plot(result)

Description

Class to handle the SIR SimInf_model.

Details

The SIR model contains three compartments; number of susceptible (S), number of infectious (I), and number of recovered (R). Moreover, it has two state transitions,

$S \stackrel{\beta S I / N}{\longrightarrow} I$

$I \stackrel{\gamma I}{\longrightarrow} R$

where $\beta$ is the transmission rate, $\gamma$ is the recovery rate, and $N=S+I+R$.

Examples

## Create an SIR model object.
model <- SIR(u0 = data.frame(S = 99, I = 1, R = 0),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the SIR model and plot the result.
set.seed(22)
result <- run(model)
plot(result)

Description

Create an SIS model to be used by the simulation framework.

Usage

SIS(u0, tspan, events = NULL, beta = NULL, gamma = NULL)

Arguments

Details

The SIS model contains two compartments; number of susceptible (S), and number of infectious (I). Moreover, it has two state transitions,

$S \stackrel{\beta S I / N}{\longrightarrow} I$

$I \stackrel{\gamma I}{\longrightarrow} S$

where $\beta$ is the transmission rate, $\gamma$ is the recovery rate, and $N=S+I$.

The argument u0 must be a data.frame with one row for each node with the following columns:

S

The number of sucsceptible in each node

I

The number of infected in each node

Value

A SimInf_model of class SIS

Examples

## Create an SIS model object.
model <- SIS(u0 = data.frame(S = 99, I = 1),
             tspan = 1:100,
             beta = 0.16,
             gamma = 0.077)

## Run the SIS model and plot the result.
set.seed(22)
result <- run(model)
plot(result)

Description

Class to handle the SIS SimInf_model.

Details

The SIS model contains two compartments; number of susceptible (S), and number of infectious (I). Moreover, it has two state transitions,

$S \stackrel{\beta S I / N}{\longrightarrow} I$

$I \stackrel{\gamma I}{\longrightarrow} S$