test_that("Vector discounting equal to non-vector for single-length elements", { #Ongoing expect_equal(disc_ongoing_v(lcldr=0.035,lclprvtime=0.5, lclcurtime=3, lclval=2500), disc_ongoing(lcldr=0.035,lclprvtime=0.5, lclcurtime=3, lclval=2500)) #Instant expect_equal(disc_instant_v(lcldr=0.035, lclcurtime=3, lclval=2500), disc_instant(lcldr=0.035, lclcurtime=3, lclval=2500)) }) test_that("Discounting works with no discounting", { #Ongoing expect_equal(disc_ongoing_v(lcldr=0,lclprvtime=1, lclcurtime=2, lclval=2500), 2500) expect_equal(disc_ongoing(lcldr=0,lclprvtime=1, lclcurtime=2, lclval=2500), 2500) #Instant expect_equal(disc_instant_v(lcldr=0, lclcurtime=2, lclval=2500), 2500) expect_equal(disc_instant(lcldr=0, lclcurtime=2, lclval=2500), 2500) #Cycle expect_equal(disc_cycle_v(lcldr=0,lclprvtime=1, lclcurtime=2, lclval=2500, cyclelength=1/12, starttime=0), 12*2500) expect_equal(disc_cycle(lcldr=0,lclprvtime=1, lclcurtime=2, lclval=2500, cyclelength=1/12, starttime=0), 12*2500) }) test_that("Discounting works with odd numbers", { #Ongoing expect_equal(disc_ongoing_v(lcldr=0,lclprvtime=0, lclcurtime=Inf, lclval=2500), Inf) expect_equal(disc_ongoing(lcldr=0,lclprvtime=0, lclcurtime=Inf, lclval=2500), Inf) expect_equal(disc_ongoing_v(lcldr=0.035,lclprvtime=0, lclcurtime=0, lclval=2500), 0) expect_equal(disc_ongoing(lcldr=0.035,lclprvtime=0, lclcurtime=0, lclval=2500), 0) expect_equal(disc_ongoing_v(lcldr=0.035,lclprvtime=5, lclcurtime=5, lclval=2500), 0) expect_equal(disc_ongoing(lcldr=0.035,lclprvtime=5, lclcurtime=5, lclval=2500), 0) #Inf*0 gives NaN, while the element-wise function just check wehterh prevtime and curtime are equal expect_equal(disc_ongoing_v(lcldr=0.035,lclprvtime=5, lclcurtime=5, lclval=Inf), NaN) expect_equal(disc_ongoing(lcldr=0.035,lclprvtime=5, lclcurtime=5, lclval=Inf), NaN) #Instant expect_equal(disc_instant_v(lcldr=0, lclcurtime=Inf, lclval=2500), 2500) expect_equal(disc_instant(lcldr=0, lclcurtime=Inf, lclval=2500), 2500) expect_equal(disc_instant_v(lcldr=0.035, lclcurtime=0, lclval=2500), 2500) expect_equal(disc_instant(lcldr=0.035, lclcurtime=0, lclval=2500), 2500) expect_equal(disc_instant_v(lcldr=0.035, lclcurtime=5, lclval=2500), 2104.93292) expect_equal(disc_instant(lcldr=0.035, lclcurtime=5, lclval=2500), 2104.93292) #Inf*0 gives NaN expect_equal(disc_instant_v(lcldr=0.035, lclcurtime=5, lclval=Inf), Inf) expect_equal(disc_instant(lcldr=0.035, lclcurtime=5, lclval=Inf), Inf) expect_equal(disc_instant_v(lcldr=5, lclcurtime=5, lclval=2500), 0.32150206) expect_equal(disc_instant(lcldr=5, lclcurtime=5, lclval=2500), 0.32150206) #Cycle expect_equal(disc_cycle_v(lcldr=0,lclprvtime=1, lclcurtime=2, lclval=2500, cyclelength=1/12, starttime=1), 12*2500) expect_equal(disc_cycle_v(lcldr=0,lclprvtime=0, lclcurtime=0, lclval=2500, cyclelength=1/12, starttime=0), 2500) expect_equal(disc_cycle_v(lcldr=0.035,lclprvtime=0, lclcurtime=0, lclval=2500, cyclelength=2, starttime=0), 2500) #Inf*0 gives NaN expect_equal(disc_cycle_v(lcldr=0.035,lclprvtime=5, lclcurtime=5, lclval=Inf, cyclelength=1/12, starttime=0), NaN) }) test_that("Vectorial discounting working as expected with vectors", { #Ongoing expect_equal(disc_ongoing_v(lcldr=0.035,lclprvtime=c(0.5,0.5,0.5), lclcurtime=c(3,3,3), lclval=c(0,1000,Inf)), c(0,2354.66015,Inf)) #Instant expect_equal(disc_instant_v(lcldr=0.035, lclcurtime=c(3,3,3), lclval=c(0,1000,Inf)), c(0,901.9427,Inf)) #Cycle expect_equal(disc_cycle_v(lcldr=0.035, lclcurtime=c(3,3,3), lclval=c(0,1000,Inf),lclprvtime=c(0.5,0.5,0.5), cyclelength=c(1/12,1/12,1/12),starttime=c(0,0,0)), c(0,28296.443022132232727,Inf)) expect_equal(disc_cycle_v(0.035,0,1,1,1000,0,5),1000) expect_equal(disc_cycle_v(0.035,1,1,1,1000,0,5),966.1835748792273) expect_equal(disc_cycle_v(0.035,0,1,0,1000,0,5),1000) expect_equal(disc_cycle_v(0.035,0,1,2,1000,0,5),1966.1835748792273) }) test_that("Create indicators works correctly",{ expect_equal( create_indicators( 2, 10, c(1,1) ), c(0,1) ) expect_equal( create_indicators( 2, 10, c(1,1), 5 ), c(0,0) ) expect_equal( create_indicators( 6, 10, c(1,1), 5 ), c(1,0) ) expect_equal( create_indicators( 9, 10, c(1,1), 5 ), c(0,0) ) expect_error( create_indicators( 12, 10, c(1,1), 5 ) ) expect_error( create_indicators( 9, 10, rep(2,20), 5 ) ) expect_error( create_indicators( 9, 10, rep(2,3), 20 ) ) }) test_that("Pick values vectorized work correctly",{ expect_equal( pick_val_v(base = list(2,3,c(2, 3, 4)), psa =sapply(1:3, function(x) eval(call( c("rnorm","rnorm","rdirichlet")[[x]], 1, c(2,3,list(c(2, 3, 4)))[[x]], c(0.1,0.1,NULL)[[x]] ))), sens = list(4,5,c(0.4,0.8,0.1)), psa_ind = FALSE, sens_ind = TRUE, indicator=list(1,2,c(3,4,5)), names_out=c("util","util2","dirichlet_vector") , indicator_sens_binary = FALSE, sens_iterator = 5, distributions = list("rnorm","rnorm","rdirichlet"), covariances = list(0.1,0.1,NULL) ), list(util = 2, util2 = 3, dirichlet_vector = c(0.36, 0.54, 0.1)) ) expect_equal( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = FALSE, sens_ind = FALSE, indicator=c(1,0) ), list(0,0) ) expect_equal( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = FALSE, sens_ind = TRUE, indicator=c(1,0) ), list(2,0) ) expect_equal( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = FALSE, sens_ind = TRUE, indicator=c(0,1) ), list(0,3) ) expect_equal( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = FALSE, sens_ind = TRUE, indicator=c(1,1) ), list(2,3) ) expect_equal( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = TRUE, sens_ind = TRUE, indicator=c(1,1) ), list(2,3) ) expect_error( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = TRUE, sens_ind = TRUE, indicator=c(1,5) ) ) expect_error( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = 5, sens_ind = TRUE, indicator=c(1,1) ) ) expect_error( pick_val_v( base = c(0,0), psa =c(rnorm(1,0,0.1),rnorm(1,0,0.1)), sens = c(2,3), psa_ind = TRUE, sens_ind = 3, indicator=c(1,1) ) ) expect_equal( pick_val_v( base = c(0,0), psa ={c(draw_tte(1,'norm',0,0.1,seed=1),draw_tte(1,'norm',0,0.1,seed=2))}, sens = c(2,3), psa_ind = TRUE, sens_ind = TRUE, indicator=c(0,0) ), {list(draw_tte(1,'norm',0,0.1,seed=1),draw_tte(1,'norm',0,0.1,seed=2))} ) expect_equal( pick_val_v( base = c(0,0), psa ={c(draw_tte(1,'norm',0,0.1,seed=1),draw_tte(1,'norm',0,0.1,seed=2))}, sens = c(2,3), psa_ind = TRUE, sens_ind = TRUE, indicator=c(1,0) ), {list(2,draw_tte(1,'norm',0,0.1,seed=2))} ) }) test_that("Conditional Multivariate normal works as expected",{ expect_equal( cond_mvn(mu = c(1, 2, 3), Sigma = matrix(c(0.2, 0.05, 0.1, 0.05, 0.3, 0.05, 0.1, 0.05, 0.4), nrow = 3), i = 1:2, xi = c(1.2,2.3), full_output = TRUE )$mean, c(1.2, 2.3, 3.1217391) ) expect_equal( cond_mvn(mu = c(1, 2, 3), Sigma = matrix(c(0.2, 0.05, 0.1, 0.05, 0.3, 0.05, 0.1, 0.05, 0.4), nrow = 3), i = 1:2, xi = c(1.2,2.3), full_output = TRUE )$covariance, structure(c(0, 0, 0, 0, 0, 0, 0, 0, 0.347826086956522), dim = c(3L, 3L)) ) expect_error( cond_mvn(mu = c(1, 2, 3), Sigma = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 3), i = 1:2, xi = c(1.2,2.3), full_output = TRUE ) ) expect_error( cond_mvn(mu = c(1, 2, 3), Sigma = matrix(c(0, 0, 0, 0, 0, 0, 0, 0, 0), nrow = 3), i = 5, xi = c(1.2,2.3), full_output = TRUE ) ) }) test_that("Model Reactions Interactivity summary can be created",{ expr <- substitute({ a <- sum(5+7) ggplot() data.frame(x=1,b=2) list(b=5) a <- list(s=7) j <- 6 if(TRUE){modify_event(list(j=5))} l <- 9 modify_item(list(afsa=ifelse(TRUE,"asda",NULL))) modify_item_seq(list( o_exn = o_exn + 1, a = NULL, b = if(a){"CZ"}else{"AW"}, rnd_prob_exn_sev = runif(1), exn_sev = rnd_prob_exn_sev <= p_sev, o_exn_mod = o_exn_mod + if(exn_sev) { 0 } else { 1 }, o_exn_sev = o_exn_sev + if(exn_sev) { 1 } else { 0 }, o_rec_time_without_exn = (o_exn == 0) * 1, o_rec_time_without_exn_sev = (o_exn_sev == 0) * 1, o_c_exn = if(exn_sev) { c_sev } else { c_mod }, o_other_c_exn_mod = if(exn_sev) { 0 } else { c_mod }, o_other_c_exn_sev = if(exn_sev) { c_sev } else { 0 }, o_qloss_exn = -if(exn_sev) { q_sev } else { q_mod }, o_other_qloss_exn_mod = -if(exn_sev) { 0 } else { q_mod }, o_other_qloss_exn_sev = -if(exn_sev) { q_sev } else { 0 }, o_qloss_cg_exn = -if(exn_sev) { q_cg_sev } else { q_cg_mod }, o_other_qloss_cg_exn_mod = -if(exn_sev) { 0 } else { q_cg_mod }, o_other_qloss_cg_exn_sev = -if(exn_sev) { q_cg_sev } else { 0 }, o_q = utility, o_other_q_gold1 = if(gold == 1) { utility } else { 0 }, o_other_q_gold2 = if(gold == 2) { utility } else { 0 }, o_other_q_gold3 = if(gold == 3) { utility } else { 0 }, o_other_q_gold4 = if(gold == 4) { utility } else { 0 }, o_other_q_on_dup = if(on_dup) { utility } else { 0 }, n_exn = n_exn + 1, n_exn_mod = n_exn_mod + (1 - exn_sev), n_exn_sev = n_exn_sev + exn_sev, u_adj_exn_lt = u_adj_exn_lt + if(exn_sev) { u_adj_sev_lt } else { u_adj_mod_lt }, utility = u_gold - u_adj_exn_lt - u_mace_lt, o_rec_utility = utility, rnd_exn = runif(1) )) if(a==1){ modify_item(list(a=list(6+b))) modify_event(list(e_exn = curtime + 14 / days_in_year + qexp(rnd_exn, r_exn))) } else{ modify_event(list(e_exn = curtime + 14 / days_in_year + qexp(rnd_exn, r_exn))) if(a>6){ modify_item(list(a=8)) } } if (sel_resp_incl == 1 & on_dup == 1) { modify_event(list(e_response = curtime, z = 6)) } }) expect_length(ast_as_list(expr),13) expect_type(ast_as_list(expr),"list") expect_equal(class(extract_elements_from_list(ast_as_list(expr))),"data.frame") expect_length(extract_elements_from_list(ast_as_list(expr)),4) #4 columns expect_equal(nrow(extract_elements_from_list(ast_as_list(expr))),43) #43 items/events changed, including assignments a <- add_reactevt(name_evt="example", input={ a <- 5 modify_item(list(w=5)) }) expect_equal(nrow(extract_from_reactions(a)),2) #2 items/events changed, including assignments expect_equal(extract_from_reactions(a), data.table(event = c("example","example"), name = c("a","w"), type = c("item","item"), conditional_flag = c(FALSE,FALSE), definition = c("5","5")) ) }) test_that("add_tte works as expected", { initial_data <- list() arm <- "control" evts <- c("start", "end") result <- add_tte(.data = initial_data, arm = arm, evts = evts, input = { start <- 0 end <- 100 }) expect_true("control" %in% names(result)) expect_equal(result$control$evts, evts) }) test_that("modify_item modifies input items correctly", { input_list_arm <- environment() input_list_arm$qaly_default_instant = 100 input_list_arm$debug = FALSE input_list_arm$accum_backwards = FALSE modify_item(list("qaly_default_instant" = 200)) expect_equal(input_list_arm$qaly_default_instant, 200) modify_item(list(new_cost = 300)) expect_equal(input_list_arm$new_cost, 300) }) test_that("modify_event modifies events correctly", { input_list_arm <- environment() input_list_arm$qaly_default_instant = 100 input_list_arm$debug = FALSE input_list_arm$cur_evtlist = FALSE input_list_arm$accum_backwards = FALSE input_list_arm$cur_evtlist = c(ae = 5, nat.death = 100) # Modify an existing event modify_event(list(ae = 10)) expect_equal(input_list_arm$cur_evtlist[["ae"]], 10) # Create new event if not exists modify_event(list(new_event = 50), create_if_null = TRUE) expect_equal(input_list_arm$cur_evtlist[["new_event"]], 50) # Ignore non-existent event expect_warning(modify_event(list(nonexistent = 20), create_if_null = FALSE)) expect_error(input_list_arm$cur_evtlist[["nonexistent"]]) }) test_that("new_event adds new events correctly", { input_list_arm <- environment() input_list_arm$qaly_default_instant = 100 input_list_arm$debug = FALSE input_list_arm$cur_evtlist = FALSE input_list_arm$accum_backwards = FALSE input_list_arm$cur_evtlist = c() new_event(list("ae" = 5)) expect_equal(input_list_arm$cur_evtlist[["ae"]], 5) expect_error(new_event(list("not_numeric" = "five")), "New event times are not all numeric, please review") }) test_that("replicate_profiles works correctly", { profiles <- data.frame(id = 1:10, age = rnorm(10, 60, 5)) # Test replication with replacement set.seed(42) result <- replicate_profiles(profiles, replications = 20, replacement = TRUE) expect_equal(nrow(result), 20) expect_true(all(result$id %in% profiles$id)) # Test replication without replacement set.seed(42) result_no_replacement <- replicate_profiles(profiles, replications = 10, replacement = FALSE) expect_equal(nrow(result_no_replacement), 10) expect_equal(sort(result_no_replacement$id), sort(profiles$id)) }) test_that("modify_item_seq works sequentially", { input_list_arm <- environment() input_list_arm$a <- 1 input_list_arm$b <- 2 input_list_arm$curtime <- 1 input_list_arm$accum_backwards <- FALSE input_list_arm$debug <- FALSE # Test sequential modification modify_item_seq(list(a = 3, b = a + 2)) expect_equal(input_list_arm$a, 3) expect_equal(input_list_arm$b, 5) # Test debug mode input_list_arm$debug <- TRUE input_list_arm$log_list <- list() modify_item_seq(list(a = 4, c = b * 2)) expect_equal(input_list_arm$a, 4) expect_equal(input_list_arm$c, 10) expect_true(length(input_list_arm$log_list) > 0) }) test_that("add_reactevt adds reactions correctly", { # Create an empty data list data_list <- list() # Add a reaction result <- add_reactevt(.data = data_list, name_evt = "start", input = { curtime <- Inf }) expect_true("start" %in% names(result)) # Test error handling for invalid event name expect_error(add_reactevt(name_evt = c("evt1", "evt2"), input = {}), "name_evt argument in add_reactevt should be a single string with at least 2 characters") }) # Luck_adj ---------------------------------------------------------------- test_that("luck_adj adjusts luck correctly", { # Test single values adj <- luck_adj(prevsurv = 0.8, cursurv = 0.6, luck = 0.9, condq = TRUE) expect_true(adj > 0 & adj < 1) # Test vectorized adjustment adj_vec <- luck_adj(prevsurv = c(0.8, 0), cursurv = c(0.6, 0.5), luck = c(0.9, 0.8), condq = TRUE) expect_equal(length(adj_vec), 2) expect_equal(adj_vec[2], 0.8) # Test conditional adjustment adj_cond <- luck_adj(prevsurv = 0.8, cursurv = 0.6, luck = 0.9, condq = FALSE) expect_true(adj_cond > 0 & adj_cond < 1) }) # Random Streams ---------------------------------------------------------- test_that("random_stream initializes correctly", { random_stream_instance <- random_stream(10) expect_type(random_stream_instance, "environment") expect_equal(length(random_stream_instance$stream), 10, info = "The stream should initialize with 10 elements") }) test_that("draw_n draws correct number of elements", { random_stream_instance <- random_stream(10) drawn_numbers <- random_stream_instance$draw_n(3) expect_equal(length(drawn_numbers), 3, info = "draw_n should draw 3 elements") expect_equal(length(random_stream_instance$stream), 7, info = "The stream should have 7 elements left after drawing 3") }) test_that("draw_n with larger n regenerates stream", { random_stream_instance <- random_stream(5) expect_warning(drawn_numbers_large <- random_stream_instance$draw_n(8), "Stream is smaller than the number of numbers drawn", info = "Should warn when trying to draw more elements than available") expect_equal(length(drawn_numbers_large), 8, info = "After regenerating, draw_n should return 8 elements") expect_equal(length(random_stream_instance$stream), 0, info = "After drawing all elements, the stream should be empty") }) test_that("generate_stream changes stream length", { random_stream_instance <- random_stream(5) random_stream_instance$generate_stream(20) expect_equal(length(random_stream_instance$stream), 20, info = "generate_stream should change stream to the correct size") }) # Time covariate tte ------------------------------------------------------ library(flexsurv) test_that("qtimecov returns numeric scalar within bounds", { param_fun_factory <- function(p0, p1, p2, p3) { function(t) p0 + p1 * t + p2 * t^2 + p3 * (floor(t) + 1) } rate_exp <- param_fun_factory(0.1, 0, 0, 0) set.seed(1) tte <- qtimecov(runif(1), a_fun = rate_exp, dist = "exp") expect_type(tte, "double") expect_length(tte, 1) expect_gt(tte, 0) expect_lt(tte, 100) }) test_that("qtimecov works for all supported distributions", { set.seed(1) param_fun_factory <- function(p0, p1, p2, p3) { function(.time) p0 + p1 * .time + p2 * .time^2 + p3 * (floor(.time) + 1) } # 1. Exponential rate_exp <- param_fun_factory(0.1, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = rate_exp, dist = "exp")) expect_equal(qtimecov(luck = 0.5,a_fun = rate_exp,dist = "exp", dt = 0.001 ),qexp(0.5,0.1), tolerance = 0.01) # 2. Gamma shape <- param_fun_factory(2, 0, 0, 0) rate <- param_fun_factory(0.2, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = shape, b_fun = rate, dist = "gamma")) expect_equal(qtimecov(luck = 0.5,a_fun = shape, b_fun = rate, dist = "gamma", dt = 0.001 ),qgamma(0.5,2,0.2), tolerance = 0.01) # 3. Lognormal meanlog <- param_fun_factory(log(10) - 0.5^2 / 2, 0, 0, 0) sdlog <- param_fun_factory(0.5, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = meanlog, b_fun = sdlog, dist = "lnorm")) expect_equal(qtimecov(0.5, a_fun = meanlog, b_fun = sdlog, dist = "lnorm",dt=0.01), qlnorm(0.5,log(10) - 0.5^2 / 2,0.5), tolerance = 0.01) # 4. Normal mean <- param_fun_factory(10, 0, 0, 0) sd <- param_fun_factory(2, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = mean, b_fun = sd, dist = "norm")) expect_equal(qtimecov(0.5, a_fun = mean, b_fun = sd, dist = "norm",dt=0.01), qnorm(0.5,10,2), tolerance = 0.01) # 5. Weibull shape <- param_fun_factory(2, 0, 0, 0) scale <- param_fun_factory(10, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = shape, b_fun = scale, dist = "weibull")) expect_equal(qtimecov(0.5, a_fun = shape, b_fun = scale, dist = "weibull",dt=0.01), qweibull(0.5,2,10), tolerance = 0.01) # 5.1 WeibullPH shape <- param_fun_factory(2, 0, 0, 0) scale <- param_fun_factory(0.01, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = shape, b_fun = scale, dist = "weibullPH")) expect_equal(qtimecov(0.5, a_fun = shape, b_fun = scale, dist = "weibullPH",dt=0.01), flexsurv::qweibullPH(0.5,2,0.01), tolerance = 0.01) # 6. Loglogistic shape <- param_fun_factory(2.5, 0, 0, 0) scale <- param_fun_factory(7.6, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = shape, b_fun = scale, dist = "llogis")) expect_equal(qtimecov(0.5, a_fun = shape, b_fun = scale, dist = "llogis",dt=0.01), flexsurv::qllogis(0.5,2.5,7.6), tolerance = 0.01) # 7. Gompertz shape <- param_fun_factory(0.01, 0, 0, 0) rate <- param_fun_factory(0.091, 0, 0, 0) expect_silent(qtimecov(runif(1), a_fun = shape, b_fun = rate, dist = "gompertz")) expect_equal(qtimecov(0.5, a_fun = shape, b_fun = rate, dist = "gompertz",dt=0.01), qgompertz(0.5,0.01,0.091), tolerance = 0.01) rate_exp <- function(t) 0.1 init_luck <- 0.95 expect_equal(qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.001),{ a <- qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.001, return_luck = TRUE,max_time = 10) qtimecov(luck = a$luck,a_fun = rate_exp,dist = "exp", dt = 0.001, start_time=a$tte)}, tolerance = 0.01) init_luck <- 0.99 expect_equal(qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.001),{ a <- qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.001, return_luck = TRUE,max_time = 5) qtimecov(luck = a$luck,a_fun = rate_exp,dist = "exp", dt = 0.001, start_time=a$tte)}, tolerance = 0.01) init_luck <- 0.3 expect_equal(qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.001),{ a <- qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.001, return_luck = TRUE,max_time = 1) qtimecov(luck = a$luck,a_fun = rate_exp,dist = "exp", dt = 0.001, start_time=a$tte)}, tolerance = 0.01) init_luck <- 0.3 expect_equal(qtimecov(luck = init_luck,a_fun = rate_exp,b_fun = function(t) 90000 ,dist = "weibull", dt = 0.001),{ a <- qtimecov(luck = init_luck,a_fun = rate_exp,b_fun = function(t) 90000,dist = "weibull", dt = 0.001, return_luck = TRUE,max_time = 1) qtimecov(luck = a$luck,a_fun = rate_exp,b_fun = function(t) 90000,dist = "weibull", dt = 0.001, start_time=a$tte)}, tolerance = 0.01) rate_exp <- function(.time) 0.1 rate_exp2 <- function(.time) 0.2 time_change <- 10 init_luck <- 0.95 expect_equal({ a <- qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.005, max_time = time_change, return_luck = TRUE) qtimecov(luck = a$luck,a_fun = rate_exp2,dist = "exp", dt = 0.005, start_time=a$tte) },{ new_luck <- luck_adj(prevsurv = 1 - pexp(q=time_change,rate_exp(1)), cursurv = 1 - pexp(q=time_change,rate_exp2(1)), luck = init_luck, condq = FALSE) #time 10 change qexp(new_luck,rate_exp2(1)) }, tolerance = 0.01) # time varying and an event ----------------------------------------------- rate_exp <- function(.time) 0.1 + 0.01*.time * 0.00001*.time^2 rate_exp2 <- function(.time) 0.2 + 0.02*.time time_change <- 8 init_luck <- 0.95 expect_equal({ a <- qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 0.005, max_time = time_change, return_luck = TRUE) qtimecov(luck = a$luck,a_fun = rate_exp2,dist = "exp", dt = 0.005, start_time=a$tte) },{# Manually reproduce what qtimecov does from a$tte t <- 0 luck <- init_luck dt <- 0.005 repeat { t <- t + dt s_prev <- 1 - pexp(t - dt, rate = rate_exp(t - dt)) s_curr <- 1 - pexp(t, rate = rate_exp(t - dt)) luck <- luck_adj(prevsurv = s_prev, cursurv = s_curr, luck = luck, condq = TRUE) res_tte <- qcond_exp(luck, rate = rate_exp(t)) total_tte <- t + res_tte if (res_tte <= dt || total_tte <= t || t >= time_change) { break } } if (total_tte <= time_change) { return(total_tte) } # Phase 2: after change repeat { t <- t + dt s_prev <- 1 - pexp(t - dt, rate = rate_exp2(t - dt)) s_curr <- 1 - pexp(t, rate = rate_exp2(t - dt)) luck <- luck_adj(prevsurv = s_prev, cursurv = s_curr, luck = luck, condq = TRUE) res_tte <- qcond_exp(luck, rate = rate_exp2(t)) total_tte <- t + res_tte if (res_tte <= dt || total_tte <= t) { break } } total_tte }, tolerance = 0.01) rate_exp <- function(.time) 0.1 + floor(.time)*0.01 init_luck <- 0.95 #30x faster expect_equal({ a <- qtimecov(luck = init_luck,a_fun = rate_exp,dist = "exp", dt = 1) a },{ time_vec <- 0:100 new_luck <- init_luck for (i in 2:length(time_vec)) { new_luck <- luck_adj(prevsurv = 1 - pexp(q=time_vec[i-1],rate_exp(time_vec[i-1])), cursurv = 1 - pexp(q=time_vec[i],rate_exp(time_vec[i])), luck = new_luck, condq = FALSE) qexp(new_luck,rate_exp(1)) } a }, tolerance = 0.01) expect_equal({ bs_rate <- 0.08 time_cov <- 0.0002271 bs_rate_f <- function(.time) bs_rate + time_cov * floor(.time) luck <- 0.403 t <- 0 # First rate and quantile (from time 0) new_rate <- bs_rate_f(0) ttdeath <- qexp(luck, rate = new_rate) # Unconditional quantile loop while (t < ttdeath) { old_rate <- new_rate new_rate <- bs_rate_f(t) # Both survival values are evaluated at time t (NOT t+1) prev_surv <- 1 - pexp(q = t, rate = old_rate) cur_surv <- 1 - pexp(q = t, rate = new_rate) # Use condq = FALSE for unconditional approach luck <- luck_adj(prevsurv = prev_surv, cursurv = cur_surv, luck = luck, condq = FALSE) # Unconditional quantile from time 0 ttdeath <- qexp(luck, rate = new_rate) t <- t + 1 } # Final time-to-event ttdeath },{ new_rate <-bs_rate <- 0.08 bs_rate_f <- function(.time){bs_rate + (time_cov * floor(.time))} luck <- 0.403 time_cov <- 0.0002271 t <- 0 qtimecov(luck = luck,a_fun = bs_rate_f,dist = "exp", dt = 1) }, tolerance = 0.01) }) test_that("qtimecov throws error for unsupported distribution", { dummy <- function(.time) 1 expect_error(qtimecov(runif(1), a_fun = dummy, dist = "beta"), "Unsupported distribution") }) test_that("qtimecov respects max_time bound", { slow_fun <- function(.time) 0.00001 tte <- qtimecov( luck = 0.999, a_fun = slow_fun, dist = "exp", max_time = 10 ) expect_lte(tte, 10) }) test_that("adj_val works as intended",{ bs_age <- 1 vec <- 1 - seq(from =0, to = 0.2, by = 0.02) expect_equal(adj_val(0,0,by=1, vec[floor(.time + bs_age)], discount = Inf), 0) expect_equal(adj_val(0,5,by=1, vec[floor(.time + bs_age)], discount = Inf), 0) expect_error(adj_val(0,20,by=1, vec[floor(.time + bs_age)])) expect_error(adj_val(0,20,by=1, vec[floor(.time + bs_age)])) expect_error(adj_val(0,-5,by=1, vec[floor(.time + bs_age)])) expect_equal(adj_val(0,0,by=1, vec[floor(.time + bs_age)]), 0) expect_equal(adj_val(0,0.1,by=1, vec[floor(.time + bs_age)], discount = 0),1) expect_equal(adj_val(0,1.1,by=1, vec[floor(.time + bs_age)] * .time, discount = 0),(1*0+0.98*0.1)/1.1) expect_equal(adj_val(0,0.1,by=1, vec[floor(.time + bs_age)] * .time, discount = 0),0) expect_equal(adj_val(0,0.1,by=1, vec[floor(.time + bs_age)] * .time, discount = 0.03),0) expect_equal(adj_val(8,9,by=0.2, vec[floor(.time + bs_age)], discount = 0),0.84) expect_equal(adj_val(8,9,by=0.5, vec[floor(.time + bs_age)], discount = Inf),0) expect_equal(adj_val(8,9,by=0.5, 1),1) expect_equal(adj_val(0,4,by=1, .time),1.5) expect_equal({ val <- 0.8 * adj_val(0,5.2,by=1, vec[floor(.time + bs_age)], discount = 0.03) disc_ongoing_v(0.03, 0, 5.2, val) },{ a <- 0 time_vec <- c(0,1,2,3,4,5,5.2) for (i in 1:(length(time_vec)-1)) { a <- a + disc_ongoing_v(0.03, time_vec[i], time_vec[i+1], 0.8 * vec[floor(time_vec[i] + bs_age)]) } a }) #17x faster expect_equal({ val <- 0.8 * adj_val(0,5.2,by=1, vec[floor(.time + bs_age)], discount = 0) disc_ongoing_v(0, 0, 5.2, val) },{ a <- 0 time_vec <- c(0,1,2,3,4,5,5.2) for (i in 1:(length(time_vec)-1)) { a <- a + disc_ongoing_v(0, time_vec[i], time_vec[i+1], 0.8 * vec[floor(time_vec[i] + bs_age)]) } a }) }) test_that("qcondweibull and qcondweibullPH can give same results", { # Set base R parameters shape <- 2 scale_base <- 10 # Convert to PH parameter scale_PH <- (1 / scale_base)^shape # = 0.01 # Set common inputs p <- 0.5 t0 <- 5 # Evaluate both expect_equal(qcond_weibull(p, shape, scale_base, lower_bound = t0), qcond_weibullPH(p, shape, scale_PH, lower_bound = t0)) } )