69: Function to compute empirical power for a list of simulated trials (

#70)

Fixes #69

DESCRIPTION

@@ -24,14 +24,14 @@ Depends:
     R (>= 3.6)
 Imports:
     checkmate,
-    stats
+    stats,
+    survival
 Suggests:
     coxphw,
     knitr,
     mvna,
     prodlim,
     rmarkdown,
-    survival,
     testthat (>= 2.0)
 VignetteBuilder:
     knitr

NAMESPACE

@@ -32,9 +32,13 @@ export(getSumPCW)
 export(getTimePoint)
 export(getWaitTimeSum)
 export(integrateVector)
+export(logRankTest)
+export(passedLogRank)
 export(piecewise_exponential)
+export(powerEmp)
 export(pwA)
 export(trackEventsPerTrial)
 export(weibull_transition)
 import(checkmate)
 importFrom(stats,acf)
+importFrom(survival,survdiff)

R/hazardFunctions.R

@@ -75,11 +75,11 @@ avgHRExpOS <- function(transitionByArm, alpha = 0.5, upper = Inf) {
   h0 <- transitionByArm[[1]]$hazard
   h1 <- transitionByArm[[2]]$hazard
 
-  num <- integrate(avgHRIntegExpOS,
+  num <- stats::integrate(avgHRIntegExpOS,
     lower = 0, upper = upper, h01 = h1$h01, h02 = h1$h02, h12 = h1$h12,
     h0 = h0, h1 = h1, alpha = alpha
   )$value
-  denom <- integrate(avgHRIntegExpOS,
+  denom <- stats::integrate(avgHRIntegExpOS,
     lower = 0, upper = upper, h01 = h0$h01, h02 = h0$h02, h12 = h0$h12,
     h0 = h0, h1 = h1, alpha = alpha
   )$value

R/package.R

@@ -6,4 +6,5 @@
 
 #' @import checkmate
 #' @importFrom stats acf
+#' @importFrom survival survdiff
 NULL

R/powerEmp.R

@@ -0,0 +1,160 @@
+#' Log-Rank Test for a Single Trial
+#'
+#' This function evaluates the significance of either PFS or OS endpoints in a trial,
+#' based on a pre-specified critical value.
+#'
+#' @param data (`data.frame`)\cr data frame containing entry and exit times of an
+#' illness-death model. See [getSimulatedData()] for details.
+#' @param typeEvent (`string`)\cr endpoint to be evaluated, possible values are `PFS` and `OS`.
+#' @param critical (positive `number`)\cr critical value of the log-rank test.
+#'
+#' @return Logical value indicating log-rank test significance.
+#' @export
+#'
+#' @examples
+#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
+#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
+#' simTrial <- getClinicalTrials(
+#'   nRep = 1, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
+#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
+#'   accrual = list(param = "intensity", value = 7)
+#' )[[1]]
+#' logRankTest(data = simTrial, typeEvent = "OS", critical = 3.4)
+logRankTest <- function(data, typeEvent = c("PFS", "OS"), critical) {
+  # Must be have the format of one row per patient (`datType` must be 1rowPatient`
+  # in getClinicalTrials()), i.e. have 10 (if censored) or 11 columns.
+  assert_data_frame(data, min.cols = 10, max.cols = 11)
+  typeEvent <- match.arg(typeEvent)
+  assert_positive_number(critical)
+
+  if (typeEvent == "OS") {
+    time <- data$OStime
+    event <- data$OSevent
+  } else if (typeEvent == "PFS") {
+    time <- data$PFStime
+    event <- data$PFSevent
+  }
+
+  logRank <- survival::survdiff(survival::Surv(time, event) ~ trt, data)
+  sqrt(logRank$chisq) > critical
+}
+
+
+#' Helper function to conduct log-rank tests for either PFS or OS
+#'
+#' This function evaluates the significance of either PFS or OS endpoints for each trial
+#' in a list of trials, based on a pre-specified critical value.
+#'
+#' @param simTrials (`list`)\cr simulated trial data sets, see [getClinicalTrials()].
+#' @param typeEvent (`string`)\cr endpoint to be evaluated, possible values are `PFS` and `OS`.
+#' @param eventNum (`integer`)\cr number of events required to trigger analysis.
+#' @param critical (positive `number`)\cr critical value of the log-rank test.
+#'
+#' @return Logical vector indicating log-rank test significance for each trial.
+#' @export
+#'
+#' @examples
+#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
+#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
+#' simTrials <- getClinicalTrials(
+#'   nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
+#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
+#'   accrual = list(param = "intensity", value = 7)
+#' )
+#' passedLogRank(simTrials = simTrials, typeEvent = "PFS", eventNum = 300, critical = 2.4)
+#' @keywords internal
+passedLogRank <- function(simTrials, typeEvent, eventNum, critical) {
+  assert_list(simTrials, null.ok = FALSE)
+  assert_positive_number(critical)
+  assert_count(eventNum, positive = TRUE)
+
+  # Censor simulated trials at time-point of OS/PFS analysis.
+  trialsAna <- lapply(
+    X = simTrials,
+    FUN = censoringByNumberEvents,
+    eventNum = eventNum,
+    typeEvent = typeEvent
+  )
+  # Compute log-rank test for all trials for OS/PFS.
+  unlist(lapply(
+    X = trialsAna,
+    FUN = logRankTest,
+    typeEvent = typeEvent,
+    critical = critical
+  ))
+}
+
+#' Empirical Power for a List of Simulated Trials
+#'
+#' This function computes four types of empirical power — PFS, OS, at-least (significant
+#' in at least one of PFS/OS), and joint (significant in both PFS and OS) — using
+#' the log-rank test. Empirical power is calculated as the proportion of significant
+#' results in simulated trials, each ending when a set number of PFS/OS events occur.
+#' Critical values for PFS and OS test significance must be specified.
+#'
+#' @param simTrials (`list`)\cr simulated trial data sets, see [getClinicalTrials()].
+#' @param criticalPFS (positive `number`)\cr critical value of the log-rank test for PFS.
+#' @param criticalOS (positive `number`)\cr critical value of the log-rank test for OS.
+#' @param eventNumPFS (`integer`)\cr number of PFS events required to trigger PFS analysis.
+#' @param eventNumOS (`integer`)\cr number of OS events required to trigger OS analysis.
+#'
+#' @return This returns values of four measures of empirical power.
+#' @export
+#'
+#' @examples
+#' transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
+#' transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
+#' simTrials <- getClinicalTrials(
+#'   nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
+#'   transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
+#'   accrual = list(param = "intensity", value = 7)
+#' )
+#' powerEmp(
+#'   simTrials = simTrials, criticalPFS = 2.4, criticalOS = 2.2,
+#'   eventNumPFS = 300, eventNumOS = 500
+#' )
+powerEmp <- function(simTrials, criticalPFS, criticalOS, eventNumPFS, eventNumOS) {
+  assert_list(simTrials, null.ok = FALSE)
+  assert_positive_number(criticalPFS)
+  assert_positive_number(criticalOS)
+  assert_count(eventNumPFS, positive = TRUE)
+  assert_count(eventNumOS, positive = TRUE)
+
+  nRep <- length(simTrials)
+  simTrials <- lapply(
+    X = simTrials,
+    FUN = function(x) if (ncol(x) == 9) getDatasetWideFormat(x) else x
+  )
+
+  # Which trials passed the log-rank test for PFS/OS?
+  passedLogRankPFS <- passedLogRank(
+    simTrials = simTrials,
+    typeEvent = "PFS",
+    eventNum = eventNumPFS,
+    critical = criticalPFS
+  )
+  passedLogRankOS <- passedLogRank(
+    simTrials = simTrials,
+    typeEvent = "OS",
+    eventNum = eventNumOS,
+    critical = criticalOS
+  )
+
+  # Empirical power is the fraction of trials with significant log-rank test:
+  powerPFS <- sum(passedLogRankPFS) / nRep
+  powerOS <- sum(passedLogRankOS) / nRep
+
+  # Derived measures of power:
+  sumPassed <- passedLogRankPFS + passedLogRankOS
+  # At-least power: At least one of PFS/OS has significant log-rank test.
+  powerAtLeast <- sum(sumPassed >= 1) / nRep
+  # Joint power: Both PFS/OS have significant log-rank tests.
+  powerJoint <- sum(sumPassed == 2) / nRep
+
+  list(
+    "powerPFS" = powerPFS,
+    "powerOS" = powerOS,
+    "powerAtLeast" = powerAtLeast,
+    "powerJoint" = powerJoint
+  )
+}

_pkgdown.yaml

@@ -68,3 +68,7 @@ reference:
       - ExpHazOS
       - avgHRIntegExpOS
       - avgHRExpOS
+  - title: Empirical power
+    contents:
+      - logRankTest
+      - powerEmp

tests/testthat/test-powerEmp.R

@@ -0,0 +1,58 @@
+# logRankTest ----
+
+test_that("logRankTest works as expected", {
+  library(survival)
+  transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
+  transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
+  simTrial <- getClinicalTrials(
+    nRep = 1, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
+    transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
+    accrual = list(param = "intensity", value = 7)
+  )[[1]]
+  actual <- logRankTest(data = simTrial, typeEvent = "OS", critical = 3.4)
+  expect_equal(actual, TRUE)
+
+  actual2 <- logRankTest(data = simTrial, typeEvent = "PFS", critical = 6)
+  expect_equal(actual2, FALSE)
+})
+
+# passedLogRank ----
+
+test_that("passedLogRank works as expected", {
+  library(survival)
+  transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
+  transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
+  simTrials <- getClinicalTrials(
+    nRep = 3, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
+    transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
+    accrual = list(param = "intensity", value = 7)
+  )
+  actual <- passedLogRank(simTrials = simTrials, typeEvent = "PFS", eventNum = 300, critical = 2.4)
+  expect_equal(actual, c(TRUE, TRUE, FALSE))
+
+  actual2 <- passedLogRank(simTrials = simTrials, typeEvent = "OS", eventNum = 300, critical = 2.4)
+  expect_equal(actual2, c(FALSE, FALSE, FALSE))
+})
+
+# powerEmp ----
+
+test_that("powerEmp works as expected", {
+  library(survival)
+  transition1 <- exponential_transition(h01 = 0.06, h02 = 0.3, h12 = 0.3)
+  transition2 <- exponential_transition(h01 = 0.1, h02 = 0.4, h12 = 0.3)
+  simTrials <- getClinicalTrials(
+    nRep = 50, nPat = c(800, 800), seed = 1234, datType = "1rowPatient",
+    transitionByArm = list(transition1, transition2), dropout = list(rate = 0.5, time = 12),
+    accrual = list(param = "intensity", value = 7)
+  )
+  actual <- powerEmp(
+    simTrials = simTrials, criticalPFS = 2.4, criticalOS = 2.2,
+    eventNumPFS = 300, eventNumOS = 500
+  )
+  expect_equal(actual, list(
+    "powerPFS" = 0.74,
+    "powerOS" = 0.52,
+    "powerAtLeast" = 0.78,
+    "powerJoint" = 0.48
+  ))
+})