Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add epidemic_peak() function #240

Merged
merged 1 commit into from
Jun 5, 2024
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -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)
75 changes: 75 additions & 0 deletions R/helpers.R
Original file line number Diff line number Diff line change
@@ -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
}
1 change: 1 addition & 0 deletions _pkgdown.yml
Original file line number Diff line number Diff line change
@@ -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.
50 changes: 50 additions & 0 deletions man/epidemic_peak.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

178 changes: 178 additions & 0 deletions tests/testthat/test-epidemic_peak.R
Original file line number Diff line number Diff line change
@@ -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]
)
)
})