Implement time varying coverage for PEV boosters: #270

giovannic · 2023-11-16T17:04:30Z

Some functionality to allow booster coverage to change over time. This "booster timed coverage" is multiplied with the previous "booster coverage" which is used to model dropouts.

Changes:

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

Fixes #269

pwinskill · 2023-11-17T10:53:00Z

R/pev_parameters.R

@@ -88,6 +90,8 @@ set_pev_epi <- function(
  min_wait,
  booster_timestep,
  booster_coverage,
+  booster_timed_coverage = NULL,


I wonder if this can be simplified given that boosters will always be associated with the primary serires.

Would the following reduced arguments work:
timesteps define the vaccine coverage change timepoints at which primary series coverage or booster coverage changes.
booster_coverage (like, coverages) is a vector of length timesteps of booster coverages
Boosters are given at booster_timestep post primary series. For clarity this could be renamed (booster_lags or booster_spacing?)

That would remove the requirement for the two addtional arguments here: booster_timed_coverage booster_timed_coverage_timestep

@giovannic @lmhaile
Any thoughts on if the above would be sensible/helpful?
Happy to re-review otherwise

Hi Pete, thanks for boosting! Yes I agree that booster_timed_coverage_timestep could be removed in favor of timesteps, as boosters will always be associated with the primary series.

Also agree that booster_timestep could be renamed to booster_spacing to avoid confusion.

We may need to keep booster_coverage and booster_timed_coverage as is to avoid breaking the code of people implementing flat booster coverage across time.

Otherwise the changes look great-- everything is working as expected now!

Hey! Sorry for the delay!

I think I get it and it makes a lot of sense. I can get on it later this week.

Just to confirm, booster_coverage would need to be a matrix? Because we have (potentially) several boosters and several timesteps.

Sorry for the delay while I was OOO Giovanni! Yes, booster_coverage would need to be a matrix to account for more than one booster.

Hope you had a nice new year and just let me know if you have other thoughts!

lmhaile · 2023-12-05T12:58:57Z

Hi Giovanni!

I'm having some trouble getting the booster functionality to work as intended-- reproducible example below. Based on the parameter changes I would expect boosters to disappear from the model output by the fourth year, but they seem to be consistent throughout the time series-- perhaps there is something I am missing?


remotes::install_github('mrc-ide/malariasimulation@feat/time_boosters')
library(malariasimulation)
plot_doses <- function(){
  cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
  
  output$month <- ceiling(output$timestep / month)
  doses <- output[, c(grep("n_pev" , names(output)), grep("month", names(output)))]
  doses <- aggregate(cbind(doses[1:4]),
                     by = list(doses$month), 
                     FUN = sum)
  doses <- as.matrix(t(doses[, -1]))
  
  barplot(doses, xlab = "Month",
          ylab = "Number of doses",
          col = cols[1:6], space = 0, axes = T,
          beside = FALSE, xaxs = "i", yaxs = "i",
          ylim = c(0, 230))
  grid(lty = 2, col = "grey80", lwd = 0.5);box()
  axis(side = 1, lty = 1, col = "black", pos = 0)
  legend("topleft", box.lty = 0, legend = c("Dose 1","Dose 2","Dose 3","Booster 1"),
         fill = cols[1:6], bg="transparent", cex = 0.8, y.intersp = 1.5)
}


year <- 365
month <- 30
sim_length <- 3 * year
human_population <- 10000
starting_EIR <- 20

simparams <- get_parameters(list(
  human_population = human_population,
  clinical_incidence_rendering_min_ages = 0,
  clinical_incidence_rendering_max_ages = 5 * year,
  individual_mosquitoes = FALSE
)
)
simparams <- set_equilibrium(parameters = simparams, init_EIR = starting_EIR)


# add changing coverage
rtssepiparams2 <- set_pev_epi(
  simparams,
  profile = rtss_profile, 
  timesteps = 1 * year, 
  coverages = 1, 
  age = 5 * month, 
  min_wait = 0, 
  booster_timestep = c(12 * month),  # booster is always implemented 1 year after third dose
  booster_coverage = c(1), # 
  booster_timed_coverage = c(1, 1, 0.0, 0.0),
  booster_timed_coverage_timestep = c(1, 2, 3, 4) * year,
  booster_profile = list(rtss_booster_profile) 
)

output <- run_simulation(timesteps = sim_length * 2, parameters = rtssepiparams2)
plot_doses()

Looking at the PR, it seems like the issue may have to do with this function, although I am not entirely sure:

malariasimulation/R/pev.R

Lines 179 to 217 in a87b4f1

    
           create_pev_booster_listener <- function( 
        
             variables, 
        
             coverage, 
        
             timed_coverage = NULL, 
        
             timed_coverage_timestep = NULL, 
        
             booster_number, 
        
             pev_profile_index, 
        
             next_booster_event, 
        
             next_booster_delay, 
        
             renderer, 
        
             strategy 
        
             ) { 
        
             render_name <- paste0("n_pev_", strategy, "_booster_", booster_number) 
        
             renderer$set_default(render_name, 0) 
        
             force(next_booster_event) # because R lazy evaluation is rubbish 
        
             force(next_booster_delay) 
        
             force(coverage) 
        
             function(timestep, target) { 
        
               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) 
        
               renderer$render(render_name, target$size(), timestep) 
        
               if (!is.null(next_booster_event)) { 
        
                 next_booster_event$schedule(target, next_booster_delay) 
        
               } 
        
             } 
        
           }

Just let me know if I can help troubleshoot!

giovannic · 2023-12-06T16:39:31Z

Ah! That was my mistake, sorry. set_pev_epi got the name of the parameter wrong, it's "timed_booster_coverage" not "booster_timed_coverage".

This should be fixed with the latest commit.

* 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

* rename booster_timesteps -> booster_spacing * make booster_coverage a matrix (timestep x booster doses) * update validation * update implementation of attach_booster_listener * fix regression tests * update documentation * update vignettes * fix correlation examples

lmhaile

Reran my same quick test with the parameter name changes and all looks well- much more intuitive to me now!


remotes::install_github('mrc-ide/malariasimulation@feat/time_boosters')
library(malariasimulation)
plot_doses <- function(){
  cols <- c("#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
  
  output$month <- ceiling(output$timestep / month)
  doses <- output[, c(grep("n_pev" , names(output)), grep("month", names(output)))]
  doses <- aggregate(cbind(doses[1:4]),
                     by = list(doses$month), 
                     FUN = sum)
  doses <- as.matrix(t(doses[, -1]))
  
  barplot(doses, xlab = "Month",
          ylab = "Number of doses",
          col = cols[1:6], space = 0, axes = T,
          beside = FALSE, xaxs = "i", yaxs = "i",
          ylim = c(0, 230))
  grid(lty = 2, col = "grey80", lwd = 0.5);box()
  axis(side = 1, lty = 1, col = "black", pos = 0)
  legend("topleft", box.lty = 0, legend = c("Dose 1","Dose 2","Dose 3","Booster 1"),
         fill = cols[1:6], bg="transparent", cex = 0.8, y.intersp = 1.5)
}
year <- 365
month <- 30
sim_length <- 3 * year
human_population <- 10000
starting_EIR <- 20

simparams <- get_parameters(list(
  human_population = human_population,
  clinical_incidence_rendering_min_ages = 0,
  clinical_incidence_rendering_max_ages = 5 * year,
  individual_mosquitoes = FALSE
)
)
simparams <- set_equilibrium(parameters = simparams, init_EIR = starting_EIR)


# add changing coverage
rtssepiparams2 <- set_pev_epi(
  simparams,
  profile = rtss_profile, 
  timesteps = c(1, 2, 3, 4) * year,
  coverages = c(1, 1, 1, 1), 
  age = 5 * month, 
  min_wait = 0, 
  booster_spacing= c(12 * month),  # booster is always implemented 1 year after third dose
  booster_coverage = matrix(c(1, 1, 0.0, 0.0), ncol = 1),
  booster_profile = list(rtss_booster_profile) 
)

output <- run_simulation(timesteps = sim_length * 2, parameters = rtssepiparams2)
plot_doses()

We may need to forecast this change to some of the other vaccine modellers (Kelly, Nora) as the booster_spacing and booster_coverage changes may break their code (a Teams post?). Thanks for the update, Giovanni!

pwinskill

This looks great - just one minor comment

pwinskill · 2024-01-05T15:59:03Z

R/pev_parameters.R

@@ -86,12 +85,13 @@ set_pev_epi <- function(
  timesteps,
  age,
  min_wait,
-  booster_timestep,
+  booster_spacing,


Should we have a check that the vector of booster spacings are monotonically increasing when more than one booster is used - setting otherwise gives quite a cryptic Error: delay must be >= 0 error

giovannic requested a review from pwinskill November 16, 2023 17:04

pwinskill reviewed Nov 17, 2023

View reviewed changes

giovannic force-pushed the feat/time_boosters branch from ddfb268 to 685dea1 Compare January 4, 2024 15:51

giovannic added 2 commits January 4, 2024 15:52

Implement time varying coverage for PEV boosters:

73cbfb4

* 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

Fix PEV EPI parameterisation bug

d49c8dd

giovannic force-pushed the feat/time_boosters branch 3 times, most recently from d13e23e to bdef3f5 Compare January 4, 2024 17:28

giovannic force-pushed the feat/time_boosters branch from bdef3f5 to 91247b4 Compare January 4, 2024 17:51

giovannic requested review from pwinskill and lmhaile January 4, 2024 17:51

lmhaile approved these changes Jan 5, 2024

View reviewed changes

pwinskill reviewed Jan 5, 2024

View reviewed changes

Check that booster_spacing is monotonically increasing

45486ad

giovannic requested a review from pwinskill January 9, 2024 10:11

pwinskill approved these changes Jan 16, 2024

View reviewed changes

giovannic merged commit b0ad081 into dev Jan 16, 2024
4 checks passed

giovannic deleted the feat/time_boosters branch January 16, 2024 10:32

giovannic mentioned this pull request Sep 11, 2024

Dev #335

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement time varying coverage for PEV boosters: #270

Implement time varying coverage for PEV boosters: #270

giovannic commented Nov 16, 2023 •

edited

Loading

pwinskill Nov 17, 2023 •

edited

Loading

pwinskill Dec 13, 2023

lmhaile Dec 13, 2023

giovannic Dec 13, 2023

giovannic Dec 21, 2023

lmhaile Jan 2, 2024

lmhaile commented Dec 5, 2023

giovannic commented Dec 6, 2023

lmhaile left a comment

pwinskill left a comment

pwinskill Jan 5, 2024

Implement time varying coverage for PEV boosters: #270

Implement time varying coverage for PEV boosters: #270

Conversation

giovannic commented Nov 16, 2023 • edited Loading

pwinskill Nov 17, 2023 • edited Loading

Choose a reason for hiding this comment

pwinskill Dec 13, 2023

Choose a reason for hiding this comment

lmhaile Dec 13, 2023

Choose a reason for hiding this comment

giovannic Dec 13, 2023

Choose a reason for hiding this comment

giovannic Dec 21, 2023

Choose a reason for hiding this comment

lmhaile Jan 2, 2024

Choose a reason for hiding this comment

lmhaile commented Dec 5, 2023

giovannic commented Dec 6, 2023

lmhaile left a comment

Choose a reason for hiding this comment

pwinskill left a comment

Choose a reason for hiding this comment

pwinskill Jan 5, 2024

Choose a reason for hiding this comment

giovannic commented Nov 16, 2023 •

edited

Loading

pwinskill Nov 17, 2023 •

edited

Loading