From a87b4f1b8510126160e6fc054090508d8f6d8438 Mon Sep 17 00:00:00 2001 From: Giovanni Charles Date: Thu, 16 Nov 2023 16:50:30 +0000 Subject: [PATCH] Implement time varying coverage for PEV boosters: * Update documentation * Update parameterisation and tests * Update booster listener to incorporate timed coverage * Update attach_events to pass timed coverage parameters * Test timed coverage sampling * Bonus: fix a biology test --- R/events.R | 4 + R/pev.R | 18 ++++- R/pev_parameters.R | 19 +++++ man/CorrelationParameters.Rd | 36 ++++----- man/set_mass_pev.Rd | 6 ++ man/set_pev_epi.Rd | 6 ++ tests/testthat/test-biology.R | 4 +- tests/testthat/test-pev-epi.R | 41 ++++++++++ tests/testthat/test-pev.R | 148 ++++++++++++++++++++++++++++++++++ 9 files changed, 261 insertions(+), 21 deletions(-) diff --git a/R/events.R b/R/events.R index 879d72bd..21207e5d 100644 --- a/R/events.R +++ b/R/events.R @@ -133,6 +133,8 @@ attach_event_listeners <- function( events$mass_pev_boosters, parameters$mass_pev_booster_timestep, parameters$mass_pev_booster_coverage, + parameters$mass_pev_timed_booster_coverage, + parameters$mass_pev_timed_booster_coverage_timestep, parameters$mass_pev_profile_indices, 'mass', renderer @@ -147,6 +149,8 @@ attach_event_listeners <- function( events$pev_epi_boosters, parameters$pev_epi_booster_timestep, parameters$pev_epi_booster_coverage, + parameters$pev_epi_timed_booster_coverage, + parameters$pev_epi_timed_booster_coverage_timestep, parameters$pev_epi_profile_indices, 'epi', renderer diff --git a/R/pev.R b/R/pev.R index 039b6044..d37ff6d3 100644 --- a/R/pev.R +++ b/R/pev.R @@ -179,6 +179,8 @@ create_pev_efficacy_listener <- function(variables, pev_profile_index) { create_pev_booster_listener <- function( variables, coverage, + timed_coverage = NULL, + timed_coverage_timestep = NULL, booster_number, pev_profile_index, next_booster_event, @@ -192,7 +194,17 @@ create_pev_booster_listener <- function( force(next_booster_delay) force(coverage) function(timestep, target) { - target <- sample_bitset(target, coverage) + if (is.null(timed_coverage)) { + t_coverage <- 1 + } else { + t_coverage <- timed_coverage[ + match_timestep(timed_coverage_timestep, timestep) + ] + } + target <- sample_bitset( + target, + coverage * t_coverage + ) variables$last_pev_timestep$queue_update(timestep, target) variables$last_eff_pev_timestep$queue_update(timestep, target) variables$pev_profile$queue_update(pev_profile_index, target) @@ -240,6 +252,8 @@ attach_pev_dose_listeners <- function( booster_events, booster_delays, booster_coverages, + booster_timed_coverage, + booster_timed_coverage_timestep, pev_profile_indices, strategy, renderer @@ -303,6 +317,8 @@ attach_pev_dose_listeners <- function( create_pev_booster_listener( variables = variables, coverage = booster_coverages[[b]], + timed_coverage = booster_timed_coverage, + timed_coverage_timestep = booster_timed_coverage_timestep, booster_number = b, pev_profile_index = pev_profile_indices[[b + 1]], next_booster_event = next_booster_event, diff --git a/R/pev_parameters.R b/R/pev_parameters.R index aeb6bd5b..c74f3a94 100644 --- a/R/pev_parameters.R +++ b/R/pev_parameters.R @@ -74,6 +74,8 @@ rtss_booster_profile <- create_pev_profile( #' an individual being vaccinated under one scheme and vaccinated under another. #' @param booster_timestep the timesteps (following the final dose) at which booster vaccinations are administered #' @param booster_coverage the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series) +#' @param booster_timed_coverage a time varying proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series), set in time with `booster_coverage_timestep` +#' @param booster_timed_coverage_timestep a vector of timesteps to change the time varying coverage specified in `booster_timed_coverage` #' @param booster_profile list of booster vaccine profiles, of type #' PEVProfile, for each timestep in booster_timeteps #' @param seasonal_boosters logical, if TRUE the first booster timestep is @@ -88,6 +90,8 @@ set_pev_epi <- function( min_wait, booster_timestep, booster_coverage, + booster_timed_coverage = NULL, + booster_timed_coverage_timestep = NULL, booster_profile, seasonal_boosters = FALSE ) { @@ -114,6 +118,9 @@ set_pev_epi <- function( if (!all(c(length(booster_coverage), length(booster_timestep), length(booster_profile)) == length(booster_timestep))) { stop('booster_timestep and booster_coverage and booster_profile does not align') } + if (length(booster_timed_coverage) != length(booster_timed_coverage_timestep)) { + stop("booster_coverage_timestep must be the same length as booster_coverage") + } # Index the new vaccine profiles profile_list <- c(list(profile), booster_profile) @@ -128,6 +135,9 @@ set_pev_epi <- function( parameters$pev_epi_booster_timestep <- booster_timestep parameters$pev_epi_min_wait <- min_wait parameters$pev_epi_booster_coverage <- booster_coverage + parameters$pev_epi_booster_timed_coverage <- booster_timed_coverage + parameters$pev_epi_booster_timed_coverage_timestep <- booster_timed_coverage_timestep + parameters$pev_epi_booster_coverage <- booster_coverage parameters$pev_epi_profile_indices <- profile_indices parameters$pev_epi_seasonal_boosters <- seasonal_boosters parameters @@ -149,6 +159,8 @@ set_pev_epi <- function( #' @param max_ages for the target population, inclusive (in timesteps) #' @param booster_timestep the timesteps (following the initial vaccination) at which booster vaccinations are administered #' @param booster_coverage the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series) +#' @param booster_timed_coverage a time varying proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series), set in time with `booster_coverage_timestep` +#' @param booster_timed_coverage_timestep a vector of timesteps to change the time varying coverage specified in `booster_timed_coverage` #' @param booster_profile list of booster vaccine profiles, of type #' PEVProfile, for each timestep in booster_timeteps #' @export @@ -162,6 +174,8 @@ set_mass_pev <- function( min_wait, booster_timestep, booster_coverage, + booster_timed_coverage = NULL, + booster_timed_coverage_timestep = NULL, booster_profile ) { stopifnot(all(timesteps >= 1)) @@ -176,6 +190,9 @@ set_mass_pev <- function( if (!all(c(length(booster_coverage), length(booster_timestep), length(booster_profile)) == length(booster_timestep))) { stop('booster_timestep, booster_coverage and booster_profile does not align') } + if (length(booster_timed_coverage) != length(booster_timed_coverage_timestep)) { + stop("booster_coverage_timestep must be the same length as booster_coverage") + } # Index the new vaccine profiles profile_list <- c(list(profile), booster_profile) @@ -191,6 +208,8 @@ set_mass_pev <- function( parameters$mass_pev_min_wait <- min_wait parameters$mass_pev_booster_timestep <- booster_timestep parameters$mass_pev_booster_coverage <- booster_coverage + parameters$mass_pev_booster_timed_coverage <- booster_timed_coverage + parameters$mass_pev_booster_timed_coverage_timestep <- booster_timed_coverage_timestep parameters$mass_pev_profile_indices <- profile_indices parameters } diff --git a/man/CorrelationParameters.Rd b/man/CorrelationParameters.Rd index c2e6ada7..b1d22578 100644 --- a/man/CorrelationParameters.Rd +++ b/man/CorrelationParameters.Rd @@ -14,17 +14,17 @@ Describes an event in the simulation \section{Methods}{ \subsection{Public methods}{ \itemize{ -\item \href{#method-CorrelationParameters-new}{\code{CorrelationParameters$new()}} -\item \href{#method-CorrelationParameters-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}} -\item \href{#method-CorrelationParameters-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}} -\item \href{#method-CorrelationParameters-sigma}{\code{CorrelationParameters$sigma()}} -\item \href{#method-CorrelationParameters-mvnorm}{\code{CorrelationParameters$mvnorm()}} -\item \href{#method-CorrelationParameters-clone}{\code{CorrelationParameters$clone()}} +\item \href{#method-new}{\code{CorrelationParameters$new()}} +\item \href{#method-inter_round_rho}{\code{CorrelationParameters$inter_round_rho()}} +\item \href{#method-inter_intervention_rho}{\code{CorrelationParameters$inter_intervention_rho()}} +\item \href{#method-sigma}{\code{CorrelationParameters$sigma()}} +\item \href{#method-mvnorm}{\code{CorrelationParameters$mvnorm()}} +\item \href{#method-clone}{\code{CorrelationParameters$clone()}} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-CorrelationParameters-new}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-new}{}}} \subsection{Method \code{new()}}{ initialise correlation parameters \subsection{Usage}{ @@ -40,8 +40,8 @@ initialise correlation parameters } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_round_rho}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-inter_round_rho}{}}} \subsection{Method \code{inter_round_rho()}}{ Add rho between rounds \subsection{Usage}{ @@ -60,8 +60,8 @@ the intervention} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-CorrelationParameters-inter_intervention_rho}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-inter_intervention_rho}{}}} \subsection{Method \code{inter_intervention_rho()}}{ Add rho between interventions \subsection{Usage}{ @@ -83,8 +83,8 @@ the intervention} } } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-CorrelationParameters-sigma}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-sigma}{}}} \subsection{Method \code{sigma()}}{ Standard deviation of each intervention between rounds \subsection{Usage}{ @@ -93,8 +93,8 @@ Standard deviation of each intervention between rounds } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-CorrelationParameters-mvnorm}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-mvnorm}{}}} \subsection{Method \code{mvnorm()}}{ multivariate norm draws for these parameters \subsection{Usage}{ @@ -103,8 +103,8 @@ multivariate norm draws for these parameters } \if{html}{\out{
}} -\if{html}{\out{}} -\if{latex}{\out{\hypertarget{method-CorrelationParameters-clone}{}}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-clone}{}}} \subsection{Method \code{clone()}}{ The objects of this class are cloneable with this method. \subsection{Usage}{ diff --git a/man/set_mass_pev.Rd b/man/set_mass_pev.Rd index 7c28547f..e088e0bb 100644 --- a/man/set_mass_pev.Rd +++ b/man/set_mass_pev.Rd @@ -14,6 +14,8 @@ set_mass_pev( min_wait, booster_timestep, booster_coverage, + booster_timed_coverage = NULL, + booster_timed_coverage_timestep = NULL, booster_profile ) } @@ -38,6 +40,10 @@ time between an individual being vaccinated under one scheme and vaccinated unde \item{booster_coverage}{the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series)} +\item{booster_timed_coverage}{a time varying proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series), set in time with \code{booster_coverage_timestep}} + +\item{booster_timed_coverage_timestep}{a vector of timesteps to change the time varying coverage specified in \code{booster_timed_coverage}} + \item{booster_profile}{list of booster vaccine profiles, of type PEVProfile, for each timestep in booster_timeteps} } diff --git a/man/set_pev_epi.Rd b/man/set_pev_epi.Rd index 4dea14ed..2f7e7ced 100644 --- a/man/set_pev_epi.Rd +++ b/man/set_pev_epi.Rd @@ -13,6 +13,8 @@ set_pev_epi( min_wait, booster_timestep, booster_coverage, + booster_timed_coverage = NULL, + booster_timed_coverage_timestep = NULL, booster_profile, seasonal_boosters = FALSE ) @@ -39,6 +41,10 @@ an individual being vaccinated under one scheme and vaccinated under another.} \item{booster_coverage}{the proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series)} +\item{booster_timed_coverage}{a time varying proportion of the vaccinated population relative to the last vaccination (whether a previous booster or the primary series), set in time with \code{booster_coverage_timestep}} + +\item{booster_timed_coverage_timestep}{a vector of timesteps to change the time varying coverage specified in \code{booster_timed_coverage}} + \item{booster_profile}{list of booster vaccine profiles, of type PEVProfile, for each timestep in booster_timeteps} diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R index 5895d45b..c739f7e7 100644 --- a/tests/testthat/test-biology.R +++ b/tests/testthat/test-biology.R @@ -92,13 +92,13 @@ test_that('FOIM is consistent with equilibrium', { psi <- unique_biting_rate(age, parameters) zeta <- variables$zeta$get_values() .pi <- human_pi(psi, zeta) - calculate_foim(a, sum(.pi * variables$infectivity$get_values())) + calculate_foim(a, sum(.pi * variables$infectivity$get_values()), 1.) } ) expect_equal( expected_foim, actual_foim, - tolerance = 1e-4 + tolerance = 1e-3 ) }) diff --git a/tests/testthat/test-pev-epi.R b/tests/testthat/test-pev-epi.R index ef9148a3..26541b9c 100644 --- a/tests/testthat/test-pev-epi.R +++ b/tests/testthat/test-pev-epi.R @@ -51,6 +51,47 @@ test_that('pev epi strategy parameterisation works', { ) }) +test_that('I can add time varying booster coverage to the pev epi strategy', { + parameters <- get_parameters() + parameters <- set_pev_epi( + parameters, + profile = rtss_profile, + coverages = c(0.1, 0.8), + timesteps = c(10, 100), + min_wait = 0, + age = 5 * 30, + booster_timestep = c(18, 36) * 30, + booster_coverage = c(.9, .8), + booster_timed_coverage = c(.5, .7), + booster_timed_coverage_timestep = c(365, 2*365), + booster_profile = list(rtss_booster_profile, rtss_booster_profile) + ) + expect_equal(parameters$pev_epi_booster_timestep, c(18, 36) * 30) + expect_equal(parameters$pev_epi_booster_coverage, c(.9, .8)) + expect_equal(parameters$pev_epi_booster_timed_coverage, c(.5, .7)) + expect_equal(parameters$pev_epi_booster_timed_coverage_timestep, c(365, 2*365)) + expect_equal(parameters$pev_profiles, list(rtss_profile, rtss_booster_profile, rtss_booster_profile)) + + expect_error( + parameters <- set_pev_epi( + parameters, + profile = rtss_profile, + coverages = c(0.1, 0.8), + timesteps = c(10, 100), + min_wait = 0, + age = 5 * 30, + booster_timestep = c(18, 36) * 30, + booster_coverage = c(.9, .8), + booster_timed_coverage = c(.5, .7), + booster_timed_coverage_timestep = 365, + booster_profile = list(rtss_booster_profile, rtss_booster_profile) + ), + "booster_coverage_timestep must be the same length as booster_coverage", + fixed = TRUE + ) +}) + + test_that('pev epi fails pre-emptively with unaligned booster parameters', { parameters <- get_parameters() expect_error( diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R index 30007cd3..16f65849 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -54,6 +54,49 @@ test_that('Mass vaccination strategy parameterisation works', { ) }) +test_that('I can add time varying booster coverage to the mass pev strategy', { + parameters <- get_parameters() + parameters <- set_mass_pev( + parameters, + profile = rtss_profile, + timesteps = 10, + coverages = 0.8, + min_wait = 0, + min_ages = 5 * 30, + max_ages = 17 * 30, + booster_timestep = c(18, 36) * 30, + booster_coverage = c(.9, .8), + booster_timed_coverage = c(.5, .7), + booster_timed_coverage_timestep = c(365, 2*365), + booster_profile = list(rtss_booster_profile, rtss_booster_profile) + ) + + expect_equal(parameters$mass_pev_booster_timestep, c(18, 36) * 30) + expect_equal(parameters$mass_pev_booster_timed_coverage, c(.5, .7)) + expect_equal(parameters$mass_pev_booster_timed_coverage_timestep, c(365, 2*365)) + expect_equal(parameters$pev_profiles, list(rtss_profile, rtss_booster_profile, rtss_booster_profile)) + + expect_error( + parameters <- set_mass_pev( + parameters, + profile = rtss_profile, + timesteps = 10, + coverages = 0.8, + min_wait = 0, + min_ages = 5 * 30, + max_ages = 17 * 30, + booster_timestep = c(18, 36) * 30, + booster_coverage = c(.9, .8), + booster_timed_coverage = c(.5, .7), + booster_timed_coverage_timestep = 365, + booster_profile = list(rtss_booster_profile, rtss_booster_profile) + ), + "booster_coverage_timestep must be the same length as booster_coverage", + fixed = TRUE + ) +}) + + test_that('Mass vaccination fails pre-emptively for unaligned booster parameters', { parameters <- get_parameters() expect_error( @@ -557,6 +600,8 @@ test_that('Mass dose events are not ruined by lazy evaluation', { booster_events = events$mass_pev_boosters, booster_delays = parameters$mass_pev_booster_timestep, booster_coverages = parameters$mass_pev_booster_coverage, + booster_timed_coverage = NULL, + booster_timed_coverage_timestep = NULL, pev_profile_indices = parameters$mass_pev_profile_indices, strategy = 'mass', renderer = mock_render(1) @@ -612,4 +657,107 @@ test_that('Efficacies are calculated correctly', { ) }) +test_that('pev timed booster coverage works with NULL', { + timestep <- 5 * 365 + parameters <- get_parameters(list(human_population = 5)) + parameters <- set_pev_epi( + parameters, + profile = rtss_profile, + timesteps = 10, + coverages = 0.8, + min_wait = 6 * 30, + age = 18 * 365, + booster_timestep = c(3, 12) * 30, + booster_coverage = c(.9, .8), + booster_timed_coverage = c(.5, .7), + booster_timed_coverage_timestep = c(timestep, timestep + 365), + booster_profile = list(rtss_booster_profile, rtss_booster_profile), + ) + events <- create_events(parameters) + + booster_event <- mock_event(events$pev_epi_boosters[[1]]) + + listener <- create_pev_booster_listener( + variables = create_variables(parameters), + coverage = .9, + timed_coverage = NULL, + timed_coverage_timestep = NULL, + booster_number = 1, + pev_profile_index = 2, + next_booster_event = booster_event, + next_booster_delay = 9 * 30, + renderer = mock_render(timestep), + strategy = 'epi' + ) + + target <- individual::Bitset$new(5)$insert(seq(5)) + + mock_sample_bitset = mockery::mock(individual::Bitset$new(5)$insert(c(1, 2))) + mockery::stub( + listener, + 'sample_bitset', + mock_sample_bitset + ) + + listener(timestep, target) + + mockery::expect_args( + mock_sample_bitset, + 1, + target, + .9 + ) +}) + + +test_that('pev boosters take into account the timed coverage', { + timestep <- 5 * 365 + parameters <- get_parameters(list(human_population = 5)) + parameters <- set_pev_epi( + parameters, + profile = rtss_profile, + timesteps = 10, + coverages = 0.8, + min_wait = 6 * 30, + age = 18 * 365, + booster_timestep = c(3, 12) * 30, + booster_coverage = c(.9, .8), + booster_timed_coverage = c(.5, .7), + booster_timed_coverage_timestep = c(timestep, timestep + 365), + booster_profile = list(rtss_booster_profile, rtss_booster_profile), + ) + events <- create_events(parameters) + + booster_event <- mock_event(events$pev_epi_boosters[[1]]) + + listener <- create_pev_booster_listener( + variables = create_variables(parameters), + coverage = .9, + timed_coverage = c(.5, .7), + timed_coverage_timestep = c(timestep, timestep + 365), + booster_number = 1, + pev_profile_index = 2, + next_booster_event = booster_event, + next_booster_delay = 9 * 30, + renderer = mock_render(timestep), + strategy = 'epi' + ) + + target <- individual::Bitset$new(5)$insert(seq(5)) + + mock_sample_bitset = mockery::mock(individual::Bitset$new(5)$insert(c(1, 2))) + mockery::stub( + listener, + 'sample_bitset', + mock_sample_bitset + ) + + listener(timestep, target) + mockery::expect_args( + mock_sample_bitset, + 1, + target, + .45 + ) +})