diff --git a/NAMESPACE b/NAMESPACE index 320867c0..1fcc5134 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/helpers.R b/R/helpers.R index 9669e7be..c7b344ff 100644 --- a/R/helpers.R +++ b/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 `` or `` of model output, typically +#' the output of a compartmental model. +#' @param compartments A character vector for the compartments of interest. +#' @return A `` 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 +} diff --git a/_pkgdown.yml b/_pkgdown.yml index 97d86a09..4f3bf7be 100644 --- a/_pkgdown.yml +++ b/_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. diff --git a/man/epidemic_peak.Rd b/man/epidemic_peak.Rd new file mode 100644 index 00000000..1b41f98c --- /dev/null +++ b/man/epidemic_peak.Rd @@ -0,0 +1,50 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/helpers.R +\name{epidemic_peak} +\alias{epidemic_peak} +\title{Get the time and size of a compartment's highest peak} +\usage{ +epidemic_peak(data, compartments = "infectious") +} +\arguments{ +\item{data}{A \verb{} or \verb{} of model output, typically +the output of a compartmental model.} + +\item{compartments}{A character vector for the compartments of interest.} +} +\value{ +A \verb{} 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 \code{compartments}. +} +\description{ +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. +} +\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")) +} diff --git a/tests/testthat/test-epidemic_peak.R b/tests/testthat/test-epidemic_peak.R new file mode 100644 index 00000000..07b3b06a --- /dev/null +++ b/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] + ) + ) +})