From e2242b26ae26cbca33486392a00cda3df5be30bb Mon Sep 17 00:00:00 2001 From: Paul Lietar Date: Wed, 10 Jul 2024 13:44:59 +0100 Subject: [PATCH] Simplify competing hazards as a single random draw. --- R/competing_hazards.R | 18 ++- tests/testthat/test-competing-hazards.R | 183 ++++-------------------- 2 files changed, 33 insertions(+), 168 deletions(-) diff --git a/R/competing_hazards.R b/R/competing_hazards.R index 8acfc387..b40ff175 100644 --- 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 index 12203679..e1240d27 100644 --- 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)) ) })