-
Notifications
You must be signed in to change notification settings - Fork 14
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
Competing hazards recovery infection resolution #300
Competing hazards recovery infection resolution #300
Conversation
6bf22cc
to
903547a
Compare
…(containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions.
91d50f4
to
b36be83
Compare
… for use in infection outcome.
…(containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions.
… for use in infection outcome.
…olution (any rendering in the infection outcome does not occur), I've split the incidence renderer into three functions: age rendering, n incidence rendering and p incidence rendering. This allows much more flexibility in when these are called in the code sequence. All age inputs are now consolidated and rendered in a single function. In addition, we intitate and populate these incidence columns with 0 when they are used.
Corrected dt variables recovery rate assignment when antimalarial resistance is modelled.
Merge branch 'competing_hazards_recovery_infection_resolution' of https://github.com/mrc-ide/malariasimulation into competing_hazards_recovery_infection_resolution # Conflicts: # R/disease_progression.R # R/human_infection.R
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a massive change. Thank you so much for the effort, it has been really thorough.
Few changes for readability and performance
R/disease_progression.R
Outdated
if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ | ||
dt_v <- dt_input | ||
} else if (isTRUE(parameters$antimalarial_resistance)) { | ||
dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We probably want to cover all cases and use boolean instead of bitwise operators.
if(isFALSE(parameters$antimalarial_resistance) & length(dt_input) == 1 & is.numeric(dt_input)){ | |
dt_v <- dt_input | |
} else if (isTRUE(parameters$antimalarial_resistance)) { | |
dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) | |
} | |
if(isFALSE(parameters$antimalarial_resistance)){ | |
stopifnot(length(dt_input) == 1 && is.numeric(dt_input)) # please move this check to pre-simulation time for efficiency | |
dt_v <- dt_input | |
} else { | |
dt_v <- dt_input$get_values(variables$state$get_index_of("Tr")) | |
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I've made this change.
R/disease_progression.R
Outdated
recovery_rates <- numeric(length = parameters$human_population) | ||
recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd | ||
recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da | ||
recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du | ||
recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v | ||
recovery_outcome$set_rates(recovery_rates) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
It's perhaps more efficient to create and update a new variable called recovery_rates
in the progression processes. Then all of this becomes.
recovery_rates <- numeric(length = parameters$human_population) | |
recovery_rates[variables$state$get_index_of("D")$to_vector()] <- 1/parameters$dd | |
recovery_rates[variables$state$get_index_of("A")$to_vector()] <- 1/parameters$da | |
recovery_rates[variables$state$get_index_of("U")$to_vector()] <- 1/parameters$du | |
recovery_rates[variables$state$get_index_of("Tr")$to_vector()] <- 1/dt_v | |
recovery_outcome$set_rates(recovery_rates) | |
recovery_outcome$set_rates(variables$recovery_rates$get_values()) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think this actually tidies up some of the issues I was having integrating with @tbreweric 's dt variable. I'm guessing that if we're storing all the recovery rates as a variable then we can replace the dt variable with a recovery_rates variable and still store all the relevant information. Let me know what you think with what I've done.
R/competing_hazards.R
Outdated
stash_source_humans = function(source_humans){ | ||
self$source_humans <- source_humans | ||
}, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
We don't need this anymore? It's ugly 😆
stash_source_humans = function(source_humans){ | |
self$source_humans <- source_humans | |
}, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yep - gone!
R/human_infection.R
Outdated
#' Updates human states and variables to represent asymptomatic/clinical/severe | ||
#' and treated malaria; and resulting boosts in immunity | ||
#' This function ends with the assignment of rates of infection to the competing | ||
#' hazard resolution object. Boosts immunity given infectious bites. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
#' hazard resolution object. Boosts immunity given infectious bites. | |
#' hazard resolution object and boosts immunity given infectious bites. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Made this change.
R/human_infection.R
Outdated
source_humans <- variables$state$get_index_of( | ||
c('S', 'A', 'U'))$and(bitten_humans) | ||
infection_outcome$stash_source_humans(source_humans) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
infection_outcome$stash_source_humans(source_humans) |
R/processes.R
Outdated
|
||
infection_outcome <- CompetingOutcome$new( | ||
targeted_process = function(timestep, target){ | ||
infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
infection_process_resolved_hazard(timestep, target, infection_outcome$source_humans, | |
infection_outcome_process(timestep, target, |
R/processes.R
Outdated
|
||
recovery_outcome <- CompetingOutcome$new( | ||
targeted_process = function(timestep, target){ | ||
recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer) | |
recovery_outcome_process(timestep, target, variables, parameters, renderer) |
R/processes.R
Outdated
create_infection_rates_process <- function(timestep){ | ||
biting_process( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
couldn't we keep the function factory in biting_process.R?
In this file we've been calling function factories from *_process.R files. But now we're making the functions in processes.R. It feels a bit inconsistent. And verbose since we have a lot of processes to create.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree - I found this a little confusing to follow
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think I understand what you're saying, and I've had a go at implementing this.
R/processes.R
Outdated
# Disease Progression | ||
# =================== | ||
|
||
create_recovery_rates_process <- function(timestep){ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Again, making the function here when it could be in disease_progression.R
R/human_infection.R
Outdated
infection_process_resolved_hazard <- function( | ||
timestep, | ||
infected_humans, | ||
source_humans, |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is not used
source_humans, |
Also if you could rebase with dev, then we can run benchmarks to get an idea of the performance tradeoff. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks @RJSheppard. Some additional comments from me.
I was looking at the plots you submitted with the initial PR. Are these still valid now? The increase in prevalence for mid-EIRs was quite marked (if I was interpreting the plots correctly), which I think we should check and be cautious about.
I think, if those changes are that apparent we should probably look at PfPr~EIR plot similar to that in Penny et al. If possible lines for current and after changes would be helpful.
R/human_infection.R
Outdated
) | ||
|
||
## capture infection rates to resolve in competing hazards | ||
infection_rates <- numeric(length = parameters$human_population) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Would the following make it more clear that rates that are not overwritten in the following line are 0?
infection_rates <- numeric(length = parameters$human_population) | |
infection_rates <- rep(0, length = parameters$human_population) |
R/processes.R
Outdated
create_infection_rates_process <- function(timestep){ | ||
biting_process( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I agree - I found this a little confusing to follow
#' @param timestep current target | ||
#' | ||
#' @noRd | ||
incidence_probability_renderer <- function( |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Maybe not for this PR - but does anyone use these probability of incidence outputs?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Whenever I look at this it always takes me a few to work out what it even means. I'd be happy to remove it, but agree it would be good to check if anyone uses it.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
These are a smoother version of the incidence outputs, these should be closer to what Griffin's model outputs.
It seems unclear to me which, if any, should be removed
R/render.R
Outdated
unique_age_combinations <- as.data.frame(unique(rbind(cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), | ||
cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), | ||
cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), | ||
cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), | ||
cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)))) | ||
|
||
ordered_unique_age_combinations <- unique_age_combinations[ | ||
with(unique_age_combinations, order(V1, V2)), | ||
] |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Suggestion for clarity (untested)
unique_age_combinations <- as.data.frame(unique(rbind(cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), | |
cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), | |
cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), | |
cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), | |
cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages)))) | |
ordered_unique_age_combinations <- unique_age_combinations[ | |
with(unique_age_combinations, order(V1, V2)), | |
] | |
age_ranges <- rbind( | |
cbind(parameters$prevalence_rendering_min_ages, parameters$prevalence_rendering_max_ages), | |
cbind(parameters$incidence_rendering_min_ages, parameters$incidence_rendering_max_ages), | |
cbind(parameters$clinical_incidence_rendering_min_ages, parameters$clinical_incidence_rendering_max_ages), | |
cbind(parameters$severe_incidence_rendering_min_ages, parameters$severe_incidence_rendering_max_ages), | |
cbind(parameters$age_group_rendering_min_ages, parameters$age_group_rendering_max_ages) | |
) | |
unique_age_combinations <- as.data.frame(unique(age_ranges)) | |
ordered_unique_age_combinations <- unique_age_combinations[order(unique_age_combinations$V1, unique_age_combinations$V2), ] | |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Added this now - all seems to be working.
for (i in seq_along(parameters$age_group_rendering_min_ages)) { | ||
lower <- parameters$age_group_rendering_min_ages[[i]] | ||
upper <- parameters$age_group_rendering_max_ages[[i]] | ||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@giovannic As these unique combinations of age ranges will not change throughout the simulation can we calculate them once when create_age_group_renderer() is called and not everytime when age_group_renderer() is called?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes! I thought I commented this, but that's gone missing. Good catch
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'm not sure I fully understand this. Is it as simple as moving the not repeating parts outside of the time looping section? The actual age rendering needs to occur at each time step, but the age rendering groups only need to be assigned once. I think I've done this now, but please let me know if it isn't what you meant!
R/utils.R
Outdated
#'@param rate rate | ||
#'@noRd | ||
rate_to_prob <- function(rate){ | ||
1- exp(-rate) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
1- exp(-rate) | |
1 - exp(-rate) |
@pwinskill I'll rebase and make some plots as requested. |
…now occurs in biting_process.R. Similarly, create_recovery_rates_process functionality is now in disease_progression.R stash_source_humans in competing hazard resolution is now redundant and has been removed from infection processes. recovery rates are now assigned to a variable which is updated as infections, recoveries (via update_infection) and deaths occur. This also takes into account recovery time of treated due to slow parasite clearance and replaced the dt variable. I've also noticed (I think) that Tr->S and U->S recoveries can be done in a single step, so have attempted to combine these at risk of making the code more complicated! Added a metapopulation incidence rendering wrapper function for initialising multiple renderers. Rewrote create_age_group_renderer for clarity following feedback. Some tests have been adjusted to reflect these changes.
…(containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions.
… for use in infection outcome.
…olution (any rendering in the infection outcome does not occur), I've split the incidence renderer into three functions: age rendering, n incidence rendering and p incidence rendering. This allows much more flexibility in when these are called in the code sequence. All age inputs are now consolidated and rendered in a single function. In addition, we intitate and populate these incidence columns with 0 when they are used. Rebasing
Corrected dt variables recovery rate assignment when antimalarial resistance is modelled.
…(containing an outcome following resolution and rates that that outcome occurs for each individual in the population) and CompetingHazard (contains all CompetingOutcomes and process function to resolve hazards). processes.R: create_biting_process now replaced with create_infection_rates_process and biting_process, that uses essentially the same functional pathway, but which stops when the rates of infection are generated, storing them in the infection_outcome object within human_infection.R Similarly, disease_progression_process(es) are replaced with new functions: create_recovery_rates_process and calculate_recovery_rates, that get the recovery rates for each individual based on disease state. We then add these processes to the process list, including the resolution function create_infection_recovery_hazard_process$resolve. This process (create_infection_recovery_hazard_process) contains infection_outcome and recovery_outcome. infection_outcome calls infection_process_resolved_hazard, which follows the rest of the infection pathway functions: assigning new infections to different human states. Note that this requires probability of infection which gets used in the incidence renderer function. recovery_outcome also updates human states and infectivities. I've adjusted some of the tests to reflect the changes in functions required to generate results. Removed white space deletions.
… for use in infection outcome.
…now occurs in biting_process.R. Similarly, create_recovery_rates_process functionality is now in disease_progression.R stash_source_humans in competing hazard resolution is now redundant and has been removed from infection processes. recovery rates are now assigned to a variable which is updated as infections, recoveries (via update_infection) and deaths occur. This also takes into account recovery time of treated due to slow parasite clearance and replaced the dt variable. I've also noticed (I think) that Tr->S and U->S recoveries can be done in a single step, so have attempted to combine these at risk of making the code more complicated! Added a metapopulation incidence rendering wrapper function for initialising multiple renderers. Rewrote create_age_group_renderer for clarity following feedback. Some tests have been adjusted to reflect these changes.
This is how benchmark results would change (along with a 95% confidence interval in relative change) if 09e6e8c is merged into dev:
Further explanation regarding interpretation and methodology can be found in the documentation. |
We only need to perform the selection for individuals that passed the occurs check, not the entire population.
/benchmark |
/benchmark |
This is how benchmark results would change (along with a 95% confidence interval in relative change) if 1f7310b is merged into dev:
Further explanation regarding interpretation and methodology can be found in the documentation. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I might try to get those missing 40% back later but it doesn't seem worth blocking merging this anymore.
7a118d1
to
4c98559
Compare
/benchmark |
/benchmark |
…very_infection_resolution
This is how benchmark results would change (along with a 95% confidence interval in relative change) if 0bd7c10 is merged into dev:
Further explanation regarding interpretation and methodology can be found in the documentation. |
Eventually this should go into individual.
/benchmark |
This is how benchmark results would change (along with a 95% confidence interval in relative change) if ded6876 is merged into dev:
Further explanation regarding interpretation and methodology can be found in the documentation. |
Any news on this? Would like to merge this |
All done from my end. There's not much more to be gained I think. |
This PR resolves the competing hazards of infection and recovery. Giovanni wrote code to generate two competing hazards objects: CompetingOutcome (in which a function describing what happens to individuals following competing hazard resolution and the rates that that outcome occurs) and CompetingHazard (which links to each of the CompetingOutcome objects that require resolution, and which contains the process function that resolves the competing hazards), found in competing_hazards.R
We first replace create_biting_process and create_progression_process(es) with create_infection_rates_process and create_recovery_rate_process.
Create_infection_rate_process uses the renamed “biting_rate” (which was previously: create_biting_rate_process) which continues as previously, passing through the infection outcome object, and now ends when the rates of infection are generated. Note that we are currently only resolving infections in SA and U compartments, NOT D or Tr (consistent with deterministic model). Note also that new infections are not rendered until after occurrence of resolution.
Similarly, all recovery/progression processes are now replaced with a single create_recovery_rates_process function, that uses calculate_recovery_rates, based on individual state rates of recovery.
The create_infection_recovery_hazard_process then uses the CompetingHazard object (to which each competing hazard outcome is assigned) to resolve competing hazards and determine which outcome occurs to which individuals.
The infection_outcome uses a new function: infection_process_resolved_hazard, which picks up from where we previously left off in the infection process, rendering incidence and assigning infections severities. Note that this function requires the infection rates to allow the incidence renderer to function.
The recovery_outcome then uses recovery_process_resolved_hazard, streamlining the progression functions into the relevant update functions. This function also passes in the dt variable for antimalarial resistance @tbreweric. (Hopefully it will still work!).
I ran simulations to compare the current default model (with mortality as the last process) with the resolved hazards. We do see something of an increase in incidence (which increases with EIR) and prevalence.
competing_hazards_mortality_process_position.pdf
competing_hazards_resolved_timesteps.pdf
Human states also see fewer susceptibles, and is possibly deviating from the equilibrium a little.
competing_hazards_resolved_states.pdf
competing_hazards_resolved_states_time.pdf
Acquired immunities seem to be increasing a bit through time - consistent with higher onward infection following higher overall prevalence, particularly at lower EIR (higher plots).
competing_hazards_resolved_immunity.pdf
I also checked through the expected infections and recoveries compared with the equilibrium solution, and found reasonable consistency, and resolution of the hazards. Let me know your thoughts!