## Introduction

In a model, parameters are the input variables used to define aspects of the system behavior. In the basic built-in SIS (Susceptible-Infected-Susceptible) model, these parameters could be the act rate, the infection probability and the recovery rate. In simple models, each of these parameters are single fixed values that do not change over the course of a simulation. In more complex models, we may want more flexibility in model parameterization.

Therefore, in this vignette, we demonstrate how to implement:

• Scenarios: sets of parameters to be changed for a simulation, either at the start or at a specific timestep.
• Random parameters: distributions of possible values rather than a single fixed value.
• Time-varying parameters and control settings: The scenarios functionality provides the most straightforward way to implement time-varying parameters, but more direct functionality is demonstrated here for advanced users seeking to implement time-varying parameters or control settings.

## Scenarios

Scenarios can be defined as a single set of parameters to be used for a particular model run. For this example we will use a simple SIS model to demonstrate how scenarios work. First, we set up the model as we would normally.

set.seed(10)

nw <- network_initialize(n = 200)
est <- netest(nw,
formation = ~edges, target.stats = 60,
coef.diss = dissolution_coefs(~offset(edges), 10, 0),
verbose = FALSE
)
#> Starting maximum pseudolikelihood estimation (MPLE):
#> Evaluating the predictor and response matrix.
#> Maximizing the pseudolikelihood.
#> Finished MPLE.
#> Stopping at the initial estimate.

param <- param.net(inf.prob = 0.9, rec.rate = 0.01, act.rate = 2)
init <- init.net(i.num = 10)
control <- control.net(type = "SIS", nsims = 1, nsteps = 250, verbose = FALSE)

### Scenario Definitions

To define the scenarios, we will make use of the EpiModel::create_scenario_list function. It takes a specially formatted data.frame as input and outputs a list of scenarios usable by EpiModel.

We will use the tibble::tribble function to create the data.frame here. But any data frame construction methods will work as long as the resulting data.frame is properly formatted.

suppressMessages(library(dplyr))

scenarios.df <- tribble(
~.scenario.id, ~.at, ~inf.prob, ~rec.rate,
"base", 0, 0.9, 0.01,
"initial_change", 0, 0.2, 0.01,
"multiple_changes", 0, 0.1, 0.04,
"multiple_changes", 100, 0.9, 0.01,
"multiple_changes", 200, 0.1, 0.1
)

knitr::kable(scenarios.df)
.scenario.id .at inf.prob rec.rate
base 0 0.9 0.01
initial_change 0 0.2 0.01
multiple_changes 0 0.1 0.04
multiple_changes 100 0.9 0.01
multiple_changes 200 0.1 0.10

A data.frame of scenarios must be formatted as follows. First, there are two required columns: - .scenario.id: an identifier for the scenario - .at: when should the changes apply during the simulation

Any number of parameters columns. They must start with a letter and contain only letters, numbers or .. Underscores _ are not accepted. If a cell is left empty (NA), EpiModel will assume that you want to set the parameter to NA. NA does not mean that the value should remain unmodified.

In the data.frame above, three rows share the same .scenario.id. This means that EpiModel will consider it a single scenario were multiple events occur during the simulation.

To go from the scenarios.df to a list of usable scenarios we run:

scenarios.list <- create_scenario_list(scenarios.df)
str(scenarios.list, max.level = 2)
#> List of 3
#>  $base :List of 2 #> ..$ id                 : chr "base"
#>   ..$.param.updater.list:List of 1 #>$ initial_change  :List of 2
#>   ..$id : chr "initial_change" #> ..$ .param.updater.list:List of 1
#>  $multiple_changes:List of 2 #> ..$ id                 : chr "multiple_changes"
#>   ..$.param.updater.list:List of 3 We now have a named list of scenarios. We will run all the scenarios, combine there results and plot the number of susceptible and infected over time to show what happened. To do so, we loop over all 3 scenarios and use EpiModel::use_scenario function to create a new sc.param object to be used for the simulation parameters. # creation of a list that will hold the result of the simulations d_list <- vector(mode = "list", length = length(scenarios.list)) names(d_list) <- names(scenarios.list) for (scenario in scenarios.list) { # the id of the scenario is stored into scenario$id
print(scenario$id) sc.param <- use_scenario(param, scenario) sim <- netsim(est, sc.param, init, control) # conversion to data.frame and storage of the scenario name in it d_sim <- as.data.frame(sim) d_sim[["scenario"]] <- scenario$id
#> [1] 10
#>
#> $param #>$param$inf.prob #> [1] 0.3 #> #>$param$act.rate #> [1] 0.5 This defines an updater that will change the value of the inf.prob parameter to 0.3 and the value to the act.rate parameter to 0.5. This change will happen at time step 10. If we want, multiple parameter changes, then we can write a list of updaters: # Create a list.of.updaters list.of.updaters <- list( # this is one updater list( at = 100, param = list( inf.prob = 0.3, act.rate = 0.3 ) ), # this is another updater list( at = 125, param = list( inf.prob = 0.01 ) ) ) In this example, we define two updaters, one that occurs at time step 100 and the other one at time step 125. As demonstrated below, these values for inf.prob and act.rate will change from the values of 0.1 established by param.netat the beginning of the time series. ### Incorporating Updaters Below we set up a complete example with a closed population SI model using the parameters and updaters defined above. Note that the parameter updaters get passed into param.net with the special .param.updater.list argument. param <- param.net( inf.prob = 0.1, act.rate = 0.1, .param.updater.list = list.of.updaters ) init <- init.net(i.num = 10) control <- control.net( type = "SI", nsims = 1, nsteps = 200, verbose = FALSE ) Next, we run the model with a simple network structure: nw <- network_initialize(n = 100) est <- netest( nw, formation = ~edges, target.stats = 50, coef.diss = dissolution_coefs(~offset(edges), 10, 0), verbose = FALSE ) #> Starting maximum pseudolikelihood estimation (MPLE): #> Evaluating the predictor and response matrix. #> Maximizing the pseudolikelihood. #> Finished MPLE. #> Stopping at the initial estimate. mod <- netsim(est, param, init, control) #> #> A MESSAGE occured in module 'epimodel.internal' at step 100 #> #> #> At timestep = 100 the following parameters were modified: #> 'inf.prob', 'act.rate' #> #> A MESSAGE occured in module 'epimodel.internal' at step 125 #> #> #> At timestep = 125 the following parameters were modified: #> 'inf.prob' The model plot here demonstrates the inflection points in the epidemic trajectory at the two parameter change points defined above. plot(mod, mean.smooth = FALSE) abline(v = c(100, 125), lty = 2, col = "grey50") #### Verbosity When creating a parameter updater, one can add an optional verbose element to the list. If TRUE, EpiModel will output a message describing what changes were performed and when they occurred. list( at = 10, param = list( inf.prob = 0.3, act.rate = 0.5 ), verbose = TRUE ) #>$at
#> [1] 10
#>
#> $param #>$param$inf.prob #> [1] 0.3 #> #>$param$act.rate #> [1] 0.5 #> #> #>$verbose
#> [1] TRUE

#### Relative Parameter Changes

It may be useful to configure the changes with respect to the current value instead of a fixed new value. This is possible as demonstrated below for inf.prob.

list(
at = 10,
param = list(
inf.prob = function(x) plogis(qlogis(x) + log(2)),
act.rate = 0.5
)
)
#> $at #> [1] 10 #> #>$param
#> $param$inf.prob
#> function(x) plogis(qlogis(x) + log(2))
#>
#> $param$act.rate
#> [1] 0.5

This updater will set the value of act.rate to 0.5 as before. But, for inf.prob we define a function (not a function call). In this case, the will apply the function to the current value of act.rate. If we consider as in the previous example that act.rate is set to 0.1 by param.net, then its new value will be obtained by adding the log odds of 2 to the original value of inf.prob: plogis(qlogis(0.1) + log(2)): 0.1818182. Whereas the previous value of inf.prob was 0.10, the updated value will be 0.18.

### Time-Varying Control Settings

Similar to time-varying parameters, we can use time-varying controls. These work in the same way. Each control updater is defined as a list of lists.

# Create a list.of.updaters
list.of.updaters <- list(
# this is one updater
list(
at = 100,
control = list(
resimulate.network = FALSE
)
),
# this is another updater
list(
at = 125,
control = list(
verbose = FALSE
)
)
)

This example sets two control updaters, one that turns off network resimulation at timestep 100 and another toggling of the model verbosity at step 125.

The list.of.updaters then gets passed into control.net through the special .control.updater.list argument.

control <- control.net(
type = "SI",
nsims = 1,
nsteps = 200,
verbose = TRUE,
.control.updater.list = list.of.updaters
)`