From e7b5c4d8db1ed586810de84045526844dd1af625 Mon Sep 17 00:00:00 2001 From: Giovanni Charles Date: Fri, 15 Sep 2023 12:07:26 +0100 Subject: [PATCH] Fix bug with pev min_wait: * Update pev-epi min_wait test to directly check sample population * Update pe-epi min_wait test to check that pev_timestep is being updated on the first dose * Copy pev-epi min_wait test to mass pev * Update pev and epi processes to calculate targets based on time of first dose * Update infection immunity function to calculate pev immunity only after 3rd dose --- NEWS.md | 2 + R/human_infection.R | 5 +- R/mortality_processes.R | 3 +- R/pev.R | 23 ++++- R/variables.R | 20 +++-- tests/testthat/test-infection-integration.R | 4 +- tests/testthat/test-pev-epi.R | 33 ++++---- tests/testthat/test-pev.R | 94 +++++++++++++++++++-- 8 files changed, 143 insertions(+), 41 deletions(-) diff --git a/NEWS.md b/NEWS.md index 538311c6..860567da 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,8 @@ * Fix bug in competing hazards between mass and EPI vaccines. Where individuals can be enrolled onto both strategies if applied on the same timestep. + * Fix bug with min_wait. Min wait was working off of the final primary dose. It + now works of of the first dose. # malariasimulation 1.6.0 diff --git a/R/human_infection.R b/R/human_infection.R index b623a2d8..1cb521c4 100644 --- a/R/human_infection.R +++ b/R/human_infection.R @@ -129,9 +129,10 @@ calculate_infections <- function( # calculate vaccine efficacy vaccine_efficacy <- rep(0, length(source_vector)) - vaccine_times <- variables$pev_timestep$get_values(source_vector) - vaccinated <- vaccine_times > -1 + vaccine_times <- variables$last_eff_pev_timestep$get_values(source_vector) pev_profile <- variables$pev_profile$get_values(source_vector) + # get vector of individuals who have received their 3rd dose + vaccinated <- vaccine_times > -1 pev_profile <- pev_profile[vaccinated] if (length(vaccinated) > 0) { antibodies <- calculate_pev_antibodies( diff --git a/R/mortality_processes.R b/R/mortality_processes.R index 65b090cb..9d95a707 100644 --- a/R/mortality_processes.R +++ b/R/mortality_processes.R @@ -106,7 +106,8 @@ reset_target <- function(variables, events, target, state, timestep) { variables$drug_time$queue_update(-1, target) # vaccination - variables$pev_timestep$queue_update(-1, target) + variables$last_pev_timestep$queue_update(-1, target) + variables$last_eff_pev_timestep$queue_update(-1, target) variables$pev_profile$queue_update(-1, target) variables$tbv_vaccinated$queue_update(-1, target) diff --git a/R/pev.R b/R/pev.R index 28eb2819..039b6044 100644 --- a/R/pev.R +++ b/R/pev.R @@ -40,7 +40,7 @@ create_epi_pev_process <- function( if (parameters$pev_epi_min_wait == 0) { target <- to_vaccinate$to_vector() } else { - not_recently_vaccinated <- variables$pev_timestep$get_index_of( + not_recently_vaccinated <- variables$last_pev_timestep$get_index_of( a = max(timestep - parameters$pev_epi_min_wait, 0), b = timestep )$not(TRUE) @@ -56,6 +56,9 @@ create_epi_pev_process <- function( ) ] + # Update the latest vaccination time + variables$last_pev_timestep$queue_update(timestep, target) + schedule_vaccination( target, events, @@ -91,7 +94,7 @@ create_mass_pev_listener <- function( if (parameters$mass_pev_min_wait == 0) { target <- in_age_group } else { - not_recently_vaccinated <- variables$pev_timestep$get_index_of( + not_recently_vaccinated <- variables$last_pev_timestep$get_index_of( a = max(timestep - parameters$mass_pev_min_wait, 0), b = timestep )$not(TRUE) @@ -116,6 +119,11 @@ create_mass_pev_listener <- function( correlations ) ] + + # Update the latest vaccination time + variables$last_pev_timestep$queue_update(timestep, target) + + # Schedule future doses schedule_vaccination( target, events, @@ -162,7 +170,7 @@ schedule_vaccination <- function( create_pev_efficacy_listener <- function(variables, pev_profile_index) { function(timestep, target) { if (target$size() > 0) { - variables$pev_timestep$queue_update(timestep, target) + variables$last_eff_pev_timestep$queue_update(timestep, target) variables$pev_profile$queue_update(pev_profile_index, target) } } @@ -185,7 +193,8 @@ create_pev_booster_listener <- function( force(coverage) function(timestep, target) { target <- sample_bitset(target, coverage) - variables$pev_timestep$queue_update(timestep, target) + 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) renderer$render(render_name, target$size(), timestep) @@ -240,6 +249,12 @@ attach_pev_dose_listeners <- function( dose_events[[d]]$add_listener( create_dosage_renderer(renderer, strategy, d) ) + # update last vaccination on every primary dose + dose_events[[d]]$add_listener( + function(t, target) { + variables$last_pev_timestep$queue_update(t, target) + } + ) if (d == length(dose_events)) { dose_events[[d]]$add_listener( create_pev_efficacy_listener( diff --git a/R/variables.R b/R/variables.R index 0931b3f7..b84074e4 100644 --- a/R/variables.R +++ b/R/variables.R @@ -1,5 +1,4 @@ -#' @title Define model variables -#' @description +#' @title Define model variables #' @description #' create_variables creates the human and mosquito variables for #' the model. Variables are used to track real data for each individual over #' time, they are read and updated by processes @@ -18,10 +17,13 @@ #' * ID - Acquired immunity to detectability #' * zeta - Heterogeneity of human individuals #' * zeta_group - Discretised heterogeneity of human individuals -#' * pev_timestep - The timestep of the last pev vaccination (-1 if there -#' haven't been any) -#' * pev_profile - The index of the profile of the last administered pev vaccine -#' (-1 if there haven't been any) +#' * last_pev_timestep - The timestep of the last pev vaccination (-1 if there +#' * last_eff_pev_timestep - The timestep of the last efficacious pev +#' vaccination, including final primary dose and booster doses (-1 if there have not been any) +#' * pev_profile - The index of the efficacy profile of any pev vaccinations. +#' Not set until the final primary dose. +#' This is only set on the final primary dose and subsequent booster doses +#' (-1 otherwise) #' * tbv_vaccinated - The timstep of the last tbv vaccination (-1 if there #' haven't been any #' * net_time - The timestep when a net was last put up (-1 if never) @@ -188,7 +190,8 @@ create_variables <- function(parameters) { drug <- individual::IntegerVariable$new(rep(0, size)) drug_time <- individual::IntegerVariable$new(rep(-1, size)) - pev_timestep <- individual::IntegerVariable$new(rep(-1, size)) + last_pev_timestep <- individual::IntegerVariable$new(rep(-1, size)) + last_eff_pev_timestep <- individual::IntegerVariable$new(rep(-1, size)) pev_profile <- individual::IntegerVariable$new(rep(-1, size)) tbv_vaccinated <- individual::DoubleVariable$new(rep(-1, size)) @@ -215,7 +218,8 @@ create_variables <- function(parameters) { infectivity = infectivity, drug = drug, drug_time = drug_time, - pev_timestep = pev_timestep, + last_pev_timestep = last_pev_timestep, + last_eff_pev_timestep = last_eff_pev_timestep, pev_profile = pev_profile, tbv_vaccinated = tbv_vaccinated, net_time = net_time, diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R index a8e6033e..238a2ae6 100644 --- a/tests/testthat/test-infection-integration.R +++ b/tests/testthat/test-infection-integration.R @@ -133,7 +133,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati ), drug = individual::DoubleVariable$new(c(1, 2, 0, 0)), drug_time = individual::DoubleVariable$new(c(20, 30, -1, -1)), - pev_timestep = individual::DoubleVariable$new(c(-1, 10, 40, -1)), + last_eff_pev_timestep = individual::DoubleVariable$new(c(-1, 10, 40, -1)), pev_profile = individual::IntegerVariable$new(c(-1, 1, 2, -1)), ib = individual::DoubleVariable$new(c(.2, .3, .5, .9)) ) @@ -366,7 +366,7 @@ test_that('prophylaxis is considered for medicated humans', { ), drug = individual::DoubleVariable$new(c(0, 2, 1, 0)), drug_time = individual::DoubleVariable$new(c(-1, 49, 40, -1)), - pev_timestep = individual::DoubleVariable$new(c(-1, -1, -1, -1)), + last_eff_pev_timestep = individual::DoubleVariable$new(c(-1, -1, -1, -1)), pev_profile = individual::IntegerVariable$new(c(-1, -1, -1, -1)), ib = individual::DoubleVariable$new(c(.2, .3, .5, .9)) ) diff --git a/tests/testthat/test-pev-epi.R b/tests/testthat/test-pev-epi.R index 3dac524d..ef9148a3 100644 --- a/tests/testthat/test-pev-epi.R +++ b/tests/testthat/test-pev-epi.R @@ -86,49 +86,48 @@ test_that('pev epi targets correct age and respects min_wait', { variables$birth <- individual::IntegerVariable$new( -c(18, 18, 2.9, 18, 18) * 365 + timestep ) - variables$pev_timestep <- mock_integer( + variables$last_pev_timestep <- mock_integer( c(50, -1, -1, 4*365, -1) ) + variables$pev_profile <- mock_integer( + c(1, -1, -1, 1, -1) + ) - events$pev_epi_doses <- lapply(events$pev_epi_doses, mock_event) - + correlations <- get_correlation_parameters(parameters) process <- create_epi_pev_process( variables, events, parameters, - get_correlation_parameters(parameters), + correlations, parameters$pev_epi_coverages, parameters$pev_epi_timesteps ) + sample_mock <- mockery::mock(c(TRUE, TRUE, FALSE)) mockery::stub( process, 'sample_intervention', - mockery::mock(c(TRUE, TRUE, FALSE)) + sample_mock ) process(timestep) mockery::expect_args( - events$pev_epi_doses[[1]]$schedule, + sample_mock, 1, - c(1, 2), - parameters$pev_doses[[1]] + c(1, 2, 5), + 'pev', + .8, + correlations ) mockery::expect_args( - events$pev_epi_doses[[2]]$schedule, + variables$last_pev_timestep$queue_update_mock(), 1, - c(1, 2), - parameters$pev_doses[[2]] + timestep, + c(1, 2) ) - mockery::expect_args( - events$pev_epi_doses[[3]]$schedule, - 1, - c(1, 2), - parameters$pev_doses[[3]] - ) }) test_that('EPI ignores individuals scheduled for mass vaccination', { diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R index c80dc457..30007cd3 100644 --- a/tests/testthat/test-pev.R +++ b/tests/testthat/test-pev.R @@ -97,7 +97,7 @@ test_that('Infection considers pev efficacy', { variables$birth <- individual::IntegerVariable$new( -c(8, 2.9, 3.2, 18.4) * 365 - 100 ) - variables$pev_timestep <- individual::IntegerVariable$new( + variables$last_eff_pev_timestep <- individual::IntegerVariable$new( c(-1, -1, 50, 50 + 30) ) variables$pev_profile <- individual::IntegerVariable$new( @@ -299,7 +299,10 @@ test_that('Mass boosters update profile params and reschedule correctly', { variables$birth <- individual::IntegerVariable$new( -c(2.9, 3.2, 18.4) * 365 + 100 ) - variables$pev_timestep <- mock_double( + variables$last_pev_timestep <- mock_double( + c(50, 50, 50) + ) + variables$last_eff_pev_timestep <- mock_double( c(50, 50, 50) ) variables$pev_profile <- mock_integer( @@ -326,7 +329,13 @@ test_that('Mass boosters update profile params and reschedule correctly', { listener(timestep, individual::Bitset$new(3)$insert(c(1, 2, 3))) expect_bitset_update( - variables$pev_timestep$queue_update_mock(), + variables$last_pev_timestep$queue_update_mock(), + timestep, + c(1, 2, 3) + ) + + expect_bitset_update( + variables$last_eff_pev_timestep$queue_update_mock(), timestep, c(1, 2, 3) ) @@ -368,7 +377,10 @@ test_that('Mass booster coverages sample subpopulations correctly', { variables$birth <- individual::IntegerVariable$new( -c(2.9, 3.2, 18.4) * 365 + 100 ) - variables$pev_timestep <- mock_double( + variables$last_pev_timestep <- mock_double( + c(50, 50, 50) + ) + variables$last_eff_pev_timestep <- mock_double( c(50, 50, 50) ) variables$pev_profile <- mock_integer( @@ -395,7 +407,13 @@ test_that('Mass booster coverages sample subpopulations correctly', { mockery::expect_args(sample_mock, 1, target, .9) expect_bitset_update( - variables$pev_timestep$queue_update_mock(), + variables$last_pev_timestep$queue_update_mock(), + timestep, + c(2, 3) + ) + + expect_bitset_update( + variables$last_eff_pev_timestep$queue_update_mock(), timestep, c(2, 3) ) @@ -413,6 +431,68 @@ test_that('Mass booster coverages sample subpopulations correctly', { ) }) +test_that('mass pev targets correct age and respects min_wait', { + timestep <- 5*365 + parameters <- get_parameters(list(human_population = 5)) + parameters <- set_mass_pev( + parameters, + profile = rtss_profile, + timesteps = c(4, 5) * 365, + coverages = c(0.8, 0.8), + min_ages = 0, + max_ages = 19 * 365, + min_wait = 2*365, + booster_timestep = c(1, 6) * 30, + booster_coverage = c(.9, .8), + booster_profile = list(rtss_booster_profile, rtss_booster_profile) + ) + events <- create_events(parameters) + variables <- create_variables(parameters) + variables$birth <- individual::IntegerVariable$new( + -c(18, 18, 30, 18, 18) * 365 + timestep + ) + variables$last_pev_timestep <- mock_integer( + c(50, -1, -1, 4*365, -1) + ) + + variables$pev_profile <- mock_integer( + c(1, -1, -1, 1, -1) + ) + + correlations <- get_correlation_parameters(parameters) + listener <- create_mass_pev_listener( + variables, + events, + parameters, + get_correlation_parameters(parameters) + ) + + sample_mock <- mockery::mock(c(TRUE, TRUE, FALSE)) + mockery::stub( + listener, + 'sample_intervention', + sample_mock + ) + + listener(timestep) + + mockery::expect_args( + sample_mock, + 1, + c(1, 2, 5), + 'pev', + .8, + correlations + ) + + mockery::expect_args( + variables$last_pev_timestep$queue_update_mock(), + 1, + timestep, + c(1, 2) + ) +}) + test_that('Mass efficacy listener works correctly', { timestep <- 50 parameters <- get_parameters() @@ -430,7 +510,7 @@ test_that('Mass efficacy listener works correctly', { ) variables <- create_variables(parameters) - variables$pev_timestep <- mock_integer(c(-1, -1, -1)) + variables$last_eff_pev_timestep <- mock_integer(c(-1, -1, -1)) variables$pev_profile <- mock_integer(c(-1, -1, -1)) listener <- create_pev_efficacy_listener(variables, 1) @@ -438,7 +518,7 @@ test_that('Mass efficacy listener works correctly', { # vaccinated time expect_bitset_update( - variables$pev_timestep$queue_update_mock(), + variables$last_eff_pev_timestep$queue_update_mock(), timestep, c(1, 2, 3) )