Add hypnozoites to the P.v model #321

RJSheppard · 2024-07-03T14:28:58Z

We must model hynozoite batch formation via bite infections, decay of hypnozoite batches through, new infections via relapse and the reseting of batch number to 0 for new births.

We create a hypnozoite variable.

The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each.

The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR). biting_rate.R

In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection.

The number of bites is considered in calculating probability of a new infection.

To calculate the infection rates for vivax, we take the biting infection rates and add them to the relase rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later.

A new function to resolve between bites and relapse infections is created. For each bite we increase hte number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches.

We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level.

When an individual dies, we reset the hypnozoite batches back to 0.

We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups.

Added a test to check relapses occur as expected.

pwinskill

Hi Richard,
A few Qs to followup on for this one, alongside some suggestions 🌟

pwinskill · 2024-07-16T08:45:05Z

R/biting_process.R

+        bitten_humans$insert(bitten)
+        if(parameters$parasite == "vivax"){
+          # p.v must pass through the number of bites each individual receives
+          n_bites_each <- table(bitten)


Does using table here work for when someone doesn't get sampled in the fast_weighted_sample() call above?

For example:

t1 <- malariasimulation:::fast_weighted_sample(20, c(0.8, 0.001, 0.3)) table(t1)

Person index 2 doesn't show up. It wasn't clear to me when n_bites_each is used downstream in human_infection.R how this is dealt with?

I think tabulate(bitten, nbins = length(lambda)) is an efficient way to do this if the above is an issue

So I was getting around that by passing through the full table - with column names indicating the individuals that had been bitten. But I think using tabulate can make lots of things a bit simpler. I've had another go with this. Let me know what you think. I've tried to use bitsets where I can.

R/human_infection.R

R/parameters.R

R/processes.R

pwinskill

...a few more changes suggested

R/human_infection.R

pwinskill · 2024-07-17T14:47:01Z

R/biting_process.R

-          n_bites_each <- table(bitten)
+          # p.v must pass through the number of bites per person
+          bites_per_person <- tabulate(bitten, nbins = length(lambda))
+          bitten_humans <- individual::IntegerVariable$new(initial_values = bites_per_person)


I see how this seems tidier, but I have a few concerns.

We now have a single variable bitten_humans that means different things for Pf and Pv and ...

that leads to a mismatch in the calculate_infections() function, which takes the argument bitten_humans() currently defined as a bitset of bitten humans.

Would it be clearer to:

Create the bitset of bitten_humans based on n_bites and lambda for both Pf and Pf and pass this through to calculate_infections() as an argument.

Additionally create a IntegerVariable for vivax called and storing n_bites. This can be accessed in calculate_infections() directly in that function as you're already passing variables as an argument.

I've given this a go (and rebased) and I think it's ready for review.

RJSheppard · 2024-10-14T09:17:44Z

/benchmark

github-actions · 2024-10-14T10:23:12Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if e2adee2 is merged into dev:

✔️large_population: 2.3m -> 2.32m [-0.2%, +2.23%]
✔️small_population: 32.7s -> 32.4s [-3.82%, +1.96%]

Further explanation regarding interpretation and methodology can be found in the documentation.
Plots and raw data are available as artifacts of the workflow run.

giovannic

Functionally looks fine to me. I have suggested some cleaner ways to organise the code and the competing hazard stuff.

giovannic · 2024-10-14T09:38:23Z

R/human_infection.R

+    if(parameters$parasite == "falciparum"){
+      source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)
+    } else if (parameters$parasite == "vivax"){
+      ## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination
+      hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)
+      source_humans <- bitten_humans$copy()$or(hypnozoites_humans)
+    }


Suggested change

if(parameters$parasite == "falciparum"){

source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)

} else if (parameters$parasite == "vivax"){

## source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination

hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)

source_humans <- bitten_humans$copy()$or(hypnozoites_humans)

}

if(parameters$parasite == "falciparum"){

source_humans <- variables$state$get_index_of(c('S','A','U'))$and(bitten_humans)

} else if (parameters$parasite == "vivax"){

# source_humans must include individuals with hypnozoites which may be impacted by prophylaxis/vaccination

hypnozoites_humans <- variables$hypnozoites$get_index_of(set = 0)$not(T)

source_humans <- bitten_humans$copy()$or(hypnozoites_humans)

}

giovannic · 2024-10-14T09:48:52Z

R/human_infection.R

@@ -48,91 +51,124 @@ simulate_infection <- function(
 #' @description Infection rates are stored in the infection outcome competing hazards object
 #' @param variables a list of all of the model variables
 #' @param bitten_humans bitset of bitten humans
+#' @param n_bites_per_person vector of number of bites each person receives (p.v only)
 #' @param parameters model parameters
 #' @param renderer model render object
 #' @param timestep current timestep
 #' @noRd
 calculate_infections <- function(


[optional] This is bordering on being cleaner as 4 functions: calculate_falciparum_infections, calculate_vivax_infections and pev.R:pev_efficacy and treatment_efficacy.

I only say this because the if (parasite==...) conditions do make this function long and trickier to follow

giovannic · 2024-10-14T10:01:48Z

R/human_infection.R

@@ -145,14 +181,16 @@ calculate_infections <- function(
 #' @param renderer model render object
 #' @param parameters model parameters
 #' @param prob vector of population probabilities of infection
+#' @param relative_rates relative rates of hypnozoite relapse relative to total infection rate
 #' @noRd
 infection_outcome_process <- function(


Again, this function is quite long. It would be cleaner as falciparum_infection_outcome and vivax_infection_outcome.

giovannic · 2024-10-14T10:06:59Z

R/processes.R

+                                prob = rate_to_prob(infection_outcome$rates),
+                                relative_rates = infection_outcome$relative_rates)


I've only just realised how ugly and genius this is...

I don't yet know how I feel about it

I think the above edit to CompetingOutcome is a better way to deal with optional arguments for the processes. That way it generalises to different competing hazard scenarios and you don't have to continue adding fields to the R6 object for more arguments.

giovannic · 2024-10-14T10:16:31Z

R/competing_hazards.R

    set_rates = function(target, rates){
      stopifnot(target$size() == length(rates))

      self$target$copy_from(target)
      self$rates <- rates
    },


I think a more readable, generalisable - and perhaps less error prone - method would be to explicity accept extra arguments for your outcome when you set rates.

Suggested change

set_rates = function(target, rates){

stopifnot(target$size() == length(rates))

self$target$copy_from(target)

self$rates <- rates

},

set_rates = function(target, rates, ...){

stopifnot(target$size() == length(rates))

self$target$copy_from(target)

self$rates <- rates

self$args <- list(...)

},

giovannic · 2024-10-14T10:17:56Z

R/competing_hazards.R

+    set_relative_rates = function(target, relative_rates){
+      stopifnot(target$size() == length(relative_rates))
+
+      self$target$copy_from(target)
+      self$relative_rates <- relative_rates
+    },


Suggested change

set_relative_rates = function(target, relative_rates){

stopifnot(target$size() == length(relative_rates))

self$target$copy_from(target)

self$relative_rates <- relative_rates

},

giovannic · 2024-10-14T10:18:19Z

R/competing_hazards.R

    },
    target = NULL,
-    rates = NULL
+    rates = NULL,
+    relative_rates = NULL


Suggested change

relative_rates = NULL

args = NULL

giovannic · 2024-10-14T10:22:03Z

R/competing_hazards.R

+
+      self$target$copy_from(target)
+      self$relative_rates <- relative_rates
+    },
    execute = function(t, target){
      private$targeted_process(t, target)


Suggested change

private$targeted_process(t, target)

do.call(private$targeted_process, c(c(t, target), self$args)

giovannic · 2024-10-14T10:25:59Z

R/processes.R

@@ -96,7 +98,8 @@ create_processes <- function(
    targeted_process = function(timestep, target){


Suggested change

targeted_process = function(timestep, target){

targeted_process = function(timestep, target, prob, relative_rates){

giovannic · 2024-10-14T10:29:11Z

R/human_infection.R

+    ## capture infection rates to resolve in competing hazards
+    infection_outcome$set_rates(
+      source_humans, 
+      infection_rates)


The edit to CompetingOutcome would allow you to add prob and relative_rates here...

giovannic · 2024-10-14T10:36:16Z

/benchmark

github-actions · 2024-10-14T11:43:15Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if e2adee2 is merged into dev:

✔️large_population: 2.36m -> 2.39m [-0.39%, +2.98%]
🚀small_population: 33.2s -> 32.5s [-3.63%, -0.59%]

Further explanation regarding interpretation and methodology can be found in the documentation.
Plots and raw data are available as artifacts of the workflow run.

RJSheppard · 2024-10-14T14:50:12Z

/benchmark

giovannic

Very nearly there for me!

giovannic · 2024-10-14T15:06:41Z

R/human_infection.R

+      # infection_outcome$set_relative_rates(
+      # hypnozoites_humans, 
+      # relative_rates)


dead code

Suggested change

# infection_outcome$set_relative_rates(

# hypnozoites_humans,

# relative_rates)

giovannic · 2024-10-14T15:11:44Z

R/human_infection.R

@@ -685,7 +944,7 @@ severe_immunity <- function(age, acquired_immunity, maternal_immunity, parameter
  parameters$theta0 * (parameters$theta1 + (1 - parameters$theta1) / (
    1 + fv * (
      (acquired_immunity + maternal_immunity) / parameters$iv0) ** parameters$kv
-    )
+  )


Suggested change

)

)

giovannic · 2024-10-14T15:23:11Z

R/competing_hazards.R

    },
    execute = function(t, target){
-      private$targeted_process(t, target)
+      do.call(private$targeted_process, list(t, target, self$args))


Oh, I think it's most intuitive to call targeted_process with the extra arguments in the same way they were passed to $set_rates. i.e.

outcome$set_rates(rates, a, b)

should call the process with:

targeted_process(target, a, b) <- this version preserves the structure of my extra arguments

not

targeted_process(target, list(a, b)) <- this version turns my extra arguments into a list.

Suggested change

do.call(private$targeted_process, list(t, target, self$args))

do.call(private$targeted_process, c((t, target), self$args))

github-actions · 2024-10-14T15:56:28Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 3b833ff is merged into dev:

✔️large_population: 2.33m -> 2.36m [-0.5%, +2.73%]
✔️small_population: 33s -> 32.5s [-4.09%, +0.68%]

Further explanation regarding interpretation and methodology can be found in the documentation.
Plots and raw data are available as artifacts of the workflow run.

RJSheppard · 2024-10-15T13:58:14Z

/benchmark

giovannic

Thanks!

github-actions · 2024-10-15T15:04:43Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if 7166cba is merged into dev:

❗🐌large_population: 2.32m -> 2.36m [+0.37%, +3.28%]
🚀small_population: 32.9s -> 32.1s [-3.66%, -1.24%]

Further explanation regarding interpretation and methodology can be found in the documentation.
Plots and raw data are available as artifacts of the workflow run.

pwinskill

Thanks for the hard work @RJSheppard!
One non-critical clarification from me

pwinskill · 2024-10-22T13:42:17Z

R/human_infection.R

+#' @param renderer model renderer
+#' @param timestep current timestep
+#' @noRd
+relapse_bite_infection_hazard_resolution <- function(


I'm probably just not following, but given the name of this function I was expecting to see some referecnce to the competing hazards object within it, but can't?

The main CH resolution is between infections and disease progression. This part of the code resolves between relapses and bite infections and doesn't use the competing hazards object. It's relatively simple and discretely resolved, so doesn't need the object.

…ctions, decay of hypnozoite batches, new infections via relapse and the reseting of batch number to 0 for new births. A hypnozoite batch variable is created. The vivax model considers the number of new bites an individual receives (in contrast to whether an individual has been bitten regardless of number of bites, as currently in the P. falciparum model) to calculate the infection rates. This is now captured in bitten_humans, which is now a list that also includes n_bites_each. The vivax model does not factor the eir by the population level biting heterogeneity (although it does use the heterogeneity to distribute biting, but retains the overall total biting according to the EIR). biting_rate.R In vivax, we track bites in all individuals, not just SAU as in falciparum. This is because new hypnozoite batches can occur in any human state, regardless of whether this results in a new infection. The number of bites is considered in calculating probability of a new infection. To calculate the infection rates for vivax, we take the biting infection rates and add them to the relapse rate. We also save the relative rates in the competing hazards object to resolve these competing hazards later. A new function to resolve between bites and relapse infections is created. For each bite we increase the number of batches - with capacity for prophylaxis to be added at a future step. We also cap the total number of batches. We add a hypnozoite batch decay process. Note that this process differs from the immunity decay processes. We work in discrete outputs (you can't have a fractional batch). We also work in the probability that an individual experiences a single loss, rather than loss at the hypnozoite level. When an individual dies, we reset the hypnozoite batches back to 0. We now render number of relapses, average hypnozoite batch number and the number of individuals with hypnozoite batches. We may also output these by specified age groups. Added a test to check relapses occur as expected and expanded the p.v vignette.

… within the function. This also removed the "set_relative_rates" part of the competing hazards, and now we have args, such that the competing hazards object can still store the relative rates it needs, but retains flexibility. Calculate_infections has been rewritten as four functions: calculate_falciparum_infections, calculate_vivax_infections, treatment_prophylaxis_efficacy and pev_efficacy. Infection_outcome_process has been split into falciparum and vivax _infection_outcome_process. Tests have been adjusted to reflect these changes.

…peting hazards object and reverting to previous competing hazard tests.

RJSheppard · 2024-10-22T15:44:58Z

/benchmark

github-actions · 2024-10-22T16:50:48Z

This is how benchmark results would change (along with a 95% confidence interval in relative change) if da381b9 is merged into dev:

✔️large_population: 2.32m -> 2.33m [-0.29%, +1.99%]
✔️small_population: 33.2s -> 32.9s [-3.77%, +2.36%]

Further explanation regarding interpretation and methodology can be found in the documentation.
Plots and raw data are available as artifacts of the workflow run.

RJSheppard requested review from giovannic, pwinskill and tbreweric July 3, 2024 14:28

pwinskill requested changes Jul 16, 2024

View reviewed changes

pwinskill requested changes Jul 17, 2024

View reviewed changes

RJSheppard force-pushed the vivax_competing_hazards_main_subpatent_infections branch from c892bc3 to 29b1fed Compare September 17, 2024 17:42

RJSheppard force-pushed the vivax_competing_hazards_main_hypnozoites branch from 9c143ac to 5f5260a Compare September 18, 2024 11:44

RJSheppard force-pushed the vivax_competing_hazards_main_subpatent_infections branch from 29b1fed to 39949ff Compare September 18, 2024 12:07

RJSheppard force-pushed the vivax_competing_hazards_main_hypnozoites branch 3 times, most recently from 4b38d43 to 1916fa1 Compare September 19, 2024 12:58

RJSheppard force-pushed the vivax_competing_hazards_main_subpatent_infections branch from 39949ff to b8cce39 Compare September 19, 2024 14:15

RJSheppard force-pushed the vivax_competing_hazards_main_hypnozoites branch from 1916fa1 to b481153 Compare September 19, 2024 14:22

RJSheppard force-pushed the vivax_competing_hazards_main_subpatent_infections branch from b8cce39 to 3472e76 Compare September 20, 2024 16:22

RJSheppard force-pushed the vivax_competing_hazards_main_hypnozoites branch from b481153 to 42083f2 Compare September 20, 2024 16:27

RJSheppard force-pushed the vivax_competing_hazards_main_subpatent_infections branch 2 times, most recently from bfe4c0a to cb35f28 Compare October 11, 2024 14:36

Base automatically changed from vivax_competing_hazards_main_subpatent_infections to dev October 14, 2024 08:47

RJSheppard force-pushed the vivax_competing_hazards_main_hypnozoites branch from 42083f2 to e2adee2 Compare October 14, 2024 09:16

giovannic requested changes Oct 14, 2024

View reviewed changes

giovannic approved these changes Oct 15, 2024

View reviewed changes

pwinskill approved these changes Oct 22, 2024

View reviewed changes

RJSheppard added 3 commits October 22, 2024 16:43

Removed code redundancies. Adjusted extra arguments inputted into com…

da381b9

…peting hazards object and reverting to previous competing hazard tests.

RJSheppard force-pushed the vivax_competing_hazards_main_hypnozoites branch from 7166cba to da381b9 Compare October 22, 2024 15:43

RJSheppard merged commit a328fb3 into dev Oct 22, 2024
4 of 5 checks passed

RJSheppard deleted the vivax_competing_hazards_main_hypnozoites branch October 22, 2024 16:15

		prob = rate_to_prob(infection_outcome$rates),
		relative_rates = infection_outcome$relative_rates)

	private$targeted_process(t, target)
	do.call(private$targeted_process, c(c(t, target), self$args)

		@@ -96,7 +98,8 @@ create_processes <- function(
		targeted_process = function(timestep, target){

	# infection_outcome$set_relative_rates(
	# hypnozoites_humans,
	# relative_rates)

	do.call(private$targeted_process, list(t, target, self$args))
	do.call(private$targeted_process, c((t, target), self$args))

Add hypnozoites to the P.v model #321

Add hypnozoites to the P.v model #321

Conversation

RJSheppard commented Jul 3, 2024

pwinskill left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

pwinskill left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

RJSheppard commented Oct 14, 2024

github-actions bot commented Oct 14, 2024

giovannic left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

giovannic commented Oct 14, 2024

github-actions bot commented Oct 14, 2024

RJSheppard commented Oct 14, 2024

giovannic left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

github-actions bot commented Oct 14, 2024

RJSheppard commented Oct 15, 2024

giovannic left a comment

Choose a reason for hiding this comment

github-actions bot commented Oct 15, 2024

pwinskill left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

RJSheppard commented Oct 22, 2024

github-actions bot commented Oct 22, 2024