Refactor `expected_event()` for better computational efficiency

[LittleBeannie]

x <- microbenchmark::microbenchmark(
  t1 = gsDesign::eEvents(lambda = c(0.05, 0.02, 0.01),  # event hazard rates
                         eta = 0.01,                    # dropout hazard rates
                         gamma = c(5, 10, 20),          # rates of entry
                         R =  c(2, 1, 2),               # duration of time periods for recruitment rates
                         S =  c(1, 1),                  # duration of piecewise constant event rates
                         T = 12,                        # time of analysis
                         Tfinal = NULL,                 # Study duration
                         minfup = 12 - 5,
                         digits = 4),
  t2 = expected_event(enroll_rate = define_enroll_rate(duration = c(2, 1, 2), rate = c(5, 10, 20)),
                      fail_rate = define_fail_rate(duration = c(1, 1, 100), fail_rate = c(0.05, 0.02, 0.01), dropout_rate = 0.01),
                      total_duration = 12,
                      simple = FALSE),
  times = 1000
)
levels(x$expr) <- c("eEvents", "expected_event")

x

[LittleBeannie]

option1: gsDesign::eEvent, will need more comments in gsDesign::eEvents and gsDesign:::eEvents1
option2: tibble -> base R / data.table (preference)

[LittleBeannie]

profile expected_event and eEevent first

[LittleBeannie]

related issue in gsDesign: keaven/gsDesign#91

[LittleBeannie]

The output of the current gsDesign:::eEvents1 can not be transformed to gsDesign::expected_event. More technical work is needed.

t1 = gsDesign:::eEvents1(lambda = c(0.05, 0.02, 0.01),  # event hazard rates
                         eta = 0.01,                    # dropout hazard rates
                         gamma = c(5, 10, 20),          # rates of entry
                         R =  c(2, 1, 2),               # duration of time periods for recruitment rates
                         S =  c(1, 1, 100),                  # duration of piecewise constant event rates
                         T = 20,                        # time of analysis
                         Tfinal = 20,                 # Study duration
                         minfup = NULL,
                         simple = FALSE)

t2 = expected_event(enroll_rate = define_enroll_rate(duration = c(2, 1, 2), rate = c(5, 10, 20)),
                    fail_rate = define_fail_rate(duration = c(1, 1, 100), fail_rate = c(0.05, 0.02, 0.01), dropout_rate = 0.01),
                    total_duration = 20,
                    simple = FALSE)
# output of eEvents1
tibble(i = 1:length(t1$A),
       t = t1$T3,
       gamma = t1$gam,
       lambda = t1$lam,
       eta = t1$et,
       q = t1$q,
       Q = t1$Q,
       p = t1$p,
       P = t1$P,
       B = t1$B,
       A = t1$A) 

# output of expected_event
t2

[jdblischak]

I've started benchmarking/profiling expected_event(). One surprising thing I noticed is that define_enroll_rate() and define_fail_rate() alone both take longer than gsDesign::eEvents():

library("gsDesign2")

enroll_rate_precomputed <- define_enroll_rate(duration = c(2, 1, 2), rate = c(5, 10, 20))
fail_rate_precomputed <- define_fail_rate(duration = c(1, 1, 100), fail_rate = c(0.05, 0.02, 0.01), dropout_rate = 0.01)

x <- microbenchmark::microbenchmark(
  t1 = gsDesign::eEvents(lambda = c(0.05, 0.02, 0.01),  # event hazard rates
                         eta = 0.01,                    # dropout hazard rates
                         gamma = c(5, 10, 20),          # rates of entry
                         R =  c(2, 1, 2),               # duration of time periods for recruitment rates
                         S =  c(1, 1),                  # duration of piecewise constant event rates
                         T = 12,                        # time of analysis
                         Tfinal = NULL,                 # Study duration
                         minfup = 12 - 5,
                         digits = 4),
  t2 = expected_event(enroll_rate = define_enroll_rate(duration = c(2, 1, 2), rate = c(5, 10, 20)),
                      fail_rate = define_fail_rate(duration = c(1, 1, 100), fail_rate = c(0.05, 0.02, 0.01), dropout_rate = 0.01),
                      total_duration = 12,
                      simple = FALSE),
  t3 = expected_event(enroll_rate = enroll_rate_precomputed,
                      fail_rate = fail_rate_precomputed,
                      total_duration = 12,
                      simple = FALSE),
  t4 = define_enroll_rate(duration = c(2, 1, 2), rate = c(5, 10, 20)),
  t5 = define_fail_rate(duration = c(1, 1, 100), fail_rate = c(0.05, 0.02, 0.01), dropout_rate = 0.01),
  times = 1000
)
levels(x$expr) <- c("eEvents", "expected_event", "expected_event_precomputed",
                    "define_enroll_rate", "define_fail_rate")
x
## Unit: microseconds
##                        expr     min       lq       mean   median       uq      max neval
##                     eEvents   166.5   234.00   282.3544   265.00   291.55   4869.4  1000
##              expected_event 24264.8 25785.00 30517.7125 28493.65 31403.75 140115.2  1000
##  expected_event_precomputed 22652.4 23914.90 27908.5403 26074.30 29175.10 100711.9  1000
##          define_enroll_rate   674.3   740.90   939.1822   797.65   899.30   9278.9  1000
##            define_fail_rate   950.5  1039.45  1233.2822  1096.15  1225.50   8289.8  1000

[LittleBeannie]

define_enroll_rate   674.3   740.90   939.1822   797.65   899.30   9278.9  1000
##            define_fail

Interesting but this is something we may expect... considering define_enroll_rate and define_fail_rate use tibble operations and there are some if statements to check inputs, while eEvents is pure base-R.

Thanks for sharing the finding and this explains our motivation to transfer expected_event from tibble to base-R.

[nanxstats]

Maybe consider replacing the tibble() calls in the define_*() functions with the low-level constructor new_tibble() calls which is much faster by not having all validations.