Working With Model Parameters

In EpiModel: Mathematical Modeling of Infectious Disease Dynamics

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(EpiModel)

Introduction

This vignette discusses mechanisms usable inside EpiModel network models with custom extension modules. More information about these may be found in the "New Network Models with EpiModel" section of the EpiModel Tutorials.

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:

• Random parameters, which have a distribution of possible values rather than a single fixed value.
• Time-varying parameters, which may change at specific time steps during the simulation.

Random Parameters

The Model

First, we define a simple SI model.

nw <- network_initialize(n = 50)

est <- netest(
  nw, formation = ~edges,
  target.stats = 25,
  coef.diss = dissolution_coefs(~offset(edges), 10, 0),
  verbose = FALSE
)

param <- param.net(
  inf.prob = 0.3,
  act.rate = 0.5,
  dummy.param = 4,
  dummy.strat.param = c(0, 1)
)

init <- init.net(i.num = 10)
control <- control.net(type = "SI", nsims = 1, nsteps = 5, verbose = FALSE)
mod <- netsim(est, param, init, control)
mod

In the parameters, we set the value for inf.prob and act.rate as fixed, but we also define dummy.param and dummy.strat.param. These latter two parameters will serve to illustrate random parameters. dummy.strat.param has two elements; this may be used, for example, to stratify a parameter by subpopulation.

The last line prints a summary of the model. In it we can see the value of the parameters under the "Fixed Parameters" section. Note the additional groups parameter defined automatically by EpiModel as part of the "SI" model definition.

Adding Random Parameters

To allow our parameters to be drawn from a distribution of random values, we use the random.params argument to param.net. There are two ways of defining which parameters are random, and the distribution of values for those random parameters to draw randomly and how to draw them.

Generator Functions

The first option is to define a generator function for each parameter we want to treat as random.

my.randoms <- list(
  act.rate = param_random(c(0.25, 0.5, 0.75)),
  dummy.param = function() rbeta(1, 1, 2),
  dummy.strat.param = function() c(
    rnorm(1, 0.05, 0.01),
    rnorm(1, 0.15, 0.03)
  )
)

param <- param.net(
  inf.prob = 0.3,
  random.params = my.randoms
)

param

Here we kept the inf.prob parameter fixed at 0.3 and defined a list object called my.randoms containing 3 elements:

• act.rate uses the param_random function factory defined by EpiModel (see ?EpiModel::param_random).
• dummy.param is a function with no argument that returns a random value from a beta distribution.
• dummy.strat.param is a function with no argument that returns 2 values sampled from normal distributions, each with different means and standard deviations.

Each element is named after the parameter it will fill and MUST BE a function taking no argument and outputting a vector of the right size for the parameter: size 1 for act.rate and dummy.param; size 2 for dummy.strat.param.

The rest of the model is run as before, although we increase the simulation count to three to demonstrate the parameter stochasticity.

control <- control.net(type = "SI", nsims = 3, nsteps = 5, verbose = FALSE)
mod <- netsim(est, param, init, control)

mod

After running 3 simulations we can see that 2 parameters are still displayed under the "Fixed Parameters" section: inf.prob and groups. The other ones are displayed under the "Random Parameters". act.rate and dummy.param now have 3 values associated with them, one per simulation. dummy.strat.param have <list> as value because each simulation has a vector of size 2.

We can inspect the values with the get_param_set function:

``` {r generators_inspect} get_param_set(mod)

#### Parameter set

The drawback of generator function approach above is that it cannot produce correlated parameters. For instance, we may want `dummy.param` and `dummy.strat.param` to be related to one another within a single simulation. We might also use another, external approach for generating parameter sets (e.g., Latin hypercube sampling of multiple parameters).

For this, we need to pre-define a `data.frame` of the possible values:

```r
n <- 5

related.param <- data.frame(
  dummy.param = rbeta(n, 1, 2)
)

related.param$dummy.strat.param_1 <- related.param$dummy.param + rnorm(n)
related.param$dummy.strat.param_2 <- related.param$dummy.param * 2 + rnorm(n)

related.param

We now have a data.frame with 5 rows and 3 columns. Each row contains a set of parameters values that will be used together in a model. This way we keep the relationship between each value.

The first column of the data.frame is named dummy.param, similar to the name of the parameter. For dummy.start.param we need two columns as the parameter contains two values. To achieve this, we name the 2 columns dummy.start.param_1 and dummy.strat.param_2. The value after the underscore informs EpiModel how it should combine these values. This in turn means that the underscore symbol is not allowed in proper parameter names.

Then we set up the rest of the parameters. related.param is saved in the my.randoms list under the special reserved name param_random_set.

my.randoms <- list(
  act.rate = param_random(c(0.25, 0.5, 0.75)),
  param_random_set = related.param

)

param <- param.net(
  inf.prob = 0.3,
  random.params = my.randoms
)

param

Notice here that we combined a generator function for act.rate and a set of correlated parameters with param_random_set.

control <- control.net(type = "SI", nsims = 3, nsteps = 5, verbose = FALSE)
mod <- netsim(est, param, init, control)

mod

The output is similar to what we saw with the generator functions.

Time-Varying Parameters

When doing modeling work, the researcher may want to simulate interventions or events that vary over time. In (our 2021 paper)[https://academic.oup.com/jid/article/223/6/1019/6122459], we simulated the effects of COVID-19–Related behavior change and clinical service interruption on HIV and STI incidence. These included simulations with systemic events with a beginning and an end. To work with these situations, EpiModel offers an optional module to programmatically handle parameter changes across time steps.

This vignette assumes a basic understanding of EpiModel custom models. If the code does not make sense, please check the EpiModel Tutorials.

Parameter Updaters

To define what parameters should change and when during the simulation, we need to define updaters. An updater is a list with to named elements: at, the time step when the change will take place, and param a named list of parameters to update.

list(
  at = 10,
  param = list(
    inf.prob = 0.3,
    act.rate = 0.5
  )
)

This example 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 during the 10th time step.

As mentioned above, we usually want to define several changes at once. EpiModel accept 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
    )
  )
)

 # The `list.of.updaters` goes into `param.net` under `param.updater.list`
 param <- param.net(
   inf.prob = 0.1,
   act.rate = 0.1,
   param.updater.list = list.of.updaters
 )

In this example, we define 2 updaters, one that occurs at time step 100 and the other one at time step 20. In practice, inf.prob and act.rate are set to 0.1 at the beginning by param.net, at time step 100 they are both updated to 0.3 and at time step 125 act.rate is reduced to 0.01.

Enabling the updater Module

This functionality is part of an optional module provided by EpiModel, updater.net. This module must then be enabled in control.net.

Below we set up a complete example with a closed population SI model using the parameters and updaters defined above. We then plot the size of the infected and susceptible population over time to see the effects of the updaters.

 control <- control.net(
   type = NULL, # must be NULL as we use a custom module
   nsims = 1,
   nsteps = 200,
   verbose = FALSE,
   updater.FUN = updater.net,
   infection.FUN = infection.net
 )

nw <- network_initialize(n = 50)
nw <- set_vertex_attribute(nw, "race", rbinom(50, 1, 0.5))
est <- netest(
  nw,
  formation = ~edges,
  target.stats = 25,
  coef.diss = dissolution_coefs(~offset(edges), 10, 0),
  verbose = FALSE
)

init <- init.net(i.num = 10)
mod <- netsim(est, param, init, control)

plot(mod)

This is the minimal example to simplify the demonstration of these time-varying parameters.

Important: when enabling the updater.net module, the researcher must pay particular attention to the order the modules are run. It is recommended to have the updater.net module run first so the changes happen at the beginning of the time step.

Some functionality of updater.net were not described above and are explained in ?EpiModel::updater.net.

Verbosity

When creating an updater, one can add an optional verbose element to the list. If TRUE, updater.net will output a message describing what changes where performed when it occurs during the simulation.

list(
  at = 10,
  param = list(
    inf.prob = 0.3,
    act.rate = 0.5
  ),
  verbose = TRUE
)

Relative Parameter Changes

It is sometime useful to configure the changes with respect to the current value instead of a fixed new value. This is possible, as demonstrated below.

list(
  at = 10,
  param = list(
    inf.prob = function(x) plogis(qlogis(x) + log(2)),
    act.rate = 0.5
  )
)

This updater will set the value of act.rate to 0.5 as we saw earlier. But, for inf.prob we put a function (not a function call). In this case, updater.net 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 an Odds Ratio of 2 to the original value plogis(qlogis(0.1) + log(2)): 0.1818182.

(Advanced) Time-Varying Control

Similarly to time-varying parameters we can use time-varying controls. They work in the same way and use the same updater module. The two differences are:

1. each updater will have a control sub-list instead of a param sub-list
2. the list of updaters goes into control$control.updater.list

# 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
    )
  )
)

 # The `list.of.updaters` goes into `control.net` under `control.updater.list`
 control <- control.net(
   type = NULL, # must be NULL as we use a custom module
   nsims = 1,
   nsteps = 200,
   verbose = TRUE,
   resimulate.network = TRUE,
   updater.FUN = updater.net,
   infection.FUN = infection.net,
   control.updater.list = list.of.updaters
 )

This example set 2 updaters, one that will turn off network resimulation at timestep 100 and another toggling of the verbosity at step 125.

EpiModel documentation built on Feb. 2, 2022, 9:06 a.m.