Simplify competing hazards as a single random draw.

mrc-ide · Jul 10, 2024 · e2242b2 · e2242b2
1 parent 1f7310b
commit e2242b2
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 168 deletions.
diff --git a/R/competing_hazards.R b/R/competing_hazards.R
@@ -52,22 +52,20 @@ CompetingHazard <- R6::R6Class(
         'cbind',
         lapply(private$outcomes, function(x) x$rates)
       )
-      occur_rates <- rowSums(event_rates)
-      occur_rng <- private$rng(private$size)
-      occurs <- which(occur_rng < rate_to_prob(occur_rates))
 
-      norm_probs <- event_rates[occurs, , drop=FALSE] / occur_rates[occurs]
-      norm_probs[is.na(norm_probs)] <- 0
+      total_rates <- rowSums(event_rates)
+      probs <- rate_to_prob(total_rates) * (event_rates / total_rates)
+      probs[is.na(probs)] <- 0
 
-      cumulative <- rep(0, length(occurs))
-      event_rng <- private$rng(length(occurs))
+      rng <- private$rng(private$size)
 
+      cumulative <- rep(0, private$size)
       for (o in seq_along(private$outcomes)) {
-        next_cumulative <- cumulative + norm_probs[,o]
-        selected <- (event_rng > cumulative) & (event_rng <= next_cumulative)
+        next_cumulative <- cumulative + probs[,o]
+        selected <- which((rng > cumulative) & (rng <= next_cumulative))
         cumulative <- next_cumulative
 
-        target <- individual::Bitset$new(private$size)$insert(occurs[selected])
+        target <- individual::Bitset$new(private$size)$insert(selected)
         if (target$size() > 0){
           private$outcomes[[o]]$execute(t, target)
         }

diff --git a/tests/testthat/test-competing-hazards.R b/tests/testthat/test-competing-hazards.R
@@ -2,77 +2,32 @@
 # 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', {
-
+test_that("hazard resolves two disjoint 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
   )
-  
+
   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
-    )
+    rng = mockery::mock(c(.05, .3, .2, .5))
   )
 
-  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_1$set_rates(c(10, 0, 10, 0))
+  outcome_2$set_rates(c(0, 10, 0, 10))
 
   hazard$resolve(0)
-
+   
   mockery::expect_args(
-    outcome_1_process,
+   outcome_1_process,
     1,
     0,
     individual::Bitset$new(size)$insert(c(1, 3))
@@ -85,164 +40,76 @@ test_that('hazard resolves two unnormalised outcomes when all events occur', {
   )
 })
 
-test_that('hazard resolves two outcomes when some events occur', {
-
+test_that("hazard resolves two competing 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
   )
-  
+
   hazard <- CompetingHazard$new(
     outcomes = list(outcome_1, outcome_2),
-    rng = mockery::mock(
-      c(0, 0, 1, 1), # some events occur
-      c(.05, .3), # event_rng
-    )
+    rng = mockery::mock(c(.7, .3, .2, .6))
   )
 
-  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_1$set_rates(c(5, 5, 5, 5))
+  outcome_2$set_rates(c(5, 5, 5, 5))
 
   hazard$resolve(0)
-
+   
   mockery::expect_args(
-    outcome_1_process,
+   outcome_1_process,
     1,
     0,
-    individual::Bitset$new(size)$insert(1)
+    individual::Bitset$new(size)$insert(c(2, 3))
   )
   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
+    individual::Bitset$new(size)$insert(c(1, 4))
   )
 })
 
-test_that('hazard resolves when rates are set to one', {
-
+test_that("hazard resolves partial 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
   )
-
-  hazard <- CompetingHazard$new(
-    outcomes = list(outcome_1, outcome_2)
-  )
-
-  outcome_1$set_rates(rep(100, 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
-    )
+    outcomes = list(outcome_1, outcome_2),
+    rng = mockery::mock(c(.8, .4, .2, .6))
   )
 
-  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))
+  outcome_1$set_rates(prob_to_rate(rep(0.5, size)))
+  outcome_2$set_rates(prob_to_rate(rep(0.5, size)))
 
   hazard$resolve(0)
-
+   
   mockery::expect_args(
-    outcome_1_process,
+   outcome_1_process,
     1,
     0,
-    individual::Bitset$new(size)$insert(c(1, 3))
+    individual::Bitset$new(size)$insert(c(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)
+    individual::Bitset$new(size)$insert(c(2, 4))
   )
 })