Programmatically Defining Scenarios

Who this is for: anyone setting up intervention scenarios for a project. Before you start: a project built from the EpiModelHIV-Template (Getting Started with EpiModelHIV), ideally with a calibrated restart point from chapter C (Making a Restart Point).

What a scenario is

A scenario is a named set of parameter changes, applied at a given time step. EpiModel's scenario API describes scenarios as a table with one row per change:

• .scenario.id: the name of the scenario.
• .at: the time step at which the changes take effect. Rows with .at below 2 change the parameters before the simulation starts; rows with .at of 2 or more become parameter updaters that take effect at that time step of the run.
• One column per parameter to change. A vector parameter takes one column per position, named _1, _2, _3 (for race-stratified parameters, 1 = Black, 2 = Hispanic, 3 = White).

Several rows with the same .scenario.id make a time-varying scenario: each row applies its changes at its own .at.

EpiModel::create_scenario_list() turns such a table into a list of scenarios, and EpiModel::use_scenario(param, scenario) applies one scenario to a param object. The EpiModel parameters vignette describes the API in full.

How the Template uses scenarios

In a project built from the EpiModelHIV-Template, scenarios live in chapter D:

• R/D-interventions/make_scenarios.R builds an example 3 by 3 grid (HIV testing rate by treatment initiation rate, each shifted by an odds ratio with apply_or()), starting at intervention_start, and writes it to data/input/scenarios.csv.
• R/D-interventions/2-scenarios_run.R (on your computer) and R/D-interventions/workflow-intervention.R (on the HPC) read that file, convert it with create_scenario_list(), and run every scenario from the restart point with EpiModelHPC::netsim_scenarios().

Every parameter named in a scenario must exist in param; otherwise use_scenario() stops with "Some of the scenario parameters are not defined in param".

A single CSV works well when every scenario changes the same parameters. When different scenarios change different parameters, or you want dozens of scenarios, build them in code as below and either use the list directly or write the CSV from code.

Building many scenarios in code

This example runs in a Template project on the default branch of EpiModelHIV-p. It uses prep.start.rate and tx.init.rate, which exist in both the current (v3.3.x) and the upcoming (v4.0.x) base model. We keep each group of scenarios in its own tibble, because each group changes a different set of parameters.

library(EpiModelHIV)
library(dplyr)

source("R/shared_variables.R", local = TRUE)
source("R/D-interventions/z-context.R", local = TRUE)
source("R/netsim_settings.R", local = TRUE)   # creates `param`

A baseline with no changes, which every analysis needs as its reference:

sc_base <- tibble(.scenario.id = "0_base", .at = 1)

A one-parameter sweep: multiply the PrEP start rate in all three race groups by 0.5, 1.5, and 2 from the start of the intervention period:

mult <- c(0.5, 1.5, 2)
sc_prep <- tibble(
  .scenario.id      = paste0("1_prep_x", mult),
  .at               = intervention_start,
  prep.start.rate_1 = param$prep.start.rate[1] * mult,
  prep.start.rate_2 = param$prep.start.rate[2] * mult,
  prep.start.rate_3 = param$prep.start.rate[3] * mult
)

Passing vectors to tibble() makes one row, and so one scenario, per multiplier. The IDs say what each scenario contains, which makes the results easy to label later.

A two-parameter grid, with every combination of two multipliers:

grid <- tidyr::expand_grid(prep = c(1, 2), tx = c(1, 1.5))
sc_grid <- tibble(
  .scenario.id      = paste0("2_prep_x", grid$prep, "_tx_x", grid$tx),
  .at               = intervention_start,
  prep.start.rate_1 = param$prep.start.rate[1] * grid$prep,
  prep.start.rate_2 = param$prep.start.rate[2] * grid$prep,
  prep.start.rate_3 = param$prep.start.rate[3] * grid$prep,
  tx.init.rate_1    = pmin(1, param$tx.init.rate[1] * grid$tx),
  tx.init.rate_2    = pmin(1, param$tx.init.rate[2] * grid$tx),
  tx.init.rate_3    = pmin(1, param$tx.init.rate[3] * grid$tx)
)

pmin(1, ...) keeps the scaled per-step probabilities at or below 1.

A time-limited intervention: two rows share one ID, so the PrEP start rate doubles at intervention_start and returns to the baseline two years later:

sc_pulse <- tibble(
  .scenario.id      = "3_prep_pulse",
  .at               = c(intervention_start, intervention_start + 2 * year_steps),
  prep.start.rate_1 = c(param$prep.start.rate[1] * 2, param$prep.start.rate[1]),
  prep.start.rate_2 = c(param$prep.start.rate[2] * 2, param$prep.start.rate[2]),
  prep.start.rate_3 = c(param$prep.start.rate[3] * 2, param$prep.start.rate[3])
)

Now convert each tibble to a scenario list separately and join the lists with purrr::reduce():

sc_tibbles <- list(sc_base, sc_prep, sc_grid, sc_pulse)
scenarios_list <- purrr::reduce(
  sc_tibbles,
  \(out, d) c(out, EpiModel::create_scenario_list(d)),
  .init = list()
)
names(scenarios_list)
#> [1] "0_base"            "1_prep_x0.5"       "1_prep_x1.5"
#> [4] "1_prep_x2"         "2_prep_x1_tx_x1"   "2_prep_x1_tx_x1.5"
#> [7] "2_prep_x2_tx_x1"   "2_prep_x2_tx_x1.5" "3_prep_pulse"

Converting the tibbles one at a time matters; see the pitfalls below.

Before you send anything to the HPC, check every scenario against your param object. use_scenario() stops with an error for an unknown parameter name:

for (sc in scenarios_list) invisible(EpiModel::use_scenario(param, sc))

scenarios_list can go straight to EpiModelHPC::netsim_scenarios(). To use the Template's scripts unchanged instead, write the scenarios to data/input/scenarios.csv, but only combine tibbles that set the same parameters:

sc_csv <- bind_rows(sc_base, sc_prep, sc_grid)   # sc_base sets nothing, so it is safe

Pitfalls

• On EpiModel 2.6.x, set every position of a vector parameter. Before EpiModel 2.7.0, a scenario that sets only prep.start.rate_1 = 0.5 turns prep.start.rate from c(0.0055, 0.0042, 0.0066) into the single value 0.5. The run then crashes later, for example with "Some posit_ids are larger than the number of nodes" in the PrEP module. EpiModel 2.7.0 fixes this: a _1 column changes only the first position and leaves the others at their baseline values. Check the EpiModel version in your project's renv.lock; on 2.6.x, set _1, _2, and _3 together, as in the examples above, even when one of them keeps its baseline value.
• Do not bind_rows() tibbles with different columns before create_scenario_list(). bind_rows() fills the missing columns with NA, and the scenario then sets those parameters to NA. For example, binding sc_prep with sc_grid gives the PrEP-only scenarios tx.init.rate = NA NA NA. Convert each tibble separately, as above.
• .at is an absolute time step, not a time since the intervention started. In chapter D runs start from the restart point, so use intervention_start (defined in R/shared_variables.R), or intervention_start + k for later changes.
• Do not use __ (two underscores) in scenario IDs. Simulation files are named sim__<scenario>__<batch>.rds and merged files df__<scenario>.rds, so a double underscore in the ID breaks the file names. Stick to letters, digits, _, ., and -.

Advanced example: the DoxyPEP project

This longer example comes from our DoxyPEP project, which implements STI prophylaxis with doxycycline. Its parameters exist only on that project's EpiModelHIV-p branch, not on the default branch, so the code does not run in a fresh Template project; read it for the techniques: groups of scenarios with different parameters, placebo scenarios, and a two-way grid for contour plots. Without any scenario, all DoxyPEP interventions are off.

Setup

library(EpiModelHIV)
library(dplyr)

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

The baseline never starts the intervention, by setting its start after the end of the run. We store each group of scenarios as a tibble in a list, and combine the lists at the end.

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

Eligibility scenarios

The intervention scenarios vary who is eligible for DoxyPEP: PrEP users (prep), people with diagnosed HIV (poz), and people with an STI diagnosis in the past 6 or 12 months (sti). The base rates of starting DoxyPEP for each group were calibrated so that 10% of the eligible population is on DoxyPEP; each scenario applies a multiplier to them.

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

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

sc_interv_list <- list()

# PrEP users only
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
)

# People with diagnosed HIV only
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
)

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

# One STI diagnosis in the past 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
)

# PrEP or diagnosed HIV, AND one STI diagnosis in the past 12 months (52 weeks)
sc_interv_list[["prep_poz_and_1sti_12m_sc"]] <- tibble(
  .scenario.id = paste0("interv_5_prep_poz_and_1sti_12m_", 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
)

In the STI scenarios, sti.dx.threshold is the number of diagnoses within the window needed to be eligible, and gc.dx.window and ct.dx.window set the window in weeks. In the last group, doxy.indic.combined = TRUE means a person needs the STI criterion and one of the other two, and the STI start probability stays at its base value while the other two are scaled.

The five tibbles have different columns, so we convert each one separately and join the results:

sc_df_list <- c(sc_no_list, sc_interv_list)

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

length(scenarios_list)
#> [1] 46
head(names(scenarios_list), 3)
#> [1] "0_no_doxy"          "interv_1_prep_0.5"  "interv_1_prep_0.75"
tail(names(scenarios_list), 2)
#> [1] "interv_5_prep_poz_and_1sti_12m_2.5" "interv_5_prep_poz_and_1sti_12m_3"

Placebo scenarios

The project also needed placebo scenarios, in which people start DoxyPEP but DoxyPEP has no effect, one for each intervention scenario. DoxyPEP acts through reduced STI transmission, additional STI screening, and STI treatment when screened positive. This tibble sets each of those to its no-effect value:

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
)

A small function adds these columns to every intervention tibble and renames the scenarios:

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, placebo_sc)
  }
  sc_placebo_list
}

sc_interv_placebo_list <- make_placebo_list(sc_interv_list, sc_placebo)

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

Scenarios for a contour plot

To plot gonorrhea incidence as a function of two parameters, we need every combination of them: doxy.eff.gc.or, the odds ratio of gonorrhea infection on DoxyPEP, and sti.screen.doxy.rate, the STI screening rate for people on DoxyPEP. We define the screening interval in months and convert it to a weekly rate.

or_sc <- c(0.01, seq(0.1, 1, 0.1))
screen_int <- c(seq(1, 12, 1), Inf)
s <- expand.grid(or = or_sc, screen = screen_int)

sc_contour_plots <- tibble(
  .scenario.id = paste0("cp_or_", s$or, "_scr_", s$screen),
  .at = 1,
  # varied for the contour plot
  doxy.eff.gc.or = s$or,
  sti.screen.doxy.rate = 1 / (s$screen * 52 / 12), # months to a weekly rate
  # shared by all scenarios
  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
)

sc_contour_plots_list <- EpiModel::create_scenario_list(sc_contour_plots)

All 143 scenarios set the same parameters, so a single tibble is enough. The odds ratio and screening interval are written into each .scenario.id, so they can be recovered when making the plot:

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

Next: Processing Intervention Results.