Compare ebola scenarios

R/scenario_comparison.R

@@ -2,7 +2,7 @@
 #'
 #' @param baseline A nested `<data.table>` of a single model outcome with
 #' parameter uncertainty. This is expected to be the output of a single call to
-#' deterministic model functions such as [model_default()].
+#' model functions such as [model_default()] or [model_ebola()].
 #' @param scenarios A nested `<data.table>` of any number of model outcomes with
 #' infection parameter sets identical to those in `baseline` for comparability,
 #' i.e., the `scenarios` must differ from the `baseline` only in any
@@ -21,11 +21,27 @@
 #' When `summarise = TRUE`, a `<data.table>` of the same number of
 #' rows as `scenarios`, with the columns `"scenario"`, `"averted_median"`,
 #' `"averted_lower"`, and `"averted_upper"`.
+#'
 #' When `summarise = FALSE`, a `<data.table>` with one row per `"scenario"`,
 #' parameter set (`"param_set"`), and demography group (`"demography_group`),
 #' with the additional column `"outcomes_averted"` giving the difference between
 #' the baseline and the comparator scenario for each parameter set for each
 #' demography group.
+#'
+#' @details
+#' Both deterministic and stochastic models (currently only the Ebola model) are
+#' supported.
+#'
+#' When comparing deterministic model scenarios, users are expected to ensure
+#' that outputs comparable in terms of demographic groups and parameters.
+#' The output is expected to have parameter uncertainty, and differences between
+#' each scenario and the baseline are calculated after matching on parameter
+#' sets.
+#'
+#' When comparing stochastic model scenarios, each scenario is matched against
+#' the baseline on the replicate number as well as the parameter set to reduce
+#' the effect of initial conditions on differences in outcomes.
+#'
 #' @export
 #'
 #' @examples
@@ -148,8 +164,8 @@ outcomes_averted <- function(baseline,
         colnames(baseline),
         must.include = c(
           "transmission_rate", "infectiousness_rate", "time_end",
-          "param_set", "population", "intervention", "vaccination",
-          "time_dependence", "increment", "scenario", "data"
+          "param_set", "population", "intervention",
+          "time_dependence", "scenario", "data"
         )
       )
   )
@@ -194,8 +210,15 @@ outcomes_averted <- function(baseline,
     )]
   )
 
+  # set flag for whether data has replicates; applies only to stochastic model
+  has_reps <- if ("replicate" %in% colnames(baseline[["data"]][[1L]])) {
+    "replicate"
+  } else {
+    NULL
+  }
+
   baseline_outcomes <- baseline_outcomes[, unlist(outcome, recursive = FALSE),
-    by = c("scenario", "param_set"),
+    by = c("scenario", "param_set")
   ]
   data.table::setnames(baseline_outcomes, "value", "baseline_value")
 
@@ -215,7 +238,7 @@ outcomes_averted <- function(baseline,
 
   # variables to merge on; # cannot use ifelse()
   merge_variables <- c(
-    "param_set", if (by_group) "demography_group" else NULL
+    "param_set", if (by_group) "demography_group" else NULL, has_reps
   )
 
   # merge and get difference in outcomes

tests/testthat/test-scenario_comparison.R

@@ -94,7 +94,7 @@ scenarios <- model_default(
   intervention = intervention_sets
 )
 
-test_that("`outcomes_averted()`: Basic expectations", {
+test_that("`outcomes_averted()`: Basic expectations, deterministic models", {
   # expect runs without conditions
   expect_no_condition(
     outcomes_averted(
@@ -215,3 +215,176 @@ test_that("`outcomes_averted()`: Errors and messages", {
     regexp = "`scenarios` must have common infection parameter sets"
   )
 })
+
+#### Scenario comparison for stochastic models ####
+# Prepare population and parameters
+demography_vector <- 67000 # small population
+contact_matrix <- matrix(1)
+
+# manual case counts divided by pop size rather than proportions as small sizes
+# introduce errors when converting to counts in the model code; extra
+# individuals may appear
+infectious <- 1
+exposed <- 10
+initial_conditions <- matrix(
+  c(demography_vector - infectious - exposed, exposed, infectious, 0, 0, 0) /
+    demography_vector,
+  nrow = 1
+)
+rownames(contact_matrix) <- "full_pop"
+pop <- population(
+  contact_matrix = contact_matrix,
+  demography_vector = demography_vector,
+  initial_conditions = initial_conditions
+)
+
+compartments <- c(
+  "susceptible", "exposed", "infectious", "hospitalised", "funeral", "removed"
+)
+
+# prepare integer values for expectations on data length
+time_end <- 100L
+replicates <- 10L
+
+# create vector of parameters
+beta <- withr::with_seed(
+  1,
+  rnorm(10, mean = 1.3 / 7, sd = 0.005)
+)
+
+baseline <- withr::with_seed(
+  1,
+  model_ebola(
+    population = pop,
+    transmission_rate = beta,
+    replicates = replicates
+  )
+)
+
+max_time <- 100
+# prepare durations as starting at 25% of the way through an epidemic
+# and ending halfway through
+time_begin <- max_time / 4
+time_end <- max_time / 2
+
+# create two interventions
+reduce_beta <- intervention(
+  type = "rate", time_begin = time_begin, time_end = time_end,
+  reduction = 0.9
+)
+
+improve_etu <- intervention(
+  type = "rate", time_begin = time_begin, time_end = time_end,
+  reduction = 0.5
+)
+
+intervention_sets <- list(
+  list(
+    transmission_rate = reduce_beta
+  ),
+  list(
+    etu_risk = improve_etu
+  ),
+  list(
+    transmission_rate = reduce_beta,
+    etu_risk = improve_etu
+  )
+)
+
+scenarios <- withr::with_seed(
+  1,
+  model_ebola(
+    population = pop,
+    transmission_rate = beta,
+    intervention = intervention_sets,
+    replicates = replicates
+  )
+)
+
+outcomes_averted(
+  baseline = baseline,
+  scenarios = scenarios, summarise = FALSE
+)
+
+test_that("`outcomes_averted()`: Basic expectations, stochastic models", {
+  # expect runs without conditions
+  expect_no_condition(
+    outcomes_averted(
+      baseline = baseline, scenarios = scenarios
+    )
+  )
+  expect_no_condition(
+    outcomes_averted(
+      baseline = baseline, scenarios = scenarios, summarise = FALSE
+    )
+  )
+  expect_no_condition(
+    outcomes_averted(
+      baseline = baseline, scenarios = scenarios, by_group = FALSE
+    )
+  )
+  expect_no_condition(
+    outcomes_averted(
+      baseline = baseline, scenarios = scenarios, by_group = FALSE,
+      summarise = FALSE
+    )
+  )
+
+  # expect return type, structure, and basic statistical correctness
+  # for default options
+  averted_default <- outcomes_averted(
+    baseline = baseline, scenarios = scenarios, by_group = TRUE
+  )
+  expect_s3_class(averted_default, "data.table")
+  expect_named(
+    averted_default,
+    c(
+      "scenario", "demography_group",
+      sprintf("averted_%s", c("median", "lower", "upper"))
+    )
+  )
+  expect_true(
+    all(averted_default$averted_median > averted_default$averted_lower &
+      averted_default$averted_median < averted_default$averted_upper)
+  )
+
+  # expectations for aggregation over full population
+  averted_aggregate <- outcomes_averted(
+    baseline = baseline, scenarios = scenarios, by_group = FALSE
+  )
+  expect_s3_class(averted_aggregate, "data.table")
+  expect_named(
+    averted_aggregate,
+    c("scenario", sprintf("averted_%s", c("median", "lower", "upper")))
+  )
+  expect_true(
+    all(averted_default$averted_median > averted_default$averted_lower &
+      averted_default$averted_median < averted_default$averted_upper)
+  )
+
+  # expectations when aggregation is not applied
+  averted_raw <- outcomes_averted(
+    baseline = baseline, scenarios = scenarios, summarise = FALSE
+  )
+  expect_s3_class(averted_aggregate, "data.table")
+  expect_named(
+    averted_raw,
+    c(
+      "scenario", "demography_group", "param_set",
+      "replicate", "outcomes_averted"
+    ),
+    ignore.order = TRUE
+  )
+
+  # expectations when aggregation is not applied
+  averted_agg_raw <- outcomes_averted(
+    baseline = baseline, scenarios = scenarios,
+    summarise = FALSE, by_group = FALSE
+  )
+  expect_s3_class(averted_agg_raw, "data.table")
+  expect_named(
+    averted_agg_raw,
+    c("scenario", "param_set", "replicate", "outcomes_averted"),
+    ignore.order = TRUE
+  )
+})