74: closed form calculations for PWCsurvOS (#97)

Co-authored-by: github-actions <41898282+github-actions[bot]@users.noreply.github.com>

Author: danielinteractive

NAMESPACE

@@ -24,7 +24,6 @@ export(ExpSurvPFS)
 export(PCWInversionMethod)
 export(PWCsurvOS)
 export(PWCsurvPFS)
-export(PwcOSInt)
 export(WeibOSInteg)
 export(WeibSurvOS)
 export(WeibSurvPFS)
@@ -50,7 +49,7 @@ export(getNumberEvents)
 export(getOneClinicalTrial)
 export(getOneToTwoRows)
 export(getPCWDistr)
-export(getPCWHazard)
+export(getPWCHazard)
 export(getResults)
 export(getSimulatedData)
 export(getSumPCW)

NEWS.md

@@ -10,8 +10,14 @@
 -   `ExpSurvOS` now returns 0 instead of NaN for large values of t.
 -   `WeibSurvOS` now does not return an error for large values of t.
 -   `PWCSurvOS` now does not return an error for large values of t.
+-   `PWCSurvOS` no longer returns values larger than 1, and is significantly faster, based on a closed form calculation instead of numerical integration.
 -   `getSimulatedData` now also works when there are no transitions from progression to death, similarly for `getOneClinicalTrial` (which now warns if there are no such transitions at all).
 
+### Miscellaneous
+
+-   Renamed piecewise constant hazards function to `getPWCHazard` (previously `getPCWHazard`).
+-   `PwcOSInt` is no longer exported, and only used for internal tests.
+
 # simIDM 0.0.5
 
 -   First CRAN version of the package.

R/corPFSOS.R

@@ -16,7 +16,7 @@ survPFS <- function(transition, t) {
   UseMethod("survPFS")
 }
 
-#' @describeIn PFS Survival Function for an exponential transition model.
+#' @describeIn survPFS Survival Function for an exponential transition model.
 #' @export
 #'
 #' @examples
@@ -26,7 +26,7 @@ survPFS.ExponentialTransition <- function(transition, t) {
   ExpSurvPFS(t = t, h01 = transition$hazards$h01, h02 = transition$hazards$h02)
 }
 
-#' @describeIn PFS Survival Function for a Weibull transition model.
+#' @describeIn survPFS Survival Function for a Weibull transition model.
 #' @export
 #'
 #' @examples
@@ -39,7 +39,7 @@ survPFS.WeibullTransition <- function(transition, t) {
   )
 }
 
-#' @describeIn PFS Survival Function for a piecewise constant transition model.
+#' @describeIn survPFS Survival Function for a piecewise constant transition model.
 #' @export
 #'
 #' @examples
@@ -73,7 +73,7 @@ survOS <- function(transition, t) {
   UseMethod("survOS")
 }
 
-#' @describeIn OS Survival Function for an exponential transition model.
+#' @describeIn survOS Survival Function for an exponential transition model.
 #' @export
 #'
 #' @examples
@@ -86,7 +86,7 @@ survOS.ExponentialTransition <- function(transition, t) {
   )
 }
 
-#' @describeIn OS Survival Function for a Weibull transition model.
+#' @describeIn survOS Survival Function for a Weibull transition model.
 #' @export
 #'
 #' @examples
@@ -100,7 +100,7 @@ survOS.WeibullTransition <- function(transition, t) {
   )
 }
 
-#' @describeIn OS Survival Function for a piecewise constant transition model.
+#' @describeIn survOS Survival Function for a piecewise constant transition model.
 #' @export
 #'
 #' @examples

R/estimateParams.R

@@ -148,7 +148,7 @@ haz.WeibullTransition <- function(transition, t, trans) {
 #' )
 #' haz(transition, 6, 2)
 haz.PWCTransition <- function(transition, t, trans) {
-  getPCWHazard(unlist(transition$hazards[trans], use.names = FALSE),
+  getPWCHazard(unlist(transition$hazards[trans], use.names = FALSE),
     unlist(transition$intervals[trans], use.names = FALSE),
     x = t
   )

R/getSimulatedData.R

@@ -204,9 +204,9 @@ getSimulatedData <- function(N,
     Sum_PCW <- getSumPCW(h$h01, h$h02, pw$pw01, pw$pw02)
     wait_time <- getPCWDistr(U, Sum_PCW$hazards, Sum_PCW$intervals, entry)
     # A binomial experiment decides on death or progression.
-    numerator <- getPCWHazard(h$h01, pw$pw01, wait_time)
-    denumerator <- getPCWHazard(h$h01, pw$pw01, wait_time) +
-      getPCWHazard(h$h02, pw$pw02, wait_time)
+    numerator <- getPWCHazard(h$h01, pw$pw01, wait_time)
+    denumerator <- getPWCHazard(h$h01, pw$pw01, wait_time) +
+      getPWCHazard(h$h02, pw$pw02, wait_time)
   }
 
   to_prob <- stats::rbinom(N, 1, numerator / denumerator)

R/piecewiseHazards.R

@@ -10,8 +10,8 @@
 #' @export
 #'
 #' @examples
-#' getPCWHazard(c(1, 1.2, 1.4), c(0, 2, 3), c(1, 4, 6))
-getPCWHazard <- function(haz, pw, x) { # nolint
+#' getPWCHazard(c(1, 1.2, 1.4), c(0, 2, 3), c(1, 4, 6))
+getPWCHazard <- function(haz, pw, x) { # nolint
   sapply(x, function(jj) {
     y <- NULL
     # Find interval and corresponding hazard value for time x[jj].
@@ -44,8 +44,8 @@ getSumPCW <- function(haz1, haz2, pw1, pw2) {
   haz_sum <- NULL
   # Get sum of hazards for all intervals.
   for (i in seq_along(cuts_sum)) {
-    haz_sum[i] <- getPCWHazard(haz1, pw1, cuts_sum[i]) +
-      getPCWHazard(haz2, pw2, cuts_sum[i])
+    haz_sum[i] <- getPWCHazard(haz1, pw1, cuts_sum[i]) +
+      getPWCHazard(haz2, pw2, cuts_sum[i])
   }
   list(
     hazards = haz_sum,

R/survivalFunctions.R

@@ -152,7 +152,7 @@ pwA <- function(t, haz, pw) {
   # Time intervals: time-points + cut points.
   pw_new <- c(pw[pw < max(t)])
   pw_time <- sort(unique(c(pw_new, t)))
-  haz_all <- getPCWHazard(haz, pw, pw_time)
+  haz_all <- getPWCHazard(haz, pw, pw_time)
 
   # Cumulative hazard.
   i <- 1:(length(pw_time) - 1)
@@ -193,13 +193,11 @@ PWCsurvPFS <- function(t, h01, h02, pw01, pw02) {
 #'
 #' @return Numeric results of the integrand used to calculate
 #' the OS survival function for piecewise constant transition hazards, see  `PWCsurvOS`.
-#' @export
 #'
-#' @examples
-#' PwcOSInt(1:5, 6, c(0.3, 0.5), c(0.5, 0.8), c(0.7, 1), c(0, 4), c(0, 8), c(0, 3))
+#' @keywords internal
 PwcOSInt <- function(x, t, h01, h02, h12, pw01, pw02, pw12) {
   PWCsurvPFS(x, h01, h02, pw01, pw02) *
-    getPCWHazard(h01, pw01, x) *
+    getPWCHazard(h01, pw01, x) *
     exp(pwA(x, h12, pw12) - pwA(t, h12, pw12))
 }
 
@@ -215,7 +213,6 @@ PwcOSInt <- function(x, t, h01, h02, h12, pw01, pw02, pw12) {
 #'
 #' @return This returns the value of OS survival function at time t.
 #' @export
-
 #'
 #' @examples
 #' PWCsurvOS(1:5, c(0.3, 0.5), c(0.5, 0.8), c(0.7, 1), c(0, 4), c(0, 8), c(0, 3))
@@ -228,19 +225,75 @@ PWCsurvOS <- function(t, h01, h02, h12, pw01, pw02, pw12) {
   assert_intervals(pw02, length(h02))
   assert_intervals(pw12, length(h12))
 
-  PWCsurvPFS(t, h01, h02, pw01, pw02) +
-    sapply(t, function(t) {
-      integrateVector(PwcOSInt,
-        upper = t,
-        t = t,
-        h01 = h01,
-        h02 = h02,
-        h12 = h12,
-        pw01 = pw01,
-        pw02 = pw02,
-        pw12 = pw12
+  # Assemble unique time points and corresponding hazard values once.
+  unique_pw_times <- sort(unique(c(pw01, pw02, pw12)))
+  h01_at_times <- getPWCHazard(h01, pw01, unique_pw_times)
+  h02_at_times <- getPWCHazard(h02, pw02, unique_pw_times)
+  h12_at_times <- getPWCHazard(h12, pw12, unique_pw_times)
+
+  # We work for the integral parts with sorted and unique time points.
+  t_sorted <- sort(unique(t))
+  t_sorted_zero <- c(0, t_sorted)
+  n_t <- length(t_sorted)
+  int_sums <- numeric(n_t)
+  cum_haz_12 <- pwA(t_sorted, h12, pw12)
+  for (i in seq_len(n_t)) {
+    t_start <- t_sorted_zero[i]
+    t_end <- t_sorted_zero[i + 1]
+    # Determine the indices of the time intervals we are working with here.
+    start_index <- findInterval(t_start, unique_pw_times)
+    # We use here left.open = TRUE, so that when t_end is identical to a change point,
+    # then we don't need an additional part.
+    end_index <- max(findInterval(t_end, unique_pw_times, left.open = TRUE), 1L)
+    index_bounds <- start_index:end_index
+    # Determine corresponding time intervals.
+    time_bounds <- if (length(index_bounds) == 1L) {
+      rbind(c(t_start, t_end))
+    } else if (length(index_bounds) == 2L) {
+      rbind(
+        c(t_start, unique_pw_times[end_index]),
+        c(unique_pw_times[end_index], t_end)
       )
-    })
+    } else {
+      n_ind_bounds <- length(index_bounds)
+      rbind(
+        c(t_start, unique_pw_times[start_index + 1L]),
+        cbind(
+          unique_pw_times[index_bounds[-c(1, n_ind_bounds)]],
+          unique_pw_times[index_bounds[-c(1, 2)]]
+        ),
+        c(unique_pw_times[end_index], t_end)
+      )
+    }
+    for (j in seq_along(index_bounds)) {
+      time_j <- time_bounds[j, 1]
+      time_jp1 <- time_bounds[j, 2]
+      intercept_a <- pwA(time_j, h12, pw12) - pwA(time_j, h01, pw01) - pwA(time_j, h02, pw02)
+      slope_b <- h12_at_times[index_bounds[j]] - h01_at_times[index_bounds[j]] - h02_at_times[index_bounds[j]]
+      h01_j <- h01_at_times[index_bounds[j]]
+      int_js <- if (slope_b != 0) {
+        h01_j / slope_b * exp(intercept_a) *
+          (exp(slope_b * (time_jp1 - time_j) - cum_haz_12[i:n_t]) - exp(-cum_haz_12[i:n_t]))
+      } else {
+        h01_j * exp(intercept_a - cum_haz_12[i:n_t]) * (time_jp1 - time_j)
+      }
+      int_sums[i:n_t] <- int_sums[i:n_t] + int_js
+    }
+  }
+
+  # Match back to all time points,
+  # which might not be unique or ordered.
+  int_sums <- int_sums[match(t, t_sorted)]
+  result <- PWCsurvPFS(t, h01, h02, pw01, pw02) + int_sums
+
+  # Cap at 1 in a safe way - i.e. first check.
+  assert_numeric(result, finite = TRUE, any.missing = FALSE)
+  above_one <- result > 1
+  if (any(above_one)) {
+    assert_true(all((result[above_one] - 1) < sqrt(.Machine$double.eps)))
+    result[above_one] <- 1
+  }
+  result
 }
 
 #' Helper Function for Single Quantile for OS Survival Function

_pkgdown.yaml

@@ -38,7 +38,7 @@ reference:
       - weibull_transition
   - title: Helper functions for piecewise or Weibull Hazards
     contents:
-      - getPCWHazard
+      - getPWCHazard
       - getSumPCW
       - getPCWDistr
       - PCWInversionMethod
@@ -72,7 +72,6 @@ reference:
       - WeibSurvOS
       - pwA
       - PWCsurvPFS
-      - PwcOSInt
       - PWCsurvOS
   - title: Hazard functions
     contents:

design/design_getSimulatedData.Rmd

@@ -261,9 +261,9 @@ getSimulatedData <- function(N,
       unique(sort(c(pw$pw01, pw$pw02))), entry
     )
     ### binomial experiment decides on death or progression
-    numerator <- getPCWHazard(h$h01, pw$pw01, wait_time)
-    denumerator <- getPCWHazard(h$h01, pw$pw01, wait_time) +
-                   getPCWHazard(h$h02, pw$pw02, wait_time)
+    numerator <- getPWCHazard(h$h01, pw$pw01, wait_time)
+    denumerator <- getPWCHazard(h$h01, pw$pw01, wait_time) +
+                   getPWCHazard(h$h02, pw$pw02, wait_time)
   }
 
   to_prob <- rbinom(N, 1, numerator / denumerator)
@@ -300,7 +300,7 @@ getSimulatedData <- function(N,
 `getsimulateData` uses the following helper functions:
 
 - `getOneTwoTransitionRows`: generates rows with transition times for patients in the transient state progression, i.e. entry and exit times for 1-> 2 or 1-> cens transitions
-- `getPCWHazard`: for piecewise-constant hazards: function to get hazard value at time x
+- `getPWCHazard`: for piecewise-constant hazards: function to get hazard value at time x
 - `getSumPCW`: calculates the sum of two piece-wise constant hazards
 - `getPCWDistr`: generates event times with a distribution resulting
             from specified piece-wise constant hazards
@@ -378,7 +378,7 @@ addStaggeredEntry <- function(simData, N, accrualParam, accrualValue) {
 ### or Weibull distributed event times
 
 ## for piecewise-constant hazards: function to get hazard value at time x
-getPCWHazard <- function(haz, pw, x) {
+getPWCHazard <- function(haz, pw, x) {
   hazVal <- sapply(x, function(jj) {
     y <- NULL
     # find interval and correponding hazard value for time x[jj]
@@ -401,8 +401,8 @@ getSumPCW <- function(haz1, haz2, pw1, pw2) {
   haz_sum <- NULL
   ## get sum of hazards for all intervals
   for (i in 1:(length(cuts_sum))) {
-    haz_sum[i] <- getPCWHazard(haz1, pw1, cuts_sum[i]) +
-                  getPCWHazard(haz2, pw2, cuts_sum[i])
+    haz_sum[i] <- getPWCHazard(haz1, pw1, cuts_sum[i]) +
+                  getPWCHazard(haz2, pw2, cuts_sum[i])
   }
   return(haz_sum)
 }

inst/WORDLIST

@@ -25,17 +25,22 @@ calender
 cdot
 censAct
 doi
+dotsb
+dotsc
+du
 du
 entryAct
 exitAct
 frac
 funder
+geq
 ik
 infty
 integrand
 lapply
 mathbb
 multistate
+neq
 pre
 recruitTime
 rightarrow