Programatically Defining Scenarios - EpiModel/EpiModeling GitHub Wiki

Introduction

This tutorial explain how to programmatically define many EpiModel scenarios using the scenario API.

For this example, we will use the DoxyPEP project that implements STI prophylaxis with doxycycline. Without any scenarios, we assume that all doxy related interventions are turned off.

We want here to define multiple scenarios that toggle on and off several elements of the doxy intervention as well as assess the effects of different intensity of intervention.

For this tutorial we do not define all the scenarios as a single "scenarios.csv" files for three reason:

Writing a 200 row CSV file by hand is a pain
More than 12 parameters in total are affected by the different scenarios. But most scenarios only touch a few at once.
We also want placebo scenarios were the interventions are turned on but do nothing. (For tracking purposes)

Setup

We just load the libraries we need and define the intervention_start and intervention_end variable. It's easier to set them here, especially if at some point these values need to change.

library("EpiModelHIV")
library("dplyr")

intervention_start <- 52 * 60
intervention_end <- intervention_start + 52 * 10

The baseline / no intervention scenario

We usually need a business as usual scenario. In this case it's simply a scenario were we never start the intervention.

sc_no_list <- list()
sc_no_list[["no_doxy_sc"]] <- tibble(
  .scenario.id = "0_no_doxy",
  .at = 1,
  doxy.start = intervention_end + 1
)

All over this tutorial, we will store the scenarios as tibbles within lists. These list will be combined at the end.

Here we have only one scenario in this "baseline" sc_no_list. But we still make it a list for consistency reasons. Also, we may need different baselines in the future.

The sc_no_list contains on element called no_doxy_sc which is a tibble formatted according to the scenario API:

.at is when should the scenario apply: 1 means "at the start of simulation"
.scenario.id is the identifier for the scenario: "0_no_doxy"
doxy.start is the parameter we want changed: we set the start of intervention after the end of simulation. (never)

This scenario is not very interesting but the process will be shared by all the others.

Intervention scenarios

For this next set of scenarios, we will vary who is eligible to DoxyPEP:

PrEP users (prep)
HIV positive individuals (poz)
Individuals with one STI diagnosis in the past 6 month (sti1_6m)

Setup

We start by storing the default parameters for our 3 groups of eligibles. These are the rate of starting DoxyPEP for an individual eligible to any of these path. These parameters were calibrated to have 10% of the eligible population on DoxyPEP.

prep_prob_base <- 0.003
poz_prob_base <- 0.0012
sti1_6_base <- 0.0068

Our scenarios will consist in turning these mechanism on and to apply a modifier to these base rates.

We'll always use the same modifiers: from 0.5 to 3.

prop_incr <- c(0.5, 0.75, 1, 1.1, 1.2, 1.5, 2, 2.5, 3)

Finally, we create a list to store the scenarios.

sc_interv_list <- list()

The PrEP eligibility scenarios

We begin with the scenarios where DoxyPEP eligibility is only for PrEP users.

Two parameters need updating:

doxy.start: this time we'll start the intervention at time step intervention_start
doxy.start.prep.prob: the probability of starting DoxyPEP if eligible through PrEp.

sc_interv_list[["prep_sc"]] <- tibble(
  .scenario.id = paste0("interv_1_prep_", prop_incr),
  .at = 1,
  doxy.start = intervention_start,
  doxy.start.prep.prob = prep_prob_base * prop_incr
)

We create a tibble that we store inside the sc_interv_list list under the name "prep_sc".

Note that for .scenario.id and doxy.start.prep.prob we pass vectors to the tibble call.

paste0("interv_1_prep_", prop_incr)
#> [1] "interv_1_prep_0.5"  "interv_1_prep_0.75" "interv_1_prep_1"
#> [4] "interv_1_prep_1.1"  "interv_1_prep_1.2"  "interv_1_prep_1.5"
#> [7] "interv_1_prep_2"    "interv_1_prep_2.5"  "interv_1_prep_3"

prep_prob_base * prop_incr
#> [1] 0.00150 0.00225 0.00300 0.00330 0.00360 0.00450 0.00600 0.00750 0.00900

This creates 9 scenarios, with increasing values of doxy.start.prep.prob. The scenarios IDs are chosen to make their content easily identifiable.

The `poz` eligibility scenarios

Similarly, we define the scenarios for eligibility through being HIV positive.

It's almost identical, except that we modify the doxy.start.poz.prob parameters and not the doxy.start.prep.prob one.

Again, 9 scenarios are created

sc_interv_list[["poz_sc"]] <- tibble(
  .scenario.id = paste0("interv_2_poz_", prop_incr),
  .at = 1,
  doxy.start = intervention_start,
  doxy.start.poz.prob = poz_prob_base * prop_incr
)

`prep` and `poz` together

This time we make both eligibility criteria active at once and create the 9 corresponding scenarios.

sc_interv_list[["prep_poz_sc"]] <- tibble(
  .scenario.id = paste0("interv_3_prep_poz_", prop_incr),
  .at = 1,
  doxy.start = intervention_start,
  doxy.start.prep.prob = prep_prob_base * prop_incr,
  doxy.start.poz.prob = poz_prob_base * prop_incr
)

The STI diagnosis eligibility criteria

For the STIs, in addition to the doxy.start.sti.prob, we need to define 3 other parameters:

sti.dx.threshold, the number of diagnosis during the time window to be eligible to DoxyPEP.
gc.dx.window and ct.dx.window which governs the window of diagnosis we are looking at, here 6 months (26 weeks).

sc_interv_list[["1sti_6m_sc"]] <- tibble(
  .scenario.id = paste0("interv_4_1sti_6m_", prop_incr),
  .at = 1,
  doxy.start = intervention_start,
  gc.dx.window = 26,
  ct.dx.window = 26,
  sti.dx.threshold = 1,
  doxy.start.sti.prob = sti1_6_base * prop_incr
)

All combined

This final set of 9 scenarios combined all of the above plus a doxy.indic.combined that is set to TRUE to indicate that a node needs the STI criterion AND either prep or poz to start DoxyPEP.

sc_interv_list[["prep_poz_scAND1sti6_sc"]] <- tibble(
  .scenario.id = paste0("interv_5_prep_poz_scAND1sti6_", prop_incr),
  .at = 1,
  doxy.start = intervention_start,
  doxy.start.prep.prob = prep_prob_base * prop_incr,
  doxy.start.poz.prob = poz_prob_base * prop_incr,
  doxy.indic.combined = TRUE,
  gc.dx.window = 52,
  ct.dx.window = 52,
  sti.dx.threshold = 1,
  doxy.start.sti.prob = sti1_6_base
)

The final list

Below is the content of the sc_interv_list list. It has 5 elements which are all tibbles with different columns.

sc_interv_list
#> $prep_sc
#> # A tibble: 9 × 4
#>   .scenario.id         .at doxy.start doxy.start.prep.prob
#>   <chr>              <dbl>      <dbl>                <dbl>
#> 1 interv_1_prep_0.5      1       3120              0.0015
#> ...
#>
#> $poz_sc
#> # A tibble: 9 × 4
#>   .scenario.id        .at doxy.start doxy.start.poz.prob
#>   <chr>             <dbl>      <dbl>               <dbl>
#> 1 interv_2_poz_0.5      1       3120             0.0006
#> ...
#>
#> $prep_poz_sc
#> # A tibble: 9 × 5
#>   .scenario.id             .at doxy.start doxy.start.prep.prob doxy.start.po
#>   <chr>                  <dbl>      <dbl>                <dbl>
#> 1 interv_3_prep_poz_0.5      1       3120              0.0015              0
#> ...
#>
#> $`1sti_6m_sc`
#> # A tibble: 9 × 7
#>   .scenario.id            .at doxy.start gc.dx.window ct.dx.window sti.dx.th
#>   <chr>                 <dbl>      <dbl>        <dbl>        <dbl>
#> 1 interv_4_1sti_6m_0.5      1       3120           26           26
#> ...
#> # ℹ 1 more variable: doxy.start.sti.prob <dbl>
#>
#> $prep_poz_scAND1sti6_sc
#> # A tibble: 9 × 10
#>   .scenario.id              .at doxy.start doxy.start.prep.prob doxy.start.p
#>   <chr>                   <dbl>      <dbl>                <dbl>
#> 1 interv_5_prep_poz_s…     1       3120              0.0015
#> ...
#> # ℹ 5 more variables: doxy.indic.combined <lgl>, gc.dx.window <dbl>,
#> #   ct.dx.window <dbl>, sti.dx.threshold <dbl>, doxy.start.sti.prob <dbl>

Creation of the scenarios list

With the scenario API, we usually call create_scenario_list to go from a tibble of scenarios to a list usable by EpiModel. However, here we have lists of tibbles.

The fix is pretty simple:

first we combine all the lists of tibbles into a single one

sc_df_list <- c(
  sc_no_list,
  sc_interv_list
)

Now sc_df_list contains all the tibbles, for the baseline as well as the intervention scenarios.

Then we use purrr::reduce to turn each tibble into a list of scenarios and to combine all in a single list.

scenarios_list <- purrr::reduce(
  sc_df_list,
  \(out, d_sc) c(out, EpiModel::create_scenario_list(d_sc)),
  .init = list()
)

names(scenarios_list)
#>  [1] "0_no_doxy"                     "interv_1_prep_0.5"
#>  [3] "interv_1_prep_0.75"            "interv_1_prep_1"
#>  [5] "interv_1_prep_1.1"             "interv_1_prep_1.2"
#>  [7] "interv_1_prep_1.5"             "interv_1_prep_2"
#>  [9] "interv_1_prep_2.5"             "interv_1_prep_3"
#> [11] "interv_2_poz_0.5"              "interv_2_poz_0.75"
#> [13] "interv_2_poz_1"                "interv_2_poz_1.1"
#> [15] "interv_2_poz_1.2"              "interv_2_poz_1.5"
#> [17] "interv_2_poz_2"                "interv_2_poz_2.5"
#> [19] "interv_2_poz_3"                "interv_3_prep_poz_0.5"
#> [21] "interv_3_prep_poz_0.75"        "interv_3_prep_poz_1"
#> [23] "interv_3_prep_poz_1.1"         "interv_3_prep_poz_1.2"
#> [25] "interv_3_prep_poz_1.5"         "interv_3_prep_poz_2"
#> [27] "interv_3_prep_poz_2.5"         "interv_3_prep_poz_3"
#> [29] "interv_4_1sti_6m_0.5"          "interv_4_1sti_6m_0.75"
#> [31] "interv_4_1sti_6m_1"            "interv_4_1sti_6m_1.1"
#> [33] "interv_4_1sti_6m_1.2"          "interv_4_1sti_6m_1.5"
#> [35] "interv_4_1sti_6m_2"            "interv_4_1sti_6m_2.5"
#> [37] "interv_4_1sti_6m_3"            "t2_5_prep_poz_scAND1sti6_0.5"
#> [39] "t2_5_prep_poz_scAND1sti6_0.75" "t2_5_prep_poz_scAND1sti6_1"
#> [41] "t2_5_prep_poz_scAND1sti6_1.1"  "t2_5_prep_poz_scAND1sti6_1.2"
#> [43] "t2_5_prep_poz_scAND1sti6_1.5"  "t2_5_prep_poz_scAND1sti6_2"
#> [45] "t2_5_prep_poz_scAND1sti6_2.5"  "t2_5_prep_poz_scAND1sti6_3"

This leaves us with a list of 46 scenarios that can be passed to EpiModelHPC::netsim_scenarios

Placebo scenarios

For this project, we needed scenarios where individuals could start DoxyPEP but where DoxyPEP would have no effect whatsoever. We call these placebo scenarios and needed one for each of the intervention scenarios defined above.

DoxyPEP effects come from several mechanisms:

STI transmission reduction
Additional STI screening
Specific STI treatment probability if screened positive

Below we make a tibble that set the relevant parameters to their no effect values.

sc_placebo <- tibble(
  sti.screen.doxy.start = intervention_end + 1,
  sti.doxy.tx.prob = 0,
  doxy.eff.gc.or = 1,
  doxy.eff.ct.or = 1
)

Then we make a function to add these columns to the scenarios we want make into placebo ones.

make_placebo_list <- function(sc_list, placebo_sc) {
  sc_placebo_list <- list()
  for (sc_name in names(sc_list)) {
    sc_d <- sc_list[[sc_name]]
    sc_d[[".scenario.id"]] <- paste0(sc_d[[".scenario.id"]], "_placebo")
    sc_placebo_list[[paste0(sc_name, "_placebo")]] <-
      bind_cols(sc_d, sc_placebo)
  }
  sc_placebo_list
}

Finally we apply this function to our sc_interv_list to create the sc_interv_placebo_list.

sc_interv_placebo_list <- make_placebo_list(sc_interv_list, sc_placebo)

Now we can reuse the code to create the full list of scenarios including the placebo ones.

sc_df_list <- c(
  sc_no_list,
  sc_interv_list,
  sc_interv_placebo_list
)

scenarios_list <- purrr::reduce(
  sc_df_list,
  \(out, d_sc) c(out, EpiModel::create_scenario_list(d_sc)),
  .init = list()
)

Contour plots scenarios

This section describe how to make the scenarios to graph a contour plot of the impact of DoxyPEP on Gonorrhea (GC) incidence as a function of the effectiveness of DoxyPEP to prevent this infection as well as the STI screening interval for individuals on DoxyPEP.

Values of the parameters

We'll look at two parameters:

doxy.gc.or: the infection risk reduction of doxyPEP on GC
sti.screen.doxy.rate: the STI screening rate for individuals on DoxyPEP

For doxy.gc.or:

or_sc <- c(0.01, seq(0.1, 1, 0.1))
#> [1] 0.01 0.10 0.20 0.30 0.40 0.50 0.60 0.70 0.80 0.90 1.00

For sti.screen.doxy.rate we first define the mean screening interval in months. The conversion to rate will be done later in the process.

screen_int <- c(seq(1, 12, 1), Inf)
#> [1]   1   2   3   4   5   6   7   8   9  10  11  12 Inf

For the contour plot we need all the combinations of these parameters.

s <- expand.grid(or = or_sc, screen = screen_int)

Making the scenarios

This time, all the scenarios share the same paramters. We can make the as a single tibble.

sc_contour_plots <- tibble(
  .scenario.id = paste0("cp_or_", s$or, "_scr_", s$screen),
  .at = 1,
  # contour plot specific
  doxy.eff.gc.or = s$or,
  sti.screen.doxy.rate = 1 / (s$screen * 52 / 12), # convert to a rates
  # shared
  doxy.start = intervention_start,
  gc.dx.window = 52,
  ct.dx.window = 52,
  sti.dx.threshold = 1,
  doxy.start.sti.prob = sti1_6_base,
  doxy.start.prep.prob = prep_prob_base,
  doxy.start.poz.prob = poz_prob_base
)

*Note: The information about the screening interval and odds ratio is embedded in the .scenario.id. This way it will be easier to get it back when making the plot.

The following snippet gives an example of how to process it: *

sc_contour_plots |>
  select(.scenario.id) |>
  tidyr::separate_wider_delim(
    .scenario.id,
    "_",
    names = c(NA, NA, "or", NA, "screen"),
    cols_remove = FALSE
  ) |>
  mutate(
    or = as.numeric(or),
    screen = as.numeric(screen)
  )
#> # A tibble: 143 × 3
#>       or screen .scenario.id
#>    <dbl>  <dbl> <chr>
#>  1  0.01      1 cp_or_0.01_scr_1
#>  2  0.1       1 cp_or_0.1_scr_1
#>  3  0.2       1 cp_or_0.2_scr_1
#>  4  0.3       1 cp_or_0.3_scr_1
#>  5  0.4       1 cp_or_0.4_scr_1
#>  6  0.5       1 cp_or_0.5_scr_1
#>  7  0.6       1 cp_or_0.6_scr_1
#>  8  0.7       1 cp_or_0.7_scr_1
#>  9  0.8       1 cp_or_0.8_scr_1
#> 10  0.9       1 cp_or_0.9_scr_1
#> # ℹ 133 more rows

Creating the scenario list

As we only have one tibble here, we can directly call EpiModel::create_scenario_list

sc_contour_plots_list <- EpiModel::create_scenario_list(sc_contour_plots)

Programatically Defining Scenarios - EpiModel/EpiModeling GitHub Wiki

Introduction

Setup

The baseline / no intervention scenario

Intervention scenarios

Setup

The PrEP eligibility scenarios

The poz eligibility scenarios

prep and poz together

The STI diagnosis eligibility criteria

All combined

The final list

Creation of the scenarios list

Placebo scenarios

Contour plots scenarios

Values of the parameters

Making the scenarios

Creating the scenario list

⚠️ **GitHub.com Fallback** ⚠️

The `poz` eligibility scenarios

`prep` and `poz` together

⚠️ GitHub.com Fallback ⚠️