diff --git a/DESCRIPTION b/DESCRIPTION index 5cb5f2c..101610f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: daedalus Title: Model health, social, and economic costs of a pandemic using _DAEDALUS_ -Version: 0.0.14 +Version: 0.0.15 Authors@R: c( person("Pratik", "Gupte", , "p.gupte24@imperial.ac.uk", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-5294-7819")), @@ -16,7 +16,7 @@ License: MIT + file LICENSE URL: https://github.com/j-idea/daedalus BugReports: https://github.com/j-idea/daedalus/issues Depends: - R (>= 2.10) + R (>= 4.3) Imports: checkmate, cli, diff --git a/NEWS.md b/NEWS.md index a1c58c5..a61b4b3 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +# daedalus 0.0.15 + +This patch fixes an issue where closures were not ended in model runs where the closure began after the epidemic had stopped growing. This mostly affected edge cases of countries with large spare hospital capacity, and relatively late `response_time`s. In such cases the closure end time was assigned as the simulation end time, inflating costs related to closures. The fix: + +- Prevents closures from being activated by the `hospital_capacity` trigger if the epidemic is not growing, even if the hospital capacity threshold is crossed; + +- Prevents closures from being activated by the `response_time` trigger if the epidemic is not growing between stage 01 and stage 02. Closures are manually turned off if the epidemic is not growing ($R_t < 1.0$). + +Tests for different response times check that the model behaves as expected. + +## Miscellaneous changes + +- The package now requires a minimum R version >= v4.3 due to the use of the new function `array2DF()`; + +- `prepare_parameters.daedalus_country()` now provides a contact matrix scaled by the leading eigenvalue but not the demography for use in `r_eff()`. + # daedalus 0.0.14 This patch adds vaccine investment scenarios to _daedalus_. All model runs must now include an assumption about the level of advance vaccine investment, defaulting to "none". diff --git a/R/class_country.R b/R/class_country.R index 2787e4c..5b53fb1 100644 --- a/R/class_country.R +++ b/R/class_country.R @@ -418,7 +418,8 @@ prepare_parameters.daedalus_country <- function(x, ...) { # scale contacts by largest real eigenvalue cm <- get_data(x, "contact_matrix") eigv <- max(Re(eigen(cm)$values)) - cm <- (cm / eigv) %*% diag(1 / demography) + cm_unscaled <- cm / eigv + cm <- (cm_unscaled) %*% diag(1 / demography) cmw <- get_data(x, "contacts_workplace") cmw <- cmw / max(cmw) # max(cmw) is leading eigenvalue of diag matrix cmw @@ -431,6 +432,7 @@ prepare_parameters.daedalus_country <- function(x, ...) { list( demography = demography, contact_matrix = cm, + cm_unscaled = cm_unscaled, # for use in Rt calculations contacts_workplace = cmw, contacts_consumer_worker = cmcw, contacts_between_sectors = get_data(x, "contacts_between_sectors") # 0s diff --git a/R/closure.R b/R/closure.R index 81fd06f..5736fb7 100644 --- a/R/closure.R +++ b/R/closure.R @@ -17,15 +17,32 @@ make_response_threshold_event <- function(response_threshold) { } event_function <- function(time, state, parameters) { - # prevent flipping switch when checkEventFunc runs - # turn closure and excess hospitalisation switch on + # NOTE: prevent flipping switch when checkEventFunc runs + # NOTE: prevent response activation if epdedemic is not growing + state <- array( + state, + c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA) + ) + + rt <- r_eff(parameters[["r0"]], state, parameters[["cm_unscaled"]]) + + # NOTE: to ensure only first hosp threshold crossing is logged is_hosp_switch_on <- rlang::env_get(parameters[["mutables"]], "hosp_switch") + if (time != parameters[["min_time"]] && !is_hosp_switch_on) { - rlang::env_bind( - parameters[["mutables"]], - switch = TRUE, hosp_switch = TRUE, - closure_time_start = time + rlang::env_poke( + parameters[["mutables"]], "hosp_switch", TRUE ) + + # NOTE: trigger response and log closure start time only if + # epidemic is growing + if (rt >= 1.0) { + rlang::env_bind( + parameters[["mutables"]], + switch = TRUE, + closure_time_start = time + ) + } } as.numeric(state) } @@ -46,11 +63,8 @@ make_rt_end_event <- function() { c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA) ) - # because contacts are divided by demography during parameter prep - cm <- parameters[["contact_matrix"]] %*% diag(parameters[["demography"]]) - # arbitrary precision, may not be hit! - r_eff(parameters[["r0"]], state, cm) - 0.99 + r_eff(parameters[["r0"]], state, parameters[["cm_unscaled"]]) - 0.99 } event_function <- function(time, state, parameters) { diff --git a/R/daedalus.R b/R/daedalus.R index 0abfc44..b27e701 100644 --- a/R/daedalus.R +++ b/R/daedalus.R @@ -275,6 +275,28 @@ daedalus <- function(country, ] initial_state <- as.numeric(initial_state) + ### cancel closures if epidemic is not growing ### + state_temp <- array( + initial_state, + c(N_AGE_GROUPS, N_MODEL_COMPARTMENTS, N_ECON_STRATA, N_VACCINE_STRATA) + ) + # close intervention IFF epidemic is not growing + # NOTE: this step placed here as a conditional on r_eff < 1.0 is not + # practical in a root-finding function (Error: 'root too near initial point') + is_epidemic_growing <- r_eff( + parameters[["r0"]], + state_temp, + parameters[["cm_unscaled"]] + ) >= 1.0 + + if (!is_epidemic_growing) { + rlang::env_bind( + parameters[["mutables"]], + switch = FALSE, + closure_time_end = response_time + ) + } + # reset min time parameters[["min_time"]] <- response_time diff --git a/R/pkg_generics.R b/R/pkg_generics.R index 228e13b..f8ff92a 100644 --- a/R/pkg_generics.R +++ b/R/pkg_generics.R @@ -102,6 +102,9 @@ set_data <- function(x, ...) { #' - `contact_matrix`: the scaled contact matrix suitable for multiplication #' with \eqn{R_0} in force of infection calculations; #' +#' - `cm_unscaled`: the contact matrix multiplied column-wise by the demography +#' vector, suitable for calculations of \eqn{R_{\text{eff}}}; see [r_eff()]. +#' #' - `contacts_workplace`: the contacts in workplaces scaled by their largest #' value (which is the leading eigenvalue of the diagonal matrix of contacts); #' diff --git a/man/prepare_parameters.Rd b/man/prepare_parameters.Rd index 596b81d..c8e4892 100644 --- a/man/prepare_parameters.Rd +++ b/man/prepare_parameters.Rd @@ -45,6 +45,8 @@ The returned parameter list consists of: \item \code{demography}: the demography vector; \item \code{contact_matrix}: the scaled contact matrix suitable for multiplication with \eqn{R_0} in force of infection calculations; +\item \code{cm_unscaled}: the contact matrix multiplied column-wise by the demography +vector, suitable for calculations of \eqn{R_{\text{eff}}}; see \code{\link[=r_eff]{r_eff()}}. \item \code{contacts_workplace}: the contacts in workplaces scaled by their largest value (which is the leading eigenvalue of the diagonal matrix of contacts); \item \code{contacts_consumer_worker}: contacts in workplaces distributed in diff --git a/tests/testthat/test-closures.R b/tests/testthat/test-closures.R index ee55f6a..e35cfb6 100644 --- a/tests/testthat/test-closures.R +++ b/tests/testthat/test-closures.R @@ -201,4 +201,43 @@ test_that("Closures: correct logging of time limits", { response_data[["closure_info"]][["closure_time_end"]], time_end ) + + # Check that closures are terminated at the start time when + # the epidemic is not growing + dvx <- daedalus_vaccination("low", vax_start_time = 90) + tend <- 100 + response_time <- seq(10, 80, 10) + + # run scenario with very few susceptibles + invisible( + lapply( + response_time, function(x) { + data <- daedalus( + "United States", "influenza_1957", + initial_state_manual = list( + p_infectious = 0.9 + ), + response_strategy = "elimination", + response_time = x, + vaccine_investment = dvx, + time_end = tend + ) + + closure_info <- data$response_data$closure_info + + expect_identical( + closure_info$closure_time_end, + closure_info$closure_time_start + ) + expect_identical( + closure_info$closure_time_end, + x + ) + expect_identical( + closure_info$closure_duration, + 0.0 + ) + } + ) + ) })