Transplant #152 to new branch, WIP #127

Co-authored-by: Banky Ahadzie <[email protected]>

NAMESPACE

@@ -8,6 +8,7 @@ S3method(print,population)
 S3method(print,vaccination)
 export(as.intervention)
 export(as.vaccination)
+export(epidemic_peak)
 export(epidemic_size)
 export(intervention)
 export(is_contacts_intervention)

R/helpers.R

@@ -359,3 +359,78 @@ new_infections <- function(data,
   # return data
   data
 }
+
+#' Get the time and size of a compartment's highest peak
+#'
+#' Get the time and size of a compartment's highest peak for all
+#' demography groups.
+#'
+#' @details
+#' This is used for epidemics with a single peak. It is useful from
+#' a public health policy point of view to determine how bad an epidemic will
+#' be and when that happens.
+#'
+#' @param data A `<data.frame>` or `<data.table>` of model output, typically
+#' the output of a compartmental model.
+#' @param compartments A character vector for the compartments of interest.
+#' @return A `<data.table>` with columns "demography_group", "compartment",
+#' "time" and "value"; these specify the name of the demography group,
+#' the epidemiological compartment, and the peak time and value for each
+#' compartment in `compartments`.
+#'
+#' @export
+#'
+#' @examples
+#' # create a population
+#' uk_population <- population(
+#'   name = "UK population",
+#'   contact_matrix = matrix(1),
+#'   demography_vector = 67e6,
+#'   initial_conditions = matrix(
+#'     c(0.9999, 0.0001, 0, 0, 0),
+#'     nrow = 1, ncol = 5L
+#'   )
+#' )
+#'
+#' # run epidemic simulation with no vaccination or intervention
+#' data <- model_default(
+#'   population = uk_population,
+#'   time_end = 600
+#' )
+#'
+#' # get the timing and peak of the exposed and infectious compartment
+#' epidemic_peak(data, c("exposed", "infectious"))
+epidemic_peak <- function(data, compartments = "infectious") {
+  # solves "no visible binding for global variable"
+  compartment <- NULL
+  value <- NULL
+  time <- NULL
+
+  # basic input checks
+  checkmate::assert_data_frame(
+    data,
+    min.cols = 4L, min.rows = 1, any.missing = FALSE
+  )
+  checkmate::assert_names(
+    colnames(data),
+    must.include = c("time", "demography_group", "compartment", "value")
+  )
+  checkmate::assert_character(
+    compartments,
+    min.len = 1,
+    any.missing = FALSE, unique = TRUE,
+    null.ok = FALSE
+  )
+
+  data.table::setDT(data)
+  out <- data[compartment %in% compartments,
+    list(
+      time = data.table::first(time[value == max(value)]),
+      value = data.table::first(max(value))
+    ),
+    by = c("demography_group", "compartment")
+  ]
+
+  # return a data.table
+  out
+}

_pkgdown.yml

@@ -56,6 +56,7 @@ reference:
     desc: Functions to help with handling package classes and outputs.
     contents:
     - epidemic_size
+    - epidemic_peak
     - new_infections
   - title: Scenario comparison
     desc: Functions to help with comparing intervention scenarios.

tests/testthat/test-epidemic_peak.R

@@ -0,0 +1,178 @@
+#' Run the default model using UK population data from polymod and
+#' demography groups (0, 20, 40)
+test_data_default_model <- function() {
+  polymod <- socialmixr::polymod
+  contact_data <- socialmixr::contact_matrix(
+    polymod,
+    countries = "United Kingdom",
+    age.limits = c(0, 20, 40),
+    symmetric = TRUE
+  )
+
+  # prepare contact matrix
+  contact_matrix <- t(contact_data[["matrix"]])
+
+  # prepare the demography vector
+  demography_vector <- contact_data[["demography"]][["population"]]
+  names(demography_vector) <- rownames(contact_matrix)
+
+
+  # initial conditions: one in every 1 million is infected
+  initial_i <- 1e-6
+  initial_conditions <- c(
+    S = 1 - initial_i, E = 0, I = initial_i, R = 0, V = 0
+  )
+
+  # build for all age groups
+  initial_conditions <- rbind(
+    initial_conditions,
+    initial_conditions,
+    initial_conditions
+  )
+  rownames(initial_conditions) <- rownames(contact_matrix)
+
+  # prepare the population to model as affected by the epidemic
+  uk_population <- population(
+    name = "UK",
+    contact_matrix = contact_matrix,
+    demography_vector = demography_vector,
+    initial_conditions = initial_conditions
+  )
+
+  # run an epidemic model using `epidemic()`
+  output <- model_default(
+    population = uk_population,
+    transmission_rate = 1.5 / 7.0,
+    infectiousness_rate = 1.0 / 3.0,
+    recovery_rate = 1.0 / 7.0,
+    # intervention = list(contacts = close_schools),
+    time_end = 600, increment = 1.0
+  )
+
+  output
+}
+
+#' Run the vacamole model with no demography groups and a double vaccination regime
+test_data_vacamole <- function() {
+  # create a population, note eleven columns for compartments
+  population <- population(
+    contact_matrix = matrix(1),
+    demography_vector = 67e6,
+    initial_conditions = matrix(
+      c(0.9999, 0, 0, 0, 0, 0.0001, 0, 0, 0, 0, 0),
+      nrow = 1, ncol = 11L
+    )
+  )
+
+  # create a vaccination regime
+  double_vax <- vaccination(
+    nu = matrix(1e-3, ncol = 2, nrow = 1),
+    time_begin = matrix(c(10, 30), nrow = 1),
+    time_end = matrix(c(50, 80), nrow = 1)
+  )
+
+  # data <-
+  model_vacamole(
+    population = population,
+    vaccination = double_vax
+  )
+}
+
+#' Run the ebola model with no demography groups and a double vaccination regime
+test_data_ebola <- function() {
+  # create a population with 6 compartments
+  population <- population(
+    contact_matrix = matrix(1),
+    demography_vector = 14e6,
+    initial_conditions = matrix(
+      c(0.999998, 0.000001, 0.000001, 0, 0, 0),
+      nrow = 1, ncol = 6L
+    )
+  )
+
+  model_ebola(
+    population = population
+  )
+}
+
+test_that("Get the epidemic peak and size for the default model with 3 demography groups", {
+  # run epidemic model, expect no condition
+  data <- test_data_default_model()
+  expect_no_condition(
+    epidemic_peak(data)
+  )
+
+  out <- epidemic_peak(data)
+
+
+  # check for output type and contents
+  expect_s3_class(out, "data.frame")
+
+
+  expect_length(data, 4L)
+  expect_named(
+    out, c("compartment", "demography_group", "value", "time"),
+    ignore.order = TRUE
+  )
+
+  expect_identical(unique(out$compartment), "infectious")
+
+  # check for all positive values within the range 0 and total population size
+  expect_true(
+    all(
+      out[["value"]] <= data[data[["compartment"]] == "susceptible" & data[["time"]] == 0, value]
+    )
+  )
+})
+
+test_that("Get the epidemic peak and size for the vacamole model", {
+  # run epidemic model, expect no condition
+  data <- test_data_vacamole()
+  expect_no_condition(
+    epidemic_peak(data)
+  )
+  out <- epidemic_peak(data)
+
+  # check for output type and contents
+  expect_s3_class(out, "data.frame")
+  expect_length(data, 4L)
+  expect_named(
+    out, c("compartment", "demography_group", "value", "time"),
+    ignore.order = TRUE
+  )
+
+  expect_identical(unique(out$compartment), "infectious")
+
+  # check for all positive values within the range 0 and total population size
+  expect_true(
+    all(
+      out[["value"]] <= data[data[["compartment"]] == "susceptible" & data[["time"]] == 0, value]
+    )
+  )
+})
+
+test_that("Get the epidemic peak and size for the ebola model", {
+  # run epidemic model, expect no condition
+  data <- test_data_ebola()
+  expect_no_condition(
+    epidemic_peak(data)
+  )
+  out <- epidemic_peak(data)
+
+  # check for output type and contents
+  expect_s3_class(out, "data.frame")
+  expect_length(data, 5L)
+  expect_named(
+    out, c("compartment", "demography_group", "value", "time"),
+    ignore.order = TRUE
+  )
+
+  expect_identical(unique(out$compartment), "infectious")
+
+  # check for all positive values within the range 0 and total population size
+  expect_true(
+    all(
+      out[["value"]] <= data[data[["compartment"]] == "susceptible" & data[["time"]] == 0, value]
+    )
+  )
+})