From 903547aeb9d8333cb36bd72bfd8aa6e1ebf92d08 Mon Sep 17 00:00:00 2001
From: RJSheppard <rjs11@ic.ac.uk>
Date: Wed, 8 May 2024 11:22:26 +0100
Subject: [PATCH] competing_hazards.R: Giovanni designed two objects:
 CompetingOutcome (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.
---
 R/biting_process.R                          | 100 +++---
 R/competing_hazards.R                       |  78 +++++
 R/disease_progression.R                     | 117 ++++---
 R/human_infection.R                         | 202 ++++++++---
 R/processes.R                               | 106 +++---
 R/render.R                                  |  54 +--
 R/utils.R                                   |  17 +-
 tests/testthat/test-biology.R               |   2 +-
 tests/testthat/test-biting-integration.R    |  33 +-
 tests/testthat/test-competing-hazards.R     | 250 ++++++++++++++
 tests/testthat/test-infection-integration.R | 364 ++++++++++++--------
 tests/testthat/test-pev.R                   |  16 +-
 12 files changed, 952 insertions(+), 387 deletions(-)
 create mode 100644 R/competing_hazards.R
 create mode 100644 tests/testthat/test-competing-hazards.R

diff --git a/R/biting_process.R b/R/biting_process.R
index 801aea20..ff78a3a3 100644
--- a/R/biting_process.R
+++ b/R/biting_process.R
@@ -15,8 +15,10 @@
 #' values (default: 1)
 #' @param mixing_index an index for this population's position in the
 #' lagged_infectivity list (default: 1)
+#' @param infection_outcome competing hazards object for infection rates
+#' @param timestep the current timestep
 #' @noRd
-create_biting_process <- function(
+biting_process <- function(
     renderer,
     solvers,
     models,
@@ -26,37 +28,38 @@ create_biting_process <- function(
     lagged_infectivity,
     lagged_eir,
     mixing = 1,
-    mixing_index = 1
+    mixing_index = 1,
+    infection_outcome,
+    timestep
 ) {
-  function(timestep) {
-    # Calculate combined EIR
-    age <- get_age(variables$birth$get_values(), timestep)
-    
-    bitten_humans <- simulate_bites(
-      renderer,
-      solvers,
-      models,
-      variables,
-      events,
-      age,
-      parameters,
-      timestep,
-      lagged_infectivity,
-      lagged_eir,
-      mixing,
-      mixing_index
-    )
-    
-    simulate_infection(
-      variables,
-      events,
-      bitten_humans,
-      age,
-      parameters,
-      timestep,
-      renderer
-    )
-  }
+
+  age <- get_age(variables$birth$get_values(), timestep)
+
+  bitten_humans <- simulate_bites(
+    renderer,
+    solvers,
+    models,
+    variables,
+    events,
+    age,
+    parameters,
+    timestep,
+    lagged_infectivity,
+    lagged_eir,
+    mixing,
+    mixing_index
+  )
+
+  simulate_infection(
+    variables,
+    events,
+    bitten_humans,
+    age,
+    parameters,
+    timestep,
+    renderer,
+    infection_outcome
+  )
 }
 
 #' @importFrom stats rpois
@@ -74,8 +77,9 @@ simulate_bites <- function(
     mixing = 1,
     mixing_index = 1
 ) {
+
   bitten_humans <- individual::Bitset$new(parameters$human_population)
-  
+
   human_infectivity <- variables$infectivity$get_values()
   if (parameters$tbv) {
     human_infectivity <- account_for_tbv(
@@ -86,21 +90,21 @@ simulate_bites <- function(
     )
   }
   renderer$render('infectivity', mean(human_infectivity), timestep)
-  
+
   # Calculate pi (the relative biting rate for each human)
   psi <- unique_biting_rate(age, parameters)
   zeta <- variables$zeta$get_values()
   .pi <- human_pi(zeta, psi)
-  
+
   # Get some indices for later
   if (parameters$individual_mosquitoes) {
     infectious_index <- variables$mosquito_state$get_index_of('Im')
     susceptible_index <- variables$mosquito_state$get_index_of('Sm')
     adult_index <- variables$mosquito_state$get_index_of('NonExistent')$not(TRUE)
   }
-  
+
   EIR <- 0
-  
+
   for (s_i in seq_along(parameters$species)) {
     species_name <- parameters$species[[s_i]]
     solver_states <- solvers[[s_i]]$get_states()
@@ -111,7 +115,7 @@ simulate_bites <- function(
     f <- blood_meal_rate(s_i, Z, parameters)
     a <- .human_blood_meal_rate(f, s_i, W, parameters)
     lambda <- effective_biting_rates(a, .pi, p_bitten)
-    
+
     if (parameters$individual_mosquitoes) {
       species_index <- variables$species$get_index_of(
         parameters$species[[s_i]]
@@ -127,10 +131,10 @@ simulate_bites <- function(
     } else {
       n_infectious <- calculate_infectious_compartmental(solver_states)
     }
-    
+
     # store the current population's EIR for later
     lagged_eir[[mixing_index]][[s_i]]$save(n_infectious * a, timestep)
-    
+
     # calculated the EIR for this timestep after mixing
     species_eir <- sum(
       vnapply(
@@ -138,7 +142,7 @@ simulate_bites <- function(
         function(l) l[[s_i]]$get(timestep - parameters$de)
       ) * mixing
     )
-    
+
     renderer$render(paste0('EIR_', species_name), species_eir, timestep)
     EIR <- EIR + species_eir
     expected_bites <- species_eir * mean(psi)
@@ -150,7 +154,7 @@ simulate_bites <- function(
         )
       }
     }
-    
+
     infectivity <- vnapply(
       lagged_infectivity,
       function(l) l$get(timestep - parameters$delay_gam)
@@ -163,7 +167,7 @@ simulate_bites <- function(
     renderer$render(paste0('FOIM_', species_name), foim, timestep)
     mu <- death_rate(f, W, Z, s_i, parameters)
     renderer$render(paste0('mu_', species_name), mu, timestep)
-    
+
     if (parameters$individual_mosquitoes) {
       # update the ODE with stats for ovoposition calculations
       aquatic_mosquito_model_update(
@@ -172,10 +176,10 @@ simulate_bites <- function(
         f,
         mu
       )
-      
+
       # update the individual mosquitoes
       susceptible_species_index <- susceptible_index$copy()$and(species_index)
-      
+
       biting_effects_individual(
         variables,
         foim,
@@ -197,7 +201,7 @@ simulate_bites <- function(
       )
     }
   }
-  
+
   renderer$render('n_bitten', bitten_humans$size(), timestep)
   bitten_humans
 }
@@ -210,9 +214,7 @@ simulate_bites <- function(
 calculate_eir <- function(species, solvers, variables, parameters, timestep) {
   a <- human_blood_meal_rate(species, variables, parameters, timestep)
   infectious <- calculate_infectious(species, solvers, variables, parameters)
-  age <- get_age(variables$birth$get_values(), timestep)
-  psi <- unique_biting_rate(age, parameters)
-  infectious * a * mean(psi)
+  infectious * a
 }
 
 effective_biting_rates <- function(a, .pi, p_bitten) {
@@ -246,7 +248,7 @@ calculate_infectious_individual <- function(
     species_index,
     parameters
 ) {
-  
+
   infectious_index$copy()$and(species_index)$size()
 }
 
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
new file mode 100644
index 00000000..ed6890a5
--- /dev/null
+++ b/R/competing_hazards.R
@@ -0,0 +1,78 @@
+## Define classes to resolve competing hazards
+CompetingOutcome <- R6::R6Class(
+  "CompetingOutcome",
+  private = list(
+    targeted_process = NULL
+  ),
+  public = list(
+    initialize = function(targeted_process, size){
+      if (!is.function(targeted_process)){
+        stop("targeted_process must be a function")
+      }
+      if (!is.numeric(size) || size <= 0){
+        stop("size must be positive integer")
+      }
+      private$targeted_process <- targeted_process
+      self$rates <- rep(0, size)
+    },
+    set_rates = function(rates){
+      self$rates <- rates
+    },
+    execute = function(t, target){
+      private$targeted_process(t, target)
+      self$rates <- rep(0, length(self$rates))
+    },
+    rates = NULL
+  )
+)
+
+CompetingHazard <- R6::R6Class(
+  "CompetingHazard",
+  private = list(
+    outcomes = list(),
+    size = NULL,
+    # RNG is passed in because mockery is not able to stub runif
+    # TODO: change when fixed
+    rng = NULL
+  ),
+  public = list(
+    initialize = function(outcomes, rng = runif){
+      if (length(outcomes) == 0){
+        stop("At least one outcome must be provided")
+      }
+      if (!all(sapply(outcomes, function(x) inherits(x, "CompetingOutcome")))){
+        stop("All outcomes must be of class CompetingOutcome")
+      }
+      private$outcomes <- outcomes
+      private$size <- length(outcomes[[1]]$rates)
+      private$rng <- rng
+    },
+    resolve = function(t){
+      event_probs <- do.call(
+        'cbind',
+        lapply(private$outcomes, function(x) x$rates)
+      )
+      occur_prob <- rowSums(event_probs)
+      occur_rng <- private$rng(private$size)
+      occurs <- occur_rng < occur_prob
+      norm_probs <- event_probs / occur_prob
+      norm_probs[is.na(norm_probs)] <- 0
+      cum_probs <- apply(norm_probs, 1, cumsum)
+      event_rng <- private$rng(private$size)
+      for(o in seq_along(private$outcomes)){
+        if (o == 1) {
+          selected <- (event_rng <= cum_probs[o,])
+        } else {
+          selected <- (event_rng > cum_probs[o - 1,]) &
+            (event_rng <= cum_probs[o,])
+        }
+        target <- individual::Bitset$new(private$size)$insert(
+          which(selected & occurs)
+        )
+        if (target$size() > 0){
+          private$outcomes[[o]]$execute(t, target)
+        }
+      }
+    }
+  )
+)
diff --git a/R/disease_progression.R b/R/disease_progression.R
index d9ac11e4..fd2cb74d 100644
--- a/R/disease_progression.R
+++ b/R/disease_progression.R
@@ -1,3 +1,71 @@
+#' @title Calculate recovery rates
+#' @description Calculates recovery rates for each individual in the population
+#' for storage in competing hazards object and future resolution
+#'
+#' @param variables the available human variables
+#' @param parameters model parameters
+#' @param recovery_outcome competing hazards object for recovery rates
+#' @noRd
+calculate_recovery_rates <- function(variables, parameters, recovery_outcome){
+  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/parameters$dt
+  recovery_outcome$set_rates(recovery_rates)
+}
+
+#' @title Disease progression outcomes (recovery)
+#' @description Following resolution of competing hazards, update state and
+#' infectivity of sampled individuals
+#'
+#' @param timestep the current timestep
+#' @param target the sampled recovering individuals
+#' @param variables the available human variables
+#' @param parameters model parameters
+#' @param renderer competing hazards object for recovery rates
+#' @noRd
+recovery_process_resolved_hazard <- function(
+    timestep,
+    target,
+    variables,
+    parameters,
+    renderer
+    ){
+
+  update_to_asymptomatic_infection(
+    variables,
+    parameters,
+    timestep,
+    variables$state$get_index_of("D")$and(target)
+  )
+
+  update_infection(
+    variables$state,
+    "U",
+    variables$infectivity,
+    parameters$cu,
+    variables$state$get_index_of("A")$and(target)
+  )
+
+  update_infection(
+    variables$state,
+    "S",
+    variables$infectivity,
+    0,
+    variables$state$get_index_of("U")$and(target)
+  )
+
+  update_infection(
+    variables$state,
+    "S",
+    variables$infectivity,
+    0,
+    variables$state$get_index_of("Tr")$and(target)
+  )
+
+}
+
 #' @title Update the state of an individual as infection events occur
 #' @description Randomly moves individuals towards the later stages of disease
 #' and updates their infectivity
@@ -18,55 +86,6 @@ update_infection <- function(
   infectivity$queue_update(new_infectivity, to_move)
 }
 
-create_progression_process <- function(
-    state,
-    from_state,
-    to_state,
-    rate,
-    infectivity,
-    new_infectivity
-) {
-  function(timestep) {
-    
-    # Retrieve the indices of all individuals in the to_move state:
-    index <- state$get_index_of(from_state)
-    
-    # If the length of rate is greater than 1 (when it's a variable):
-    if (inherits(rate, "DoubleVariable")) {
-      rate <- rate$get_values(index)
-    }
-    
-    # Sample the individuals to be moved into a new Bitset using the transition rate(s):
-    to_move <- index$sample(1/rate)
-    
-    # Update the infection status of those individuals who are moving:
-    update_infection(
-      state,
-      to_state,
-      infectivity,
-      new_infectivity,
-      to_move
-    )
-  }
-}
-
-create_asymptomatic_progression_process <- function(
-  state,
-  rate,
-  variables,
-  parameters
-  ) {
-  function(timestep) {
-    to_move <- state$get_index_of('D')$sample(1/rate)
-    update_to_asymptomatic_infection(
-      variables,
-      parameters,
-      timestep,
-      to_move
-    )
-  }
-}
-
 #' @title Modelling the progression to asymptomatic disease
 #' @description Randomly moves individuals to asymptomatic disease and
 #' calculates the infectivity for their age and immunity
diff --git a/R/human_infection.R b/R/human_infection.R
index d62b16f0..eedd4029 100644
--- a/R/human_infection.R
+++ b/R/human_infection.R
@@ -1,23 +1,26 @@
 #' @title Simulate malaria infection in humans
 #' @description
-#' 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.
 #' @param variables a list of all of the model variables
 #' @param events a list of all of the model events
 #' @param bitten_humans a bitset of bitten humans
 #' @param age of each human (timesteps)
 #' @param parameters of the model
 #' @param timestep current timestep
+#' @param renderer the model renderer object
+#' @param infection_outcome competing hazards object for infection rates
 #' @noRd
 simulate_infection <- function(
-  variables,
-  events,
-  bitten_humans,
-  age,
-  parameters,
-  timestep,
-  renderer
-  ) {
+    variables,
+    events,
+    bitten_humans,
+    age,
+    parameters,
+    timestep,
+    renderer,
+    infection_outcome
+) {
   if (bitten_humans$size() > 0) {
     boost_immunity(
       variables$ib,
@@ -29,11 +32,118 @@ simulate_infection <- function(
   }
 
   # Calculate Infected
-  infected_humans <- calculate_infections(
+  calculate_infections(
+    variables,
+    bitten_humans,
+    parameters,
+    renderer,
+    timestep,
+    infection_outcome
+  )
+
+}
+
+#' @title Calculate overall infections for bitten humans
+#' @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 parameters model parameters
+#' @param renderer model render object
+#' @param timestep current timestep
+#' @noRd
+calculate_infections <- function(
     variables,
     bitten_humans,
     parameters,
     renderer,
+    timestep,
+    infection_outcome
+) {
+  source_humans <- variables$state$get_index_of(
+    c('S', 'A', 'U'))$and(bitten_humans)
+
+  b <- blood_immunity(variables$ib$get_values(source_humans), parameters)
+
+  source_vector <- source_humans$to_vector()
+
+  # calculate prophylaxis
+  prophylaxis <- rep(0, length(source_vector))
+  drug <- variables$drug$get_values(source_vector)
+  medicated <- (drug > 0)
+  if (any(medicated)) {
+    drug <- drug[medicated]
+    drug_time <- variables$drug_time$get_values(source_vector[medicated])
+    prophylaxis[medicated] <- weibull_survival(
+      timestep - drug_time,
+      parameters$drug_prophylaxis_shape[drug],
+      parameters$drug_prophylaxis_scale[drug]
+    )
+  }
+
+  # calculate vaccine efficacy
+  vaccine_efficacy <- rep(0, length(source_vector))
+  vaccine_times <- variables$pev_timestep$get_values(source_vector)
+  vaccinated <- vaccine_times > -1
+  pev_profile <- variables$pev_profile$get_values(source_vector)
+  pev_profile <- pev_profile[vaccinated]
+  if (length(vaccinated) > 0) {
+    antibodies <- calculate_pev_antibodies(
+      timestep - vaccine_times[vaccinated],
+      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'cs')),
+      invlogit(sample_pev_param(pev_profile, parameters$pev_profiles, 'rho')),
+      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'ds')),
+      exp(sample_pev_param(pev_profile, parameters$pev_profiles, 'dl')),
+      parameters
+    )
+    vmax <- vnapply(parameters$pev_profiles, function(p) p$vmax)
+    beta <- vnapply(parameters$pev_profiles, function(p) p$beta)
+    alpha <- vnapply(parameters$pev_profiles, function(p) p$alpha)
+    vaccine_efficacy[vaccinated] <- calculate_pev_efficacy(
+      antibodies,
+      vmax[pev_profile],
+      beta[pev_profile],
+      alpha[pev_profile]
+    )
+  }
+
+  prob <- b * (1 - prophylaxis) * (1 - vaccine_efficacy)
+
+  ## Capture infection rates to resolve in competing hazards
+  infection_rates <- numeric(length = parameters$human_population)
+  infection_rates[source_vector] <- prob_to_rate(prob)
+  infection_outcome$set_rates(infection_rates)
+}
+
+#' @title Assigns infections to appropriate human states
+#' @description
+#' Updates human states and variables to represent asymptomatic/clinical/severe
+#' and treated malaria; and resulting boosts in immunity
+#' @param timestep current timestep
+#' @param infected_humans bitset of infected humans
+#' @param variables a list of all of the model variables
+#' @param renderer model render object
+#' @param parameters model parameters
+#' @param prob vector of population probabilities of infection
+#' @noRd
+infection_process_resolved_hazard <- function(
+    timestep,
+    infected_humans,
+    variables,
+    renderer,
+    parameters,
+    prob){
+
+  source_humans <- variables$state$get_index_of(values = c("S","A","U"))
+
+  incidence_renderer(
+    variables$birth,
+    renderer,
+    infected_humans,
+    source_humans,
+    prob[source_humans$to_vector()],
+    'inc_',
+    parameters$incidence_rendering_min_ages,
+    parameters$incidence_rendering_max_ages,
     timestep
   )
 
@@ -90,6 +200,7 @@ simulate_infection <- function(
   )
 }
 
+<<<<<<< HEAD
 #' @title Calculate overall infections for bitten humans
 #' @description
 #' Sample infected humans given prophylaxis and vaccination
@@ -171,6 +282,8 @@ calculate_infections <- function(
 
   infected
 }
+=======
+>>>>>>> d694b08 (Infection and recovery competing hazard resolution.)
 
 #' @title Calculate clinical infections
 #' @description
@@ -182,12 +295,12 @@ calculate_infections <- function(
 #' @param timestep current timestep
 #' @noRd
 calculate_clinical_infections <- function(
-  variables,
-  infections,
-  parameters,
-  renderer,
-  timestep
-  ) {
+    variables,
+    infections,
+    parameters,
+    renderer,
+    timestep
+) {
   ica <- variables$ica$get_values(infections)
   icm <- variables$icm$get_values(infections)
   phi <- clinical_immunity(ica, icm, parameters)
@@ -216,12 +329,12 @@ calculate_clinical_infections <- function(
 #' @param renderer model outputs
 #' @noRd
 update_severe_disease <- function(
-  timestep,
-  infections,
-  variables,
-  parameters,
-  renderer
-  ) {
+    timestep,
+    infections,
+    variables,
+    parameters,
+    renderer
+) {
   age <- get_age(variables$birth$get_values(infections), timestep)
   iva <- variables$iva$get_values(infections)
   ivm <- variables$ivm$get_values(infections)
@@ -269,11 +382,14 @@ calculate_treated <- function(
     timestep,
     renderer
 ) {
+<<<<<<< HEAD
   
   if(clinical_infections$size() == 0) {
     return(individual::Bitset$new(parameters$human_population))
   }
   
+=======
+>>>>>>> d694b08 (Infection and recovery competing hazard resolution.)
   treatment_coverages <- get_treatment_coverages(parameters, timestep)
   ft <- sum(treatment_coverages)
   
@@ -284,8 +400,12 @@ calculate_treated <- function(
   renderer$render('ft', ft, timestep)
   seek_treatment <- sample_bitset(clinical_infections, ft)
   n_treat <- seek_treatment$size()
+<<<<<<< HEAD
+=======
+
+>>>>>>> d694b08 (Infection and recovery competing hazard resolution.)
   renderer$render('n_treated', n_treat, timestep)
-  
+
   drugs <- as.numeric(parameters$clinical_treatment_drugs[
     sample.int(
       length(parameters$clinical_treatment_drugs),
@@ -380,13 +500,13 @@ calculate_treated <- function(
 #' @param parameters model parameters
 #' @noRd
 schedule_infections <- function(
-  variables,
-  clinical_infections,
-  treated,
-  infections,
-  parameters,
-  timestep
-  ) {
+    variables,
+    clinical_infections,
+    treated,
+    infections,
+    parameters,
+    timestep
+) {
   included <- treated$not(TRUE)
 
   to_infect <- clinical_infections$and(included)
@@ -418,12 +538,12 @@ schedule_infections <- function(
 # Utility functions
 # =================
 boost_immunity <- function(
-  immunity_variable,
-  exposed_index,
-  last_boosted_variable,
-  timestep,
-  delay
-  ) {
+    immunity_variable,
+    exposed_index,
+    last_boosted_variable,
+    timestep,
+    delay
+) {
   # record who can be boosted
   exposed_index_vector <- exposed_index$to_vector()
   last_boosted <- last_boosted_variable$get_values(exposed_index)
@@ -446,7 +566,7 @@ boost_immunity <- function(
 # Implemented from Winskill 2017 - Supplementary Information page 4
 clinical_immunity <- function(acquired_immunity, maternal_immunity, parameters) {
   acquired_immunity[acquired_immunity > 0] <- acquired_immunity[acquired_immunity > 0] + 0.5
-  
+
   parameters$phi0 * (
     parameters$phi1 +
       (1 - parameters$phi1) /
@@ -461,14 +581,14 @@ clinical_immunity <- function(acquired_immunity, maternal_immunity, parameters)
 # Implemented from Winskill 2017 - Supplementary Information page 5
 severe_immunity <- function(age, acquired_immunity, maternal_immunity, parameters) {
   acquired_immunity[acquired_immunity > 0] <- acquired_immunity[acquired_immunity > 0] + 0.5
-  
+
   fv <- 1 - (1 - parameters$fv0) / (
     1 + (age / parameters$av) ** parameters$gammav
   )
   parameters$theta0 * (parameters$theta1 + (1 - parameters$theta1) / (
     1 + fv * (
       (acquired_immunity + maternal_immunity) / parameters$iv0) ** parameters$kv
-    )
+  )
   )
 }
 
@@ -495,7 +615,7 @@ asymptomatic_infectivity <- function(age, immunity, parameters) {
 # Implemented from Winskill 2017 - Supplementary Information page 4
 blood_immunity <- function(ib, parameters) {
   ib[ib > 0] <- ib[ib > 0] + 0.5
-  
+
   parameters$b0 * (
     parameters$b1 +
       (1 - parameters$b1) /
diff --git a/R/processes.R b/R/processes.R
index 12f8c3db..f4c18ab3 100644
--- a/R/processes.R
+++ b/R/processes.R
@@ -30,7 +30,7 @@ create_processes <- function(
     correlations,
     lagged_eir,
     lagged_infectivity,
-    timesteps, 
+    timesteps,
     mixing = 1,
     mixing_index = 1
 ) {
@@ -62,15 +62,37 @@ create_processes <- function(
     )
   }
 
+  # =====================================================
+  # Competing Hazard Outcomes (Infections and Recoveries)
+  # =====================================================
+
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters, prob = rate_to_prob(infection_outcome$rates))
+    },
+    size = parameters$human_population
+  )
+
+  recovery_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      recovery_process_resolved_hazard(timestep, target, variables, parameters, renderer)
+    },
+    size = parameters$human_population
+  )
+
+  create_infection_recovery_hazard_process <- CompetingHazard$new(
+    outcomes = list(infection_outcome, recovery_outcome)
+  )
+
   # ==============================
   # Biting and mortality processes
   # ==============================
-  # schedule infections for humans and set last_boosted_*
+  # simulate infections for humans and set last_boosted_*
   # move mosquitoes into incubating state
   # kill mosquitoes caught in vector control
-  processes <- c(
-    processes,
-    create_biting_process(
+
+  create_infection_rates_process <- function(timestep){
+    biting_process(
       renderer,
       solvers,
       models,
@@ -80,59 +102,29 @@ create_processes <- function(
       lagged_infectivity,
       lagged_eir,
       mixing,
-      mixing_index
-    ),
-    create_asymptomatic_progression_process(
-      variables$state,
-      parameters$dd,
-      variables,
-      parameters
-    ),
-    create_progression_process(
-      variables$state,
-      'A',
-      'U',
-      parameters$da,
-      variables$infectivity,
-      parameters$cu
-    ),
-    create_progression_process(
-      variables$state,
-      'U',
-      'S',
-      parameters$du,
-      variables$infectivity,
-      0
+      mixing_index,
+      infection_outcome,
+      timestep
     )
-  )
-  
-  # =======================
-  # Antimalarial Resistance
-  # =======================
-  # Add an a new process which governs the transition from Tr to S when
-  # antimalarial resistance is simulated. The rate of transition switches
-  # from a parameter to a variable when antimalarial resistance == TRUE.
-  
-  # Assign the dt input to a separate object with the default single parameter value:
-  dt_input <- parameters$dt
-  
-  # If antimalarial resistance is switched on, assign dt variable values to the 
-  if(parameters$antimalarial_resistance) {
-    dt_input <- variables$dt
   }
-  
-  # Create the progression process for Tr --> S specifying dt_input as the rate:
-  processes <- c(
-    processes,
-    create_progression_process(
-      variables$state,
-      'Tr',
-      'S',
-      dt_input,
-      variables$infectivity,
-      0
+
+  # ===================
+  # Disease Progression
+  # ===================
+
+  create_recovery_rates_process <- function(timestep){
+    calculate_recovery_rates(
+      variables,
+      parameters,
+      recovery_outcome
     )
-  )
+  }
+
+  processes <- c(processes,
+                 create_infection_rates_process,
+                 create_recovery_rates_process,
+                 create_infection_recovery_hazard_process$resolve)
+
 
   # ===============
   # ODE integration
@@ -264,7 +256,11 @@ create_processes <- function(
     )
   }
 
+  # ======================
   # Mortality step
+  # ======================
+  # Mortality is not resolved as a competing hazard
+
   processes <- c(
     processes,
     create_mortality_process(variables, events, renderer, parameters))
diff --git a/R/render.R b/R/render.R
index 60cb5b73..3024597b 100644
--- a/R/render.R
+++ b/R/render.R
@@ -3,16 +3,16 @@ in_age_range <- function(birth, timestep, lower, upper) {
 }
 
 #' @title Render prevalence statistics
-#' 
+#'
 #' @description renders prevalence numerators and denominators for indivduals
 #' detected by microscopy and with severe malaria
-#' 
+#'
 #' @param state human infection state
 #' @param birth variable for birth of the individual
 #' @param immunity to detection
 #' @param parameters model parameters
 #' @param renderer model renderer
-#' 
+#'
 #' @noRd
 create_prevelance_renderer <- function(
   state,
@@ -41,7 +41,7 @@ create_prevelance_renderer <- function(
         paste0('n_', lower, '_', upper),
         in_age$size(),
         timestep
-      ) 
+      )
       renderer$render(
         paste0('n_detect_', lower, '_', upper),
         in_age$copy()$and(detected)$size(),
@@ -59,9 +59,9 @@ create_prevelance_renderer <- function(
 }
 
 #' @title Render incidence statistics
-#' 
+#'
 #' @description renders incidence (new for this timestep) for indivduals
-#' 
+#'
 #' @param birth variable for birth of the individual
 #' @param renderer object for model outputs
 #' @param target incidence population
@@ -71,19 +71,19 @@ create_prevelance_renderer <- function(
 #' @param lowers age bounds
 #' @param uppers age bounds
 #' @param timestep current target
-#' 
+#'
 #' @noRd
 incidence_renderer <- function(
-  birth,
-  renderer,
-  target,
-  source_pop,
-  prob,
-  prefix,
-  lowers,
-  uppers,
-  timestep
-  ) {
+    birth,
+    renderer,
+    target,
+    source_pop,
+    prob,
+    prefix,
+    lowers,
+    uppers,
+    timestep
+) {
   for (i in seq_along(lowers)) {
     lower <- lowers[[i]]
     upper <- uppers[[i]]
@@ -104,9 +104,9 @@ incidence_renderer <- function(
 }
 
 create_variable_mean_renderer_process <- function(
-  renderer,
-  names,
-  variables
+    renderer,
+    names,
+    variables
 ) {
   function(timestep) {
     for (i in seq_along(variables)) {
@@ -120,12 +120,12 @@ create_variable_mean_renderer_process <- function(
 }
 
 create_vector_count_renderer_individual <- function(
-  mosquito_state,
-  species,
-  state,
-  renderer,
-  parameters
-  ) {
+    mosquito_state,
+    species,
+    state,
+    renderer,
+    parameters
+) {
   function(timestep) {
     adult <- mosquito_state$get_index_of('NonExistent')$not(TRUE)
     for (i in seq_along(parameters$species)) {
@@ -168,7 +168,7 @@ create_age_group_renderer <- function(
         paste0('n_age_', lower, '_', upper),
         in_age$size(),
         timestep
-      ) 
+      )
     }
   }
 }
diff --git a/R/utils.R b/R/utils.R
index 2262c939..02c4d265 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -2,7 +2,7 @@ vnapply <- function(X, FUN, ...) vapply(X, FUN, ..., numeric(1))
 vlapply <- function(X, FUN, ...) vapply(X, FUN, ..., logical(1))
 vcapply <- function(X, FUN, ...) vapply(X, FUN, ..., character(1))
 
-#' @importFrom stats rbinom 
+#' @importFrom stats rbinom
 bernoulli <- function(size, p) sample.int(size, rbinom(1, size, min(p, 1)))
 
 sample_bitset <- function(b, rate) {
@@ -88,3 +88,18 @@ RandomState <- R6::R6Class(
     }
   )
 )
+
+#'@title Convert probability to a rate
+#'@param prob probability
+#'@noRd
+prob_to_rate <- function(prob){
+  -log(1 - prob)
+}
+
+
+#'@title Convert rate to a probability
+#'@param rate rate
+#'@noRd
+rate_to_prob <- function(rate){
+  1- exp(-rate)
+}
diff --git a/tests/testthat/test-biology.R b/tests/testthat/test-biology.R
index c739f7e7..eb04116b 100644
--- a/tests/testthat/test-biology.R
+++ b/tests/testthat/test-biology.R
@@ -92,7 +92,7 @@ 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()), 1.)
+      calculate_foim(a, sum(.pi * variables$infectivity$get_values()), mixing = 1)
     }
   )
   expect_equal(
diff --git a/tests/testthat/test-biting-integration.R b/tests/testthat/test-biting-integration.R
index 7da79b66..65ffc011 100644
--- a/tests/testthat/test-biting-integration.R
+++ b/tests/testthat/test-biting-integration.R
@@ -14,15 +14,11 @@ test_that('biting_process integrates mosquito effects and human infection', {
   models <- parameterise_mosquito_models(parameters, timestep)
   solvers <- parameterise_solvers(models, parameters)
 
-  biting_process <- create_biting_process(
-    renderer,
-    solvers,
-    models,
-    variables,
-    events,
-    parameters,
-    lagged_foim,
-    lagged_eir
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
   )
 
   bitten <- individual::Bitset$new(parameters$human_population)
@@ -31,7 +27,21 @@ test_that('biting_process integrates mosquito effects and human infection', {
 
   mockery::stub(biting_process, 'simulate_bites', bites_mock)
   mockery::stub(biting_process, 'simulate_infection', infection_mock)
-  biting_process(timestep)
+
+  biting_process(
+    renderer,
+    solvers,
+    models,
+    variables,
+    events,
+    parameters,
+    lagged_foim,
+    lagged_eir,
+    mixing = 1,
+    mixing_index = 1,
+    infection_outcome,
+    timestep
+  )
 
   mockery::expect_args(
     bites_mock,
@@ -59,7 +69,8 @@ test_that('biting_process integrates mosquito effects and human infection', {
     age,
     parameters,
     timestep,
-    renderer
+    renderer,
+    infection_outcome
   )
 })
 
diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R
new file mode 100644
index 00000000..13738f53
--- /dev/null
+++ b/tests/testthat/test-competing-hazards.R
@@ -0,0 +1,250 @@
+# ==============================
+# Test the CompetingHazard class
+# ==============================
+
+test_that('hazard resolves two normalised outcomes when all events occur', {
+
+    size <- 4
+    outcome_1_process <- mockery::mock()
+    outcome_1 <- CompetingOutcome$new(
+      targeted_process = outcome_1_process,
+      size = size
+    )
+    outcome_2_process <- mockery::mock()
+    outcome_2 <- CompetingOutcome$new(
+      targeted_process = outcome_2_process,
+      size = size
+    )
+
+    hazard <- CompetingHazard$new(
+      outcomes = list(outcome_1, outcome_2),
+      rng = mockery::mock(
+        c(0, 0, 0, 0), # all events occur
+        c(.05, .3, .2, .5) # event_rng
+      )
+    )
+
+    outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4))
+    outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6))
+
+    hazard$resolve(0)
+
+    mockery::expect_args(
+      outcome_1_process,
+      1,
+      0,
+      individual::Bitset$new(size)$insert(c(1, 3))
+    )
+    mockery::expect_args(
+      outcome_2_process,
+      1,
+      0,
+      individual::Bitset$new(size)$insert(c(2, 4))
+    )
+})
+
+test_that('hazard resolves two unnormalised outcomes when all events occur', {
+
+  size <- 4
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    outcomes = list(outcome_1, outcome_2),
+    rng = mockery::mock(
+      c(0, 0, 0, 0), # all events occur
+      c(.05, .3, .2, .5) # event_rng
+    )
+  )
+
+  outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2)
+  outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2)
+
+  hazard$resolve(0)
+
+  mockery::expect_args(
+    outcome_1_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(c(1, 3))
+  )
+  mockery::expect_args(
+    outcome_2_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(c(2, 4))
+  )
+})
+
+test_that('hazard resolves two outcomes when some events occur', {
+
+  size <- 4
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    outcomes = list(outcome_1, outcome_2),
+    rng = mockery::mock(
+      c(0, 0, 1, 1), # some events occur
+      c(.05, .3, .2, .5), # event_rng
+    )
+  )
+
+  outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2)
+  outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2)
+
+  hazard$resolve(0)
+
+  mockery::expect_args(
+    outcome_1_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(1)
+  )
+  mockery::expect_args(
+    outcome_2_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(2)
+  )
+})
+
+
+test_that('hazard resolves when no rates are set', {
+
+  size <- 4
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    outcomes = list(outcome_1, outcome_2)
+  )
+
+  hazard$resolve(0)
+
+  mockery::expect_called(
+    outcome_1_process,
+    0
+  )
+  mockery::expect_called(
+    outcome_2_process,
+    0
+  )
+})
+
+test_that('hazard resolves when rates are set to one', {
+
+  size <- 4
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    outcomes = list(outcome_1, outcome_2)
+  )
+
+  outcome_1$set_rates(rep(1, size))
+  hazard$resolve(0)
+
+  mockery::expect_args(
+    outcome_1_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(c(1, 2, 3, 4))
+  )
+  mockery::expect_called(
+    outcome_2_process,
+    0
+  )
+})
+
+test_that('hazard resolves three outcomes', {
+
+  size <- 4
+  outcome_1_process <- mockery::mock()
+  outcome_1 <- CompetingOutcome$new(
+    targeted_process = outcome_1_process,
+    size = size
+  )
+
+  outcome_2_process <- mockery::mock()
+  outcome_2 <- CompetingOutcome$new(
+    targeted_process = outcome_2_process,
+    size = size
+  )
+
+  outcome_3_process <- mockery::mock()
+  outcome_3 <- CompetingOutcome$new(
+    targeted_process = outcome_3_process,
+    size = size
+  )
+
+  hazard <- CompetingHazard$new(
+    outcomes = list(outcome_1, outcome_2, outcome_3),
+    rng = mockery::mock(
+      c(0, 0, 0, 0), # all events occur
+      c(.04, .3, .14, .6), # event_rng
+    )
+  )
+
+  outcome_1$set_rates(c(0.1, 0.2, 0.3, 0.4) / 2)
+  outcome_2$set_rates(c(0.9, 0.8, 0.7, 0.6) / 2)
+  outcome_3$set_rates(c(0.5, 0.5, 0.5, 0.5))
+
+  hazard$resolve(0)
+
+  mockery::expect_args(
+    outcome_1_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(c(1, 3))
+  )
+  mockery::expect_args(
+    outcome_2_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(2)
+  )
+  mockery::expect_args(
+    outcome_3_process,
+    1,
+    0,
+    individual::Bitset$new(size)$insert(4)
+  )
+})
+
diff --git a/tests/testthat/test-infection-integration.R b/tests/testthat/test-infection-integration.R
index 83d7def9..2e9d5c9c 100644
--- a/tests/testthat/test-infection-integration.R
+++ b/tests/testthat/test-infection-integration.R
@@ -6,7 +6,7 @@ test_that('simulate_infection integrates different types of infection and schedu
   ))
   events <- create_events(parameters)
   renderer <- mock_render(timestep)
-  
+
   age <- c(20, 24, 5, 39, 20, 24, 5, 39) * 365
   immunity <- c(.2, .3, .5, .9, .2, .3, .5, .9)
   asymptomatics <- mockery::mock()
@@ -15,27 +15,22 @@ test_that('simulate_infection integrates different types of infection and schedu
     id = individual::DoubleVariable$new(immunity),
     state = list(get_index_of = mockery::mock(asymptomatics))
   )
-  
+
   bitten <- individual::Bitset$new(population)$insert(c(1, 3, 5, 7))
   boost_immunity_mock <- mockery::mock()
   infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
   infection_mock <- mockery::mock(infected)
-  clinical <- individual::Bitset$new(population)$insert(c(1, 3))
-  clinical_infection_mock <- mockery::mock(clinical)
-  severe <- individual::Bitset$new(population)$insert(1)
-  severe_infection_mock <- mockery::mock(severe)
-  treated <- individual::Bitset$new(population)$insert(3)
-  treated_mock <- mockery::mock(treated)
-  schedule_mock <- mockery::mock()
-  
+
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
+  )
+
   mockery::stub(simulate_infection, 'boost_immunity', boost_immunity_mock)
   mockery::stub(simulate_infection, 'calculate_infections', infection_mock)
-  mockery::stub(simulate_infection, 'calculate_clinical_infections', clinical_infection_mock)
-  mockery::stub(simulate_infection, 'update_severe_disease', severe_infection_mock)
-  mockery::stub(simulate_infection, 'calculate_treated', treated_mock)
-  mockery::stub(simulate_infection, 'schedule_infections', schedule_mock)
-  mockery::stub(simulate_infection, 'incidence_renderer', mockery::mock())
-  mockery::stub(simulate_infection, 'clinical_incidence_renderer', mockery::mock())
+
   simulate_infection(
     variables,
     events,
@@ -43,9 +38,10 @@ test_that('simulate_infection integrates different types of infection and schedu
     age,
     parameters,
     timestep,
-    renderer
+    renderer,
+    infection_outcome
   )
-  
+
   mockery::expect_args(
     boost_immunity_mock,
     1,
@@ -55,7 +51,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     5,
     parameters$ub
   )
-  
+
   mockery::expect_args(
     infection_mock,
     1,
@@ -63,9 +59,69 @@ test_that('simulate_infection integrates different types of infection and schedu
     bitten,
     parameters,
     renderer,
-    timestep
+    timestep,
+    infection_outcome
+  )
+})
+
+
+test_that('simulate_infection integrates different types of infection and scheduling', {
+  population <- 8
+  timestep <- 5
+  parameters <- get_parameters(list(
+    human_population = population
+  ))
+  events <- create_events(parameters)
+  renderer <- mock_render(timestep)
+
+  age <- c(20, 24, 5, 39, 20, 24, 5, 39) * 365
+  immunity <- c(.2, .3, .5, .9, .2, .3, .5, .9)
+  # asymptomatics <- mockery::mock()
+  variables <- list(
+    ib = individual::DoubleVariable$new(immunity),
+    id = individual::DoubleVariable$new(immunity),
+    # state = list(get_index_of = mockery::mock(asymptomatics, cycle = T))
+    state = individual::CategoricalVariable$new(categories = c("S","A","U","D","Tr"), initial_values = rep("S", population))#list(get_index_of = mockery::mock(asymptomatics, cycle = T))
+  )
+  prob <- rep(0.5,population)
+
+  infected <- individual::Bitset$new(population)$insert(c(1, 3, 5))
+  clinical <- individual::Bitset$new(population)$insert(c(1, 3))
+  clinical_infection_mock <- mockery::mock(clinical)
+  boost_immunity_mock <- mockery::mock()
+  severe <- individual::Bitset$new(population)$insert(1)
+  severe_infection_mock <- mockery::mock(severe)
+  treated <- individual::Bitset$new(population)$insert(3)
+  treated_mock <- mockery::mock(treated)
+  schedule_mock <- mockery::mock()
+
+
+  mockery::stub(infection_process_resolved_hazard, 'incidence_renderer', mockery::mock())
+  mockery::stub(infection_process_resolved_hazard, 'boost_immunity', boost_immunity_mock)
+  mockery::stub(infection_process_resolved_hazard, 'calculate_clinical_infections', clinical_infection_mock)
+  mockery::stub(infection_process_resolved_hazard, 'update_severe_disease', severe_infection_mock)
+  mockery::stub(infection_process_resolved_hazard, 'calculate_treated', treated_mock)
+  mockery::stub(infection_process_resolved_hazard, 'schedule_infections', schedule_mock)
+
+  infection_process_resolved_hazard(
+    timestep,
+    infected,
+    variables,
+    renderer,
+    parameters,
+    prob)
+
+
+  mockery::expect_args(
+    boost_immunity_mock,
+    1,
+    variables$ica,
+    infected,
+    variables$last_boosted_ica,
+    5,
+    parameters$uc
   )
-  
+
   mockery::expect_args(
     clinical_infection_mock,
     1,
@@ -75,7 +131,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     renderer,
     timestep
   )
-  
+
   mockery::expect_args(
     severe_infection_mock,
     1,
@@ -85,7 +141,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     parameters,
     renderer
   )
-  
+
   mockery::expect_args(
     treated_mock,
     1,
@@ -95,7 +151,7 @@ test_that('simulate_infection integrates different types of infection and schedu
     timestep,
     renderer
   )
-  
+
   mockery::expect_args(
     schedule_mock,
     1,
@@ -125,7 +181,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     booster_coverage = matrix(1),
     booster_profile = list(rtss_booster_profile)
   )
-  
+
   variables <- list(
     state = individual::CategoricalVariable$new(
       c('D', 'S', 'A', 'U', 'Tr'),
@@ -137,7 +193,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     pev_profile = individual::IntegerVariable$new(c(-1, 1, 2, -1)),
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
-  
+
   immunity_mock <- mockery::mock(c(.2, .3, .4))
   weibull_mock <- mockery::mock(.2)
   vaccine_antibodies_mock <- mockery::mock(c(2, 3))
@@ -148,7 +204,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
   mockery::stub(calculate_infections, 'calculate_pev_antibodies', vaccine_antibodies_mock)
   mockery::stub(calculate_infections, 'calculate_pev_efficacy', vaccine_efficacy_mock)
   mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
-  
+
   # remove randomness from vaccine parameters
   mockery::stub(
     calculate_infections,
@@ -158,19 +214,27 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     },
     depth = 4
   )
-  
+
   bitten_humans <- individual::Bitset$new(4)$insert(c(1, 2, 3, 4))
-  
+
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
+  )
+
   infections <- calculate_infections(
     variables,
-    bitten_humans, 
+    bitten_humans,
     parameters,
     mock_render(timestep),
-    timestep
+    timestep,
+    infection_outcome = infection_outcome
   )
-  
-  expect_equal(infections$to_vector(), 3)
-  
+
+  expect_equal(sum(infections!=0), 3)
+
   mockery::expect_args(immunity_mock, 1, c(.3, .5, .9), parameters)
   mockery::expect_args(
     weibull_mock,
@@ -179,7 +243,7 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     parameters$drug_prophylaxis_shape[[2]],
     parameters$drug_prophylaxis_scale[[2]]
   )
-  
+
   mockery::expect_args(
     vaccine_antibodies_mock,
     1,
@@ -198,18 +262,12 @@ test_that('calculate_infections works various combinations of drug and vaccinati
     c(rtss_profile$beta, rtss_booster_profile$beta),
     c(rtss_profile$alpha, rtss_booster_profile$alpha)
   )
-  mockery::expect_args(
-    bernoulli_mock,
-    1,
-    c(.2 * .8 * .8, .3 * .7, .4)
-  )
-  
 })
 
 test_that('calculate_clinical_infections correctly samples clinically infected', {
   timestep <- 5
   parameters <- get_parameters()
-  
+
   variables <- list(
     ica = individual::DoubleVariable$new(c(.2, .3, .5, .9)),
     icm = individual::DoubleVariable$new(c(.2, .3, .5, .9)),
@@ -217,25 +275,25 @@ test_that('calculate_clinical_infections correctly samples clinically infected',
     last_boosted_ica = individual::DoubleVariable$new(c(-1, -1, 1, -1)),
     last_boosted_id = individual::DoubleVariable$new(c(-1, -1, 1, -1))
   )
-  
+
   immunity_mock <- mockery::mock(c(.2, .3, .4))
   boost_mock <- mockery::mock()
   mockery::stub(calculate_clinical_infections, 'boost_immunity', boost_mock)
-  
+
   mockery::stub(calculate_clinical_infections, 'clinical_immunity', immunity_mock)
   bernoulli_mock <- mockery::mock(c(1, 3))
   mockery::stub(calculate_clinical_infections, 'bernoulli_multi_p', bernoulli_mock)
-  
+
   infections <- individual::Bitset$new(4)$insert(c(2, 3, 4))
-  
+
   clinical_infections <- calculate_clinical_infections(
     variables,
     infections,
     parameters
   )
-  
+
   expect_equal(clinical_infections$to_vector(), c(2, 4))
-  
+
   mockery::expect_args(
     immunity_mock,
     1,
@@ -243,7 +301,7 @@ test_that('calculate_clinical_infections correctly samples clinically infected',
     c(.3, .5, .9),
     parameters
   )
-  
+
   mockery::expect_args(
     bernoulli_mock,
     1,
@@ -264,10 +322,10 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     drug = list(queue_update = mockery::mock()),
     drug_time = list(queue_update = mockery::mock())
   )
-  
+
   recovery_mock <- mockery::mock()
   mockery::stub(calculate_treated, 'recovery$schedule', recovery_mock)
-  
+
   seek_treatment <- individual::Bitset$new(4)$insert(c(1, 2, 4))
   mockery::stub(
     calculate_treated,
@@ -279,7 +337,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
   bernoulli_mock <- mockery::mock(c(1, 3))
   mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock)
   mockery::stub(calculate_treated, 'log_uniform', mockery::mock(c(3, 4)))
-  
+
   clinical_infections <- individual::Bitset$new(4)
   clinical_infections$insert(c(1, 2, 3, 4))
   calculate_treated(
@@ -289,7 +347,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     timestep,
     mock_render(timestep)
   )
-  
+
   mockery::expect_args(
     sample_mock,
     1,
@@ -298,13 +356,13 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     c(.25, .25),
     TRUE
   )
-  
+
   mockery::expect_args(
     bernoulli_mock,
     1,
     parameters$drug_efficacy[c(2, 1, 1, 1)]
   )
-  
+
   expect_bitset_update(variables$state$queue_update, 'Tr', c(1, 4))
   expect_bitset_update(
     variables$infectivity$queue_update,
@@ -316,7 +374,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
 })
 
 test_that('calculate_treated correctly samples treated and updates the drug state when resistance set', {
-  
+
   parameters <- get_parameters()
   parameters <- set_drugs(parameters = parameters, drugs = list(AL_params, SP_AQ_params))
   parameters <- set_clinical_treatment(parameters = parameters, drug = 1, timesteps = 1, coverages = 0.25)
@@ -330,7 +388,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
                                             early_treatment_failure_probability = 0.2,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 10)
   parameters <- set_antimalarial_resistance(parameters = parameters,
                                             drug = 2,
@@ -341,9 +399,9 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
                                             early_treatment_failure_probability = 0.9,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 15)
-  
+
   clinical_infections <- individual::Bitset$new(20)$insert(1:20)
   timestep <- 5
   events <- create_events(parameters)
@@ -355,22 +413,22 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     dt = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = timestep)
-  
+
   # Set up seek_treatment mock and instruct calculate_treated() to return it when sample_bitset() called:
   seek_treatment <- individual::Bitset$new(20)$insert(c(1:10))
   seek_treatment_mock <- mockery::mock(seek_treatment)
   mockery::stub(where = calculate_treated, what = 'sample_bitset', how = seek_treatment_mock)
-  
+
   # Set up drugs mock and instruct it to return it when sample.int() called:
   mock_drugs <- mockery::mock(c(2, 1, 1, 1, 2, 2, 2, 1, 2, 1))
   mockery::stub(calculate_treated, 'sample.int', mock_drugs)
-  
+
   # Set up bernoulli mock and instruct calculate_treated to return it when bernoulli_multi_p() called:
   bernoulli_mock <- mockery::mock(c(1, 2, 3, 4, 5, 6, 7, 8, 9),
                                   c(1, 2, 3, 4, 5, 6, 7),
                                   c(1))
   mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock)
-  
+
   calculate_treated(
     variables,
     clinical_infections,
@@ -378,7 +436,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     timestep,
     renderer
   )
-  
+
   mockery::expect_args(
     mock_drugs,
     1,
@@ -387,38 +445,38 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     c(.25, .25),
     TRUE
   )
-  
+
   mockery::expect_args(
     seek_treatment_mock,
     1,
     clinical_infections,
     0.5
   )
-  
+
   mockery::expect_args(
     bernoulli_mock,
     1,
     c(0.9, 0.95, 0.95, 0.95, 0.9, 0.9, 0.9, 0.95, 0.9, 0.95)
   )
-  
+
   mockery::expect_args(
     bernoulli_mock,
     2,
     c(0.73, 0.9, 0.9, 0.9, 0.73, 0.73, 0.73, 0.9, 0.73)
   )
-  
+
   mockery::expect_args(
     bernoulli_mock,
     2,
     1 - (unlist(parameters$artemisinin_resistance_proportion[c(2, 1, 1, 1, 2, 2, 2, 1, 2)]) * unlist(parameters$early_treatment_failure_probability[c(2, 1, 1, 1, 2, 2, 2, 1, 2)]))
   )
-  
+
   mockery::expect_args(
     bernoulli_mock,
     3,
     unlist(parameters$artemisinin_resistance_proportion[c(2, 1, 1, 1, 2, 2, 2)]) * unlist(parameters$slow_parasite_clearance_probability[c(2, 1, 1, 1, 2, 2, 2)])
   )
-  
+
   expect_bitset_update(variables$state$queue_update, 'Tr', c(1, 2, 3, 4, 5, 6, 7))
   expect_bitset_update(
     variables$infectivity$queue_update,
@@ -432,7 +490,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
 })
 
 test_that('calculate_treated correctly samples treated and updates the drug state when resistance not set for all drugs', {
-  
+
   # Establish the parameters
   parameters <- get_parameters()
   parameters <- set_drugs(parameters = parameters, drugs = list(AL_params, SP_AQ_params))
@@ -447,18 +505,18 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
                                             early_treatment_failure_probability = 0.3,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 20)
-  
+
   # Establish Bitset of clinically infected individuals
   clinical_infections <- individual::Bitset$new(20)$insert(1:20)
-  
+
   # Set the timestep to 5:
   timestep <- 5
-  
+
   # Establish the events:
   events <- create_events(parameters)
-  
+
   # Establish list of variables used in calculate_treated() using mocks:
   variables <- list(
     state = list(queue_update = mockery::mock()),
@@ -467,30 +525,30 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     drug_time = list(queue_update = mockery::mock()),
     dt = list(queue_update = mockery::mock())
   )
-  
+
   # Create a Bitset of individuals seeking treatment individuals:
   seek_treatment <- individual::Bitset$new(20)$insert(c(1:10))
-  
+
   # Create a mock of seek_treatment:
   seek_treatment_mock <- mockery::mock(seek_treatment)
-  
+
   # Specify that, when calculate_treated() calls sample_bitset(), return the seek_treatment_mock:
   mockery::stub(where = calculate_treated, what = 'sample_bitset', how = seek_treatment_mock)
-  
+
   # Create a mock_drugs object (5 of each drug):
   mock_drugs <- mockery::mock(c(2, 1, 1, 1, 2, 2, 2, 1, 2, 1))
-  
+
   # Specify that when calculate_treated() calls sample.int(), it returns mock_drugs:
   mockery::stub(calculate_treated, 'sample.int', mock_drugs)
-  
+
   # Create a bernoulli_mock of i) individuals susceptible, and ii) individuals successfully treated:
   bernoulli_mock <- mockery::mock(c(1, 2, 3, 4, 5, 6, 7, 8, 9),
                                   c(1, 2, 3, 4, 5, 6, 7),
                                   c(1))
-  
+
   # Specify that when calculate_treated() calls bernoulli_multi_p() it returns the bernoulli_mock:
   mockery::stub(calculate_treated, 'bernoulli_multi_p', bernoulli_mock)
-  
+
   # Run the calculate_treated() function now the mocks and stubs are established:
   calculate_treated(
     variables,
@@ -499,7 +557,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     timestep,
     mock_render(timestep)
   )
-  
+
   # Check that mock_drugs was called only once, and that the arguments used in the function call
   # mock_drugs() was used in (sample.int()) match those expected:
   mockery::expect_args(
@@ -510,7 +568,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     c(.25, .25),
     TRUE
   )
-  
+
   # Check that seek_treatment_mock was called only once, and that the arguments used in the function
   # call mock_drugs() was used in (sample_bitset()) match those expected:
   mockery::expect_args(
@@ -519,7 +577,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     clinical_infections,
     0.5
   )
-  
+
   # Check that the first time bernoulli_mock was called the arguments used in the function
   # call bernoulli_mock was involved in (bernoulli_multi_p()) match those expected:
   mockery::expect_args(
@@ -527,7 +585,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     1,
     parameters$drug_efficacy[c(2, 1, 1, 1, 2, 2, 2, 1, 2, 1)]
   )
-  
+
   # Check that the secnd time bernoulli_mock was called (bernoulli_multi_p()) the arguments used in
   # the function it was called in are as expected:
   mockery::expect_args(
@@ -535,14 +593,14 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     2,
     c(0.76, 1, 1, 1, 0.76, 0.76, 0.76, 1, 0.76)
   )
-  
+
   # Check that update queued that updates the state of successfully treated individuals to "Tr"
   expect_bitset_update(
     variables$state$queue_update,
     'Tr',
     c(1, 2, 3, 4, 5, 6, 7)
   )
-  
+
   # Check that update queued that updates the infectivity of successfully treated individuals to "Tr"
   # to their new infectivity (drug concentration x infectivity of "D" compartment)
   expect_bitset_update(
@@ -550,7 +608,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     parameters$cd * parameters$drug_rel_c[c(2, 1, 1, 1, 2, 2, 2)],
     c(1, 2, 3, 4, 5, 6, 7)
   )
-  
+
   # Check that update queued that updates the drug of successfully treated individuals to the drug
   # they took:
   expect_bitset_update(
@@ -558,7 +616,7 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     c(2, 1, 1, 1, 2, 2, 2),
     c(1, 2, 3, 4, 5, 6, 7)
   )
-  
+
   # Check that update queued that updates the drug time of successfully treated individuals to the
   # simulated/mocked time step (5)
   expect_bitset_update(
@@ -566,54 +624,54 @@ test_that('calculate_treated correctly samples treated and updates the drug stat
     5,
     c(1, 2, 3, 4, 5, 6, 7)
   )
-  
+
   # Check that update queued for dt for the non-slow parasite clearance individuals is correct:
   expect_bitset_update(
     variables$dt$queue_update,
     parameters$dt,
-    c(2, 3, 4, 5, 6, 7), 
+    c(2, 3, 4, 5, 6, 7),
     1)
-  
+
   # Check that update queued for dt for the slow parasite clearance individuals is correct:
   expect_bitset_update(
     variables$dt$queue_update,
     unlist(parameters$dt_slow_parasite_clearance),
-    c(1), 
+    c(1),
     2)
-  
+
 })
 
 test_that('schedule_infections correctly schedules new infections', {
   parameters <- get_parameters(list(human_population = 20))
   variables <- create_variables(parameters)
-  
+
   infections <- individual::Bitset$new(20)$insert(1:20)
   clinical_infections <- individual::Bitset$new(20)$insert(5:15)
   treated <- individual::Bitset$new(20)$insert(7:12)
-  
+
   infection_mock <- mockery::mock()
   asymp_mock <- mockery::mock()
-  
+
   mockery::stub(schedule_infections, 'update_infection', infection_mock)
   mockery::stub(schedule_infections, 'update_to_asymptomatic_infection', asymp_mock)
-  
+
   schedule_infections(
     variables,
     clinical_infections,
     treated,
     infections,
     parameters,
-    42 
+    42
   )
-  
+
   actual_infected <- mockery::mock_args(infection_mock)[[1]][[5]]$to_vector()
   actual_asymp_infected <- mockery::mock_args(asymp_mock)[[1]][[4]]$to_vector()
-  
+
   expect_equal(
     actual_infected,
     c(5, 6, 13, 14, 15)
   )
-  
+
   expect_equal(
     actual_asymp_infected,
     c(1, 2, 3, 4, 16, 17, 18, 19, 20)
@@ -625,7 +683,7 @@ test_that('prophylaxis is considered for medicated humans', {
   parameters <- set_drugs(parameters, list(AL_params, DHA_PQP_params))
   events <- create_events(parameters)
   timestep <- 50
-  
+
   variables = list(
     state = individual::CategoricalVariable$new(
       c('D', 'S', 'A', 'U', 'Tr'),
@@ -637,21 +695,29 @@ test_that('prophylaxis is considered for medicated humans', {
     pev_profile = individual::IntegerVariable$new(c(-1, -1, -1, -1)),
     ib = individual::DoubleVariable$new(c(.2, .3, .5, .9))
   )
-  
+
   bitten <- individual::Bitset$new(4)$insert(seq(4))
   m <- mockery::mock(seq(3))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', m)
-  
-  calculate_infections(
+
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
+  )
+
+  infection_rates <- calculate_infections(
     variables,
     bitten,
     parameters,
     mock_render(timestep),
-    timestep
+    timestep,
+    infection_outcome
   )
-  
+
   expect_equal(
-    mockery::mock_args(m)[[1]][[1]],
+    rate_to_prob(infection_rates[infection_rates!=0]),
     c(2.491951e-07, 2.384032e-01, 5.899334e-01),
     tolerance = 1e-3
   )
@@ -661,17 +727,17 @@ test_that('boost_immunity respects the delay period', {
   level <- c(2.4, 1.2, 0., 4.)
   immunity <- individual::DoubleVariable$new(level)
   last_boosted <- individual::DoubleVariable$new(c(11, 5, 1, 13))
-  
+
   level_mock <- mockery::mock()
   mockery::stub(boost_immunity, 'immunity_variable$queue_update', level_mock)
-  
+
   last_mock <- mockery::mock()
   mockery::stub(boost_immunity, 'last_boosted_variable$queue_update', last_mock)
-  
+
   index <- individual::Bitset$new(4)$insert(seq(4))
   timestep <- 15
   delay <- 4
-  
+
   boost_immunity(
     immunity,
     index,
@@ -679,14 +745,14 @@ test_that('boost_immunity respects the delay period', {
     timestep,
     delay
   )
-  
+
   mockery::expect_args(
     level_mock,
     1,
     c(3.4, 2.2, 1),
     seq(3)
   )
-  
+
   mockery::expect_args(
     last_mock,
     1,
@@ -699,18 +765,18 @@ test_that('boost_immunity respects the delay period', {
   level <- c(2.4, 1.2, 0., 4., 0.)
   immunity <- individual::DoubleVariable$new(level)
   last_boosted <- individual::DoubleVariable$new(c(11, 5, 1, 13, -1))
-  
+
   index <- individual::Bitset$new(5)
   index$insert(seq(5))
   timestep <- 15
   delay <- 4
-  
+
   level_mock <- mockery::mock()
   mockery::stub(boost_immunity, 'immunity_variable$queue_update', level_mock)
-  
+
   last_mock <- mockery::mock()
   mockery::stub(boost_immunity, 'last_boosted_variable$queue_update', last_mock)
-  
+
   boost_immunity(
     immunity,
     index,
@@ -718,14 +784,14 @@ test_that('boost_immunity respects the delay period', {
     timestep,
     delay
   )
-  
+
   mockery::expect_args(
     level_mock,
     1,
     c(3.4, 2.2, 1, 1),
     c(seq(3), 5)
   )
-  
+
   mockery::expect_args(
     last_mock,
     1,
@@ -738,18 +804,18 @@ test_that('boost_immunity does not update when there is no-one to update', {
   level <- c(2.4, 1.2, 0., 4., 0.)
   immunity <- individual::DoubleVariable$new(level)
   last_boosted <- individual::DoubleVariable$new(c(11, 5, 1, 13, -1))
-  
+
   index <- individual::Bitset$new(5)
   index$insert(seq(5))
   timestep <- 15
   delay <- 4
-  
+
   level_mock <- mockery::mock()
   mockery::stub(boost_immunity, 'immunity$queue_update', level_mock)
-  
+
   last_mock <- mockery::mock()
   mockery::stub(boost_immunity, 'last_boosted$queue_update', last_mock)
-  
+
   boost_immunity(
     immunity,
     index,
@@ -771,11 +837,11 @@ test_that('update_severe_disease renders with no infections', {
     severe_incidence_rendering_max_ages = 5 * 365
   ))
   variables <- create_variables(parameters)
-  
+
   render_function <- mockery::mock()
   mockery::stub(update_severe_disease, 'incidence_renderer', render_function)
   empty <- individual::Bitset$new(population)
-  
+
   update_severe_disease(
     timestep,
     empty,
@@ -783,7 +849,7 @@ test_that('update_severe_disease renders with no infections', {
     parameters,
     renderer
   )
-  
+
   mockery::expect_args(
     render_function,
     1,
@@ -812,7 +878,7 @@ test_that('calculate_treated returns empty Bitset when there is no clinical trea
                                             early_treatment_failure_probability = 0.2,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 10)
   clinical_infections <- individual::Bitset$new(20)$insert(1:20)
   timestep <- 5
@@ -824,13 +890,13 @@ test_that('calculate_treated returns empty Bitset when there is no clinical trea
     drug_time = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = 10)
-  
+
   treated <- calculate_treated(variables = variables,
                                clinical_infections = clinical_infections,
                                parameters = parameters,
                                timestep = timestep,
                                renderer = renderer)
-  
+
   expect_identical(object = treated$size(), expected = 0, info = "Error: calculate_treated() returning non-zero number of treated individuals
                  in the absence of clinical treatment")
 })
@@ -848,7 +914,7 @@ test_that('calculate_treated returns empty Bitset when the clinically_infected i
                                             early_treatment_failure_probability = 0.2,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 10)
   clinical_infections <- individual::Bitset$new(20)
   timestep <- 5
@@ -860,13 +926,13 @@ test_that('calculate_treated returns empty Bitset when the clinically_infected i
     drug_time = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = 10)
-  
+
   treated <- calculate_treated(variables = variables,
                                clinical_infections = clinical_infections,
                                parameters = parameters,
                                timestep = timestep,
                                renderer = renderer)
-  
+
   expect_identical(object = treated$size(), expected = 0, info = "Error: calculate_treated() returning non-zero number of treated individuals
                  in the absence of clinically infected individuals")
 })
@@ -884,13 +950,13 @@ test_that('calculate_treated() returns an empty Bitset when the parameter list c
     drug_time = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = 10)
-  
+
   treated <- calculate_treated(variables = variables,
                                clinical_infections = clinical_infections,
                                parameters = parameters,
                                timestep = timestep,
                                renderer = renderer)
-  
+
   expect_identical(object = treated$size(), expected = 0, info = "Error: calculate_treated() returning non-zero number of treated individuals
                  in the absence of clinical treatment or resistance parameters")
 })
@@ -918,7 +984,7 @@ test_that('Number of treatment failures matches number of individuals treated wh
                                             early_treatment_failure_probability = 1,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 10)
   parameters <- set_antimalarial_resistance(parameters = parameters,
                                             drug = 2,
@@ -929,25 +995,25 @@ test_that('Number of treatment failures matches number of individuals treated wh
                                             early_treatment_failure_probability = 1,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 20)
-  
+
   clinical_infections <- individual::Bitset$new(100)
   clinical_infections$insert(sample.int(n = 100, size = round(runif(n = 1, min = 10, max = 100)), replace = FALSE))
   timestep <- 5
   events <- create_events(parameters)
   variables <- create_variables(parameters = parameters)
   renderer <- individual::Render$new(timesteps = 10)
-  
+
   treated <- calculate_treated(variables = variables,
                                clinical_infections = clinical_infections,
                                parameters = parameters,
                                timestep = timestep,
                                renderer = renderer)
-  
+
   expect_identical(renderer$to_dataframe()[timestep,'n_early_treatment_failure'], renderer$to_dataframe()[timestep,'n_treated'] - renderer$to_dataframe()[timestep,'n_drug_efficacy_failures'], info = "Error: Number of
                  early treatment failures does not match number of treated individuals (minus drug efficacy failures) when artemisinin resistance proportion and
-                 and early treatment failure probability both equal 1") 
+                 and early treatment failure probability both equal 1")
 })
 
 test_that('calculate_treated() successfully adds additional resistance columns to the renderer', {
@@ -963,9 +1029,9 @@ test_that('calculate_treated() successfully adds additional resistance columns t
                                             early_treatment_failure_probability = 0.5,
                                             late_clinical_failure_probability = 0,
                                             late_parasitological_failure_probability = 0,
-                                            reinfection_during_prophylaxis_probability = 0, 
+                                            reinfection_during_prophylaxis_probability = 0,
                                             slow_parasite_clearance_time = 10)
-  
+
   clinical_infections <- individual::Bitset$new(20)$insert(1:20)
   timestep <- 5
   events <- create_events(parameters)
@@ -977,13 +1043,13 @@ test_that('calculate_treated() successfully adds additional resistance columns t
     dt = list(queue_update = mockery::mock())
   )
   renderer <- individual::Render$new(timesteps = 10)
-  
+
   treated <- calculate_treated(variables = variables,
                                clinical_infections = clinical_infections,
                                parameters = parameters,
                                timestep = timestep,
                                renderer = renderer)
-  
+
   calculate_treated_column_names <- c("ft",
                                       "n_treated",
                                       "n_drug_efficacy_failures",
diff --git a/tests/testthat/test-pev.R b/tests/testthat/test-pev.R
index 09c4cf7d..6e114b1b 100644
--- a/tests/testthat/test-pev.R
+++ b/tests/testthat/test-pev.R
@@ -36,7 +36,7 @@ test_that('Mass vaccination strategy parameterisation works', {
     ), "all(coverages >= 0) && all(coverages <= 1) is not TRUE",
     fixed = TRUE
   )
-  
+
   expect_error(
     parameters <- set_mass_pev(
       parameters,
@@ -150,6 +150,13 @@ test_that('Infection considers pev efficacy', {
     rep(.2, 4)
   )
 
+  infection_outcome <- CompetingOutcome$new(
+    targeted_process = function(timestep, target){
+      infection_process_resolved_hazard(timestep, target, variables, renderer, parameters)
+    },
+    size = parameters$human_population
+  )
+
   # remove randomness from infection sampling
   bernoulli_mock <- mockery::mock(c(1, 2))
   mockery::stub(calculate_infections, 'bernoulli_multi_p', bernoulli_mock)
@@ -164,16 +171,17 @@ test_that('Infection considers pev efficacy', {
     depth = 4
   )
 
-  calculate_infections(
+  infection_rates <- calculate_infections(
     variables = variables,
     bitten_humans = individual::Bitset$new(4)$insert(seq(4)),
     parameters = parameters,
     renderer = mock_render(timestep),
-    timestep = timestep
+    timestep = timestep,
+    infection_outcome
   )
 
   expect_equal(
-    mockery::mock_args(bernoulli_mock)[[1]][[1]],
+    rate_to_prob(infection_rates[infection_rates!=0]),
     c(0.590, 0.590, 0.215, 0.244),
     tolerance=1e-3
   )