Merge pull request Merck#216 from jdblischak/rmst-formula

Create formula interface for RMST test

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: simtrial
 Type: Package
 Title: Clinical Trial Simulation
-Version: 0.3.2.10
+Version: 0.3.2.13
 Authors@R: c(
     person("Keaven", "Anderson", email = "[email protected]", role = c("aut")),
     person("Yilong", "Zhang", email = "[email protected]", role =  c("aut")),

R/maxcombo.R

@@ -36,7 +36,7 @@
 #' sim_pw_surv(n = 200) |>
 #'   cut_data_by_event(150) |>
 #'   maxcombo(rho = c(0, 0), gamma = c(0, 0.5))
-maxcombo <- function(data, rho, gamma){
+maxcombo <- function(data, rho, gamma) {
   stopifnot(
     is.numeric(rho), is.numeric(gamma),
     rho >= 0, gamma >= 0,
@@ -50,7 +50,7 @@ maxcombo <- function(data, rho, gamma){
       return_corr = TRUE
     )
 
-  ans <-  data.frame(p_value = pvalue_maxcombo(ans))
+  ans <- data.frame(p_value = pvalue_maxcombo(ans))
 
   return(ans)
 }

R/rmst.R

@@ -25,9 +25,25 @@
 #' @param var_label_tte Column name of the TTE variable.
 #' @param var_label_event Column name of the event variable.
 #' @param var_label_group Column name of the grouping variable.
+#' @param formula (default: `NULL`) A formula that indicates the TTE, event, and
+#'   group variables (in that exact order; see Details below). This is an
+#'   alternative to specifying the variables as strings. If a formula is
+#'   provided, the values passed to `var_label_tte`, `var_label_event`, and
+#'   `var_label_group` are ignored.
 #' @param reference A group label indicating the reference group.
 #' @param alpha Type I error.
 #'
+#' @details
+#' The argument `formula` is provided as a convenience to easily specify the TTE,
+#' event, and grouping variables. Note however that only the order of the three
+#' variables is actually used by the underlying function. Any functions applied
+#' in the formula are ignored, and thus should only be used for documenting your
+#' intent. For example, you can use the syntax from the survival package
+#' `Surv(tte | event) ~ group` to highlight the relation between the TTE and
+#' event variables, but the function `Surv()` is never actually executed.
+#' Importantly, you shouldn't apply any transformation functions such as `log()`
+#' since these will also be ignored.
+#'
 #' @return The z statistics.
 #'
 #' @export
@@ -42,14 +58,55 @@
 #'   tau = 10,
 #'   reference = "0"
 #' )
+#'
+#' # Formula interface
+#' library("survival")
+#'
+#' rmst(
+#'   data = ex1_delayed_effect,
+#'   formula = Surv(month | evntd) ~ trt,
+#'   tau = 10,
+#'   reference = "0"
+#' )
+#'
+#' # alternative
+#' rmst(
+#'   data = ex1_delayed_effect,
+#'   formula = ~ Surv(month, evntd, trt),
+#'   tau = 10,
+#'   reference = "0"
+#' )
+#'
+#' # This example doesn't make statistical sense, but demonstrates that only the
+#' # order of the 3 variables actually matters
+#' rmst(
+#'   data = ex1_delayed_effect,
+#'   formula = month ~ evntd + trt,
+#'   tau = 10,
+#'   reference = "0"
+#' )
 rmst <- function(
     data,
     tau = 10,
     var_label_tte = "tte",
     var_label_event = "event",
     var_label_group = "treatment",
+    formula = NULL,
     reference = "control",
     alpha = 0.05) {
+  stopifnot(is.data.frame(data))
+
+  if (!is.null(formula)) {
+    stopifnot(inherits(formula, "formula"))
+    variables <- colnames(stats::get_all_vars(formula = formula, data = data))
+    if (length(variables) != 3) {
+      stop("The formula interface requires exactly 3 variables specified")
+    }
+    var_label_tte <- variables[1]
+    var_label_event <- variables[2]
+    var_label_group <- variables[3]
+  }
+
   res <- rmst_two_arm(
     time_var = data[[var_label_tte]],
     event_var = data[[var_label_event]],

R/sim_gs_n.R

@@ -285,7 +285,7 @@ sim_gs_n <- function(
 
     for (i_analysis in seq_len(n_analysis)) {
       # Get cut date
-      cut_date[i_analysis] <- cutting[[i_analysis]](data = simu_data)
+      cut_date[i_analysis] <- cutting[[i_analysis]](simu_data)
 
       # Cut the data
       simu_data_cut <- simu_data |> cut_data_by_date(cut_date[i_analysis])

R/wlr.R

@@ -40,35 +40,35 @@
 #'   cut_data_by_event(150) |>
 #'   wlr(weight = early_zero(early_period = 4))
 #'
-wlr <- function(data, weight){
-
+wlr <- function(data, weight) {
   if (inherits(weight, "fh")) {
     ans <- data |>
       counting_process(arm = "experimental") |>
       fh_weight(rho_gamma = data.frame(rho = weight$rho, gamma = weight$gamma))
-
   } else if (inherits(weight, "mb")) {
     ans <- data |>
       counting_process(arm = "experimental") |>
       mb_weight(delay = weight$delay, w_max = weight$w_max)
     setDT(ans)
-    ans <- ans[,
-               .(
-                 s = sum(o_minus_e * mb_weight),
-                 v = sum(var_o_minus_e * mb_weight^2)
-               )
+    ans <- ans[
+      ,
+      .(
+        s = sum(o_minus_e * mb_weight),
+        v = sum(var_o_minus_e * mb_weight^2)
+      )
     ][, .(z = s / sqrt(v))]
     setDF(ans)
   } else if (inherits(weight, "early_period")) {
     ans <- data |>
       counting_process(arm = "experimental") |>
       early_zero_weight(early_period = weight$early_period)
     setDT(ans)
-    ans <- ans[,
-               .(
-                 s = sum(o_minus_e * weight),
-                 v = sum(var_o_minus_e * weight^2)
-               )
+    ans <- ans[
+      ,
+      .(
+        s = sum(o_minus_e * weight),
+        v = sum(var_o_minus_e * weight^2)
+      )
     ][, .(z = s / sqrt(v))]
     setDF(ans)
   }

R/wlr_weight.R

@@ -25,7 +25,7 @@
 #' @return A list of parameters of the Fleming-Harrington weighting function
 #' @examples
 #' fh(rho = 0, gamma = 0.5)
-fh <- function(rho = 0, gamma = 0){
+fh <- function(rho = 0, gamma = 0) {
   structure(list(rho = rho, gamma = gamma), class = c("list", "fh", "wlr"))
 }
 
@@ -42,7 +42,7 @@ fh <- function(rho = 0, gamma = 0){
 #'
 #' @examples
 #' mb(delay = 6, w_max = 2)
-mb <- function(delay = 4, w_max = Inf){
+mb <- function(delay = 4, w_max = Inf) {
   structure(list(delay = delay, w_max = w_max), class = c("list", "mb", "wlr"))
 }
 
@@ -59,6 +59,6 @@ mb <- function(delay = 4, w_max = Inf){
 #'
 #' @examples
 #' early_zero(6)
-early_zero <- function(early_period){
+early_zero <- function(early_period) {
   structure(list(early_period = early_period), class = c("list", "early_period", "wlr"))
 }

tests/testthat/helper-pmvnorm.R

@@ -137,7 +137,7 @@ ptvnorm <- function(h, r, ro) {
   hp2 <- (h2 * rr133 - h1 * f3 - h3 * f2) / fac / sqrt(rr133)
   TV <- w * exp((rr12 * h12 - h122) / rr122) / sqrt(rr122) * pnorm(hp1) *
     r12 + w * exp((rr13 * h13 - h132) / rr133) / sqrt(rr133) * pnorm(hp2) *
-    r13
+      r13
   TV <- sum(TV)
   rho <- matrix(c(1, r23, r23, 1), 2, 2)
   p2 <- mvtnorm::pmvnorm(-Inf, c(h2, h3), c(0, 0), rho)

tests/testthat/helper-sim_gs_n.R

@@ -3,24 +3,26 @@
 test_enroll_rate <- function() {
   # parameters for enrollment
   enroll_rampup_duration <- 4 # duration for enrollment ramp up
-  enroll_duration <- 16       # total enrollment duration
+  enroll_duration <- 16 # total enrollment duration
   enroll_rate <- gsDesign2::define_enroll_rate(
-    duration = c(enroll_rampup_duration,
-                 enroll_duration - enroll_rampup_duration),
+    duration = c(
+      enroll_rampup_duration,
+      enroll_duration - enroll_rampup_duration
+    ),
     rate = c(10, 30)
   )
   return(enroll_rate)
 }
 
 test_fail_rate <- function() {
   # parameters for treatment effect
-  delay_effect_duration <- 3  # delay treatment effect in months
-  median_ctrl <- 9            # survival median of the control arm
-  median_exp <- c(9, 14)      # survival median of the experimental arm
+  delay_effect_duration <- 3 # delay treatment effect in months
+  median_ctrl <- 9 # survival median of the control arm
+  median_exp <- c(9, 14) # survival median of the experimental arm
   dropout_rate <- 0.001
   fail_rate <- gsDesign2::define_fail_rate(
     duration = c(delay_effect_duration, 100),
-    fail_rate = log(2) /  median_ctrl,
+    fail_rate = log(2) / median_ctrl,
     hr = median_ctrl / median_exp,
     dropout_rate = dropout_rate
   )
@@ -29,9 +31,9 @@ test_fail_rate <- function() {
 
 test_cutting <- function() {
   # other related parameters
-  alpha <- 0.025              # type I error
-  beta <- 0.1                 # type II error
-  ratio <- 1                  # randomization ratio (exp:ctrl)
+  alpha <- 0.025 # type I error
+  beta <- 0.1 # type II error
+  ratio <- 1 # randomization ratio (exp:ctrl)
   # Define cuttings of 2 IAs and 1 FA
   # IA1
   # The 1st interim analysis will occur at the later of the following 3 conditions:

tests/testthat/helper-simfix.R

@@ -1,6 +1,6 @@
 # Helper functions used by test-double_programming_simfix.R
 
-test_simfix <- function () {
+test_simfix <- function() {
   # Study design using gsDesign
   alpha <- 0.025
   gamma <- c(5, 5, 47)