# Working With Model Parameters

## 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
    )
    #> 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.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
    #> EpiModel Simulation
    #> =======================
    #> Model class: netsim
    #>
    #> Simulation Summary
    #> -----------------------
    #> Model type: SI
    #> No. simulations: 1
    #> No. time steps: 5
    #> No. NW groups: 1
    #>
    #> Fixed Parameters
    #> ---------------------------
    #> inf.prob = 0.3
    #> act.rate = 0.5
    #> dummy.param = 4
    #> dummy.strat.param = 0 1
    #> groups = 1
    #>
    #> Model Output
    #> -----------------------
    #> Variables: s.num i.num num si.flow
    #> Networks: sim1
    #> Transmissions: sim1
    #>
    #> Formation Diagnostics
    #> -----------------------
    #>       Target Sim Mean Pct Diff Sim SD
    #> edges     25     28.8     15.2  2.387
    #>
    #>
    #> Dissolution Diagnostics
    #> -----------------------
    #>                Target Sim Mean Pct Diff Sim SD
    #> Edge Duration    10.0    9.242   -7.578     NA
    #> Pct Edges Diss    0.1    0.072  -28.276     NA

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
    #> Fixed Parameters
    #> ---------------------------
    #> inf.prob = 0.3
    #>
    #> Random Parameters
    #> (Not drawn yet)
    #> ---------------------------
    #> act.rate = <function>
    #> dummy.param = <function>
    #> dummy.strat.param = <function>

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
    #> EpiModel Simulation
    #> =======================
    #> Model class: netsim
    #>
    #> Simulation Summary
    #> -----------------------
    #> Model type: SI
    #> No. simulations: 3
    #> No. time steps: 5
    #> No. NW groups: 1
    #>
    #> Fixed Parameters
    #> ---------------------------
    #> inf.prob = 0.3
    #> groups = 1
    #>
    #> Random Parameters
    #> ---------------------------
    #> act.rate = 0.5 0.25 0.5
    #> dummy.param = 0.02531928 0.251341 0.03892299
    #> dummy.strat.param = <list>
    #>
    #> Model Output
    #> -----------------------
    #> Variables: s.num i.num num si.flow
    #> Networks: sim1 ... sim3
    #> Transmissions: sim1 ... sim3
    #>
    #> Formation Diagnostics
    #> -----------------------
    #>       Target Sim Mean Pct Diff Sim SD
    #> edges     25   25.733    2.933  3.845
    #>
    #>
    #> Dissolution Diagnostics
    #> -----------------------
    #>                Target Sim Mean Pct Diff Sim SD
    #> Edge Duration    10.0    9.318   -6.819  0.085
    #> Pct Edges Diss    0.1    0.075  -24.740  0.023

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:

    get_param_set(mod)
    #>   sim inf.prob vital groups act.rate dummy.param dummy.strat.param_1
    #> 1   1      0.3 FALSE      1     0.50  0.02531928          0.04814017
    #> 2   2      0.3 FALSE      1     0.25  0.25134096          0.04426900
    #> 3   3      0.3 FALSE      1     0.50  0.03892299          0.05038573
    #>   dummy.strat.param_2
    #> 1           0.2378170
    #> 2           0.1604752
    #> 3           0.1277672

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

    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
    #>   dummy.param dummy.strat.param_1 dummy.strat.param_2
    #> 1  0.30583678           0.5932591         0.214550981
    #> 2  0.53237743           1.9773455        -0.009634672
    #> 3  0.14786943           3.4883784         0.442632331
    #> 4  0.07227561          -0.3417346         2.659206204
    #> 5  0.24006401          -0.8567142         0.795942272

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
    #> Fixed Parameters
    #> ---------------------------
    #> inf.prob = 0.3
    #>
    #> Random Parameters
    #> (Not drawn yet)
    #> ---------------------------
    #> act.rate = <function>
    #> param_random_set = <data.frame> ( dimensions: 5 3 )

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
    #> EpiModel Simulation
    #> =======================
    #> Model class: netsim
    #>
    #> Simulation Summary
    #> -----------------------
    #> Model type: SI
    #> No. simulations: 3
    #> No. time steps: 5
    #> No. NW groups: 1
    #>
    #> Fixed Parameters
    #> ---------------------------
    #> inf.prob = 0.3
    #> groups = 1
    #>
    #> Random Parameters
    #> ---------------------------
    #> dummy.param = 0.240064 0.07227561 0.1478694
    #> dummy.strat.param = <list>
    #> act.rate = 0.75 0.5 0.25
    #>
    #> Model Output
    #> -----------------------
    #> Variables: s.num i.num num si.flow
    #> Networks: sim1 ... sim3
    #> Transmissions: sim1 ... sim3
    #>
    #> Formation Diagnostics
    #> -----------------------
    #>       Target Sim Mean Pct Diff Sim SD
    #> edges     25     24.8     -0.8  2.597
    #>
    #>
    #> Dissolution Diagnostics
    #> -----------------------
    #>                Target Sim Mean Pct Diff Sim SD
    #> Edge Duration    10.0    9.214   -7.859  0.088
    #> Pct Edges Diss    0.1    0.079  -20.984  0.031

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
    )
    #> Starting maximum pseudolikelihood estimation (MPLE):
    #> Evaluating the predictor and response matrix.
    #> Maximizing the pseudolikelihood.
    #> Finished MPLE.
    #> Stopping at the initial estimate.

    init <- init.net(i.num = 10)
    mod <- netsim(est, param, init, control)
    #>
    #>  A MESSAGE occured in module 'updater.FUN' at step 100
    #>
    #>
    #> At timestep = 100 the following parameters where modified:
    #> 'inf.prob', 'act.rate'
    #>
    #>  A MESSAGE occured in module 'updater.FUN' at step 125
    #>
    #>
    #> At timestep = 125 the following parameters where modified:
    #> 'inf.prob'

    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.

    #> $at
    #> [1] 10
    #>
    #> $param
    #> $param$inf.prob
    #> [1] 0.3
    #>
    #> $param$act.rate
    #> [1] 0.5
    #>
    #>
    #> $verbose
    #> [1] 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.

    #> $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 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`

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