Processing Intervention Results

Who this is for: anyone turning intervention scenario runs into results tables. Before you start: scenarios defined as in Programmatically Defining Scenarios and run with R/D-interventions/2-scenarios_run.R (locally) or workflow-intervention.R (on the HPC).

The EpiModelHIV-Template's chapter D contains a complete processing pipeline: outcomes.R computes outcomes, labels.R names and formats them, and 3-process_tables.R runs both and writes the table. These files are scaffolds: you adapt them to your research question. This page explains the ideas behind them so you can change them with confidence. The excerpts are from the Template's default branch; your project's copies are the ones to edit.

Where the data come from

Each scenario run saves one file per batch of simulations, sim__<scenario>__<batch>.rds. EpiModelHPC::merge_netsim_scenarios_tibble() then merges them into one tibble per scenario, df__<scenario>.rds in data/run/scenarios/merged_tibbles/, keeping only the last steps_to_keep time steps (in the Template, the intervention period). The HPC workflow runs this merge as one of its steps, and 2-scenarios_run.R does it locally.

Each merged tibble has one row per simulation per time step, with a column for every epidemic tracker (num.B, hiv.inf.B, hiv.incid.B, and so on). A simulation is identified by sim, which is unique within a scenario across batches (the original within-batch number is kept as sim_number).

EpiModelHPC::get_scenarios_tibble_infos() lists the merged files with their scenario names:

scenarios_tibble_dir <- fs::path(scenarios_dir, "merged_tibbles")
scenarios_info <- EpiModelHPC::get_scenarios_tibble_infos(scenarios_tibble_dir)

From time series to one row per simulation

Results tables summarize each simulation by a few numbers, so the first goal is a tibble with one row per simulation and one column per outcome:

scenario_name	sim	lst_ir100_B	cml_incid_B	...
test_1_treat_1	1	value	value	...
test_1_treat_1	2	value	value	...
...	...	...	...	...

We compute outcomes after the run, from the stored trackers, rather than inside the model. That keeps the model's output small and lets you define new outcomes without rerunning anything.

outcomes.R follows a naming convention that the rest of the pipeline relies on:

• lst_ outcomes are averaged over the last year of the run, for rates and prevalences at the end of the intervention period.
• cml_ outcomes are summed over the intervention period, for counts such as cumulative infections.
• _nia and _pia are the number and proportion of infections averted compared with the reference scenario.
• _B, _H, _W mark the race group.

mutate_outcomes() adds the per-time-step outcome columns. For example, HIV incidence per 100 person-years among Black MSM divides new infections by the number at risk:

ir100_m <- 100 * year_steps

mutate_outcomes <- function(d) {
  d |>
    mutate(
      lst_ir100_B = hiv.incid.B / (num.B - hiv.inf.B - hiv.incid.B) * ir100_m,
      cml_incid_B = hiv.incid.B
      # ... and the same for H and W
    )
}

Two functions then collapse the time series, one per prefix:

make_last_year_outcomes <- function(d) {
  d |>
    filter(time >= max(time) - year_steps) |>
    summarise(
      across(starts_with("lst_"), \(x) mean(x, na.rm = TRUE)),
      .by = c("scenario_name", "sim")
    )
}

make_cumulative_outcomes <- function(d) {
  d |>
    filter(time >= intervention_start) |>
    summarise(
      across(starts_with("cml_"), \(x) sum(x, na.rm = TRUE)),
      .by = c("scenario_name", "sim")
    )
}

Because they select columns by prefix, adding an outcome usually means adding one line to mutate_outcomes() with the right prefix.

To add an outcome, check the tracker names in your own output first (glimpse(readRDS(scenarios_info$file_path[1]))): they change between versions of EpiModelHIV, and code written for an older version (for example with i.num or incid.gc) fails on current output.

Infections averted

The number of infections averted (NIA) compares each simulation of an intervention scenario with the reference scenario, usually the one without the intervention. make_d_ref() reads the reference scenario's merged file and summarizes its cml_ outcomes as the median across simulations; mutate_nia_pia() then computes, for each simulation, NIA as the reference value minus the simulation's value, and the proportion of infections averted (PIA) as NIA divided by the reference value.

In 3-process_tables.R, the reference is named directly:

d_ref <- make_d_ref(
  fs::path(scenarios_tibble_dir, "df__test_1_treat_1.rds")
)

When you rename your scenarios, update this file name to your reference scenario. A wrong name is one of the most common mistakes in this chapter.

process_one_scenario() puts the steps together for one scenario: read the file, add outcomes, collapse to one row per simulation, join the lst_ and cml_ results, and add NIA and PIA.

Processing many scenarios in parallel

3-process_tables.R applies process_one_scenario() to every scenario and stacks the results:

d_ls <- future.apply::future_lapply(
  seq_len(nrow(scenarios_info)),
  function(i) process_one_scenario(scenarios_info[i, ], d_ref)
)
d_sc_raw <- dplyr::bind_rows(d_ls)

future_lapply() runs sequentially unless you choose a parallel plan, for example future::plan("multisession", workers = 4) before the call. With many large scenario files, running this step on the HPC as part of the workflow is faster than copying the merged files to your computer.

Formatted tables

The final table has one row per scenario and one column per outcome. Each cell shows the median across simulations followed by the 95% simulation interval (the 2.5th and 97.5th percentiles), for example 4.39 (4.29, 4.52).

labels.R controls the table with two objects:

• var_labels maps outcome names to column labels. Only outcomes listed here appear in the table, in this order:

  var_labels <- c(
    "lst_ir100_B" = "HIV IR100 Black (ly)",
    "cml_incid_B" = "HIV Cumulative Incidence Black (10y)",
    "cml_pia_B"   = "HIV PIA Black (10y)"
    # ...
  )

• format_patterns chooses a number format for each outcome by matching regular expressions against its name. The first pattern that matches wins, so specific patterns come first and a catch-all default comes last. The formatters come from the scales package:

  format_patterns <- list(
    small_num = list(patterns = c("lst_ir100", "cml_nnt"), fun = scales::label_number(0.01)),
    perc      = list(patterns = c("cml_pia"), fun = scales::label_percent(0.1)),
    default   = list(patterns = ".*", fun = scales::label_number(1))
  )

format_table(), defined in the same file, computes the quantiles by scenario, formats each cell, and orders the columns as in var_labels. 3-process_tables.R writes the result:

source("R/D-interventions/labels.R", local = TRUE)

format_table(d_sc_raw, var_labels, format_patterns) |>
  write.csv(fs::path(output_dir, "table.csv"), row.names = FALSE)

Tips

• Several tables. format_table() keeps only the outcomes in var_labels, so to make several tables (for example one for incidence and one for averted infections), call it once per subset of var_labels.
• Publication format. The CSV is for checking results, not for publication. A good way to bridge the gap is to build a formatted table skeleton once in Excel or Google Sheets and paste the values from the CSV into it, so the formatting survives each update of the results.
• Stale results. Delete the old contents of data/run/scenarios/ before rerunning scenarios, or old and new simulations end up mixed in the merged tibbles.
• The time column changed in EpiModel 2.7.0. In a run restarted from a restart point, EpiModel 2.6.x labeled each time step one step late, counting from start (2) rather than 1. EpiModel 2.7.0 labels them from 1, so time is the step the modules saw as at. A filter on an absolute step, such as time >= intervention_start in make_cumulative_outcomes(), therefore also kept the last step before the intervention on 2.6.x; filters relative to max(time) are unaffected. Expect small differences in cumulative outcomes when you compare results from before and after the upgrade.
• Plots. 4-process_plots.R reuses the same outcomes for figures; process_one_scenario_plots() in outcomes.R shows how to recover numeric settings from scenario names.