Merge pull request #174 from Merck/rpwexp

Use inverse CDF method C++ implementation in `rpwexp()`

DESCRIPTION

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

NAMESPACE

@@ -12,9 +12,7 @@ export(mb_weight)
 export(pvalue_maxcombo)
 export(randomize_by_fixed_block)
 export(rpwexp)
-export(rpwexpRcpp)
 export(rpwexp_enroll)
-export(rpwexpinvRcpp)
 export(sim_fixed_n)
 export(sim_pw_surv)
 export(to_sim_pw_surv)

R/RcppExports.R

@@ -1,23 +1,25 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
-#' The piecewise exponential distribution in C++
+#' The piecewise exponential distribution using the inverse CDF method
 #'
 #' @param n Number of observations to be generated.
 #' @param fail_rate A data frame containing `duration` and `rate` variables.
 #'
-#' @export
-rpwexpRcpp <- function(n, fail_rate) {
-    .Call(`_simtrial_rpwexpRcpp`, n, fail_rate)
+#' @noRd
+#'
+rpwexp_inverse_cdf_cpp <- function(n, fail_rate) {
+    .Call(`_simtrial_rpwexp_inverse_cdf_cpp`, n, fail_rate)
 }
 
-#' The piecewise exponential distribution using inverse CDF method in C++
+#' The piecewise exponential distribution using the naive method
 #'
 #' @param n Number of observations to be generated.
 #' @param fail_rate A data frame containing `duration` and `rate` variables.
 #'
-#' @export
-rpwexpinvRcpp <- function(n, fail_rate) {
-    .Call(`_simtrial_rpwexpinvRcpp`, n, fail_rate)
+#' @noRd
+#'
+rpwexp_naive_cpp <- function(n, fail_rate) {
+    .Call(`_simtrial_rpwexp_naive_cpp`, n, fail_rate)
 }
 

R/rpwexp.R

@@ -70,51 +70,5 @@
 rpwexp <- function(
     n = 100,
     fail_rate = data.frame(duration = c(1, 1), rate = c(10, 20))) {
-  n_rate <- nrow(fail_rate)
-
-  if (n_rate == 1) {
-    # Set failure time to Inf if 0 failure rate
-    if (fail_rate$rate == 0) {
-      ans <- rep(Inf, n)
-    } else {
-      # Generate exponential failure time if non-0 failure rate
-      ans <- stats::rexp(n, fail_rate$rate)
-    }
-  } else {
-    # Start of first failure rate interval
-    start_time <- 0
-    # Ends of failure rate interval
-    end_time <- cumsum(fail_rate$duration)
-    # Initiate vector for failure times
-    ans <- rep(0, n)
-    # Index for event times not yet reached
-    indx <- rep(TRUE, n)
-
-    for (i in 1:n_rate) {
-      # Number of event times left to generate
-      n_event_left <- sum(indx)
-
-      # Stop if event is arrived
-      if (n_event_left == 0) {
-        break
-      }
-
-      # Set failure time to Inf for interval i if 0 fail rate
-      if (fail_rate$rate[i] == 0) {
-        ans[indx] <- start_time + rep(Inf, n_event_left)
-      } else {
-        # Generate exponential failure time for interval i if non-0 failure rate
-        ans[indx] <- start_time + stats::rexp(n_event_left, fail_rate$rate[i])
-      }
-
-      # Skip this for last interval as all remaining times are generated there
-      if (i < n_rate) {
-        start_time <- end_time[i]
-        # Update index of event times not yet reached
-        indx <- (ans > end_time[i])
-      }
-    }
-  }
-
-  ans
+  rpwexp_inverse_cdf_cpp(n = n, fail_rate = fail_rate)
 }

R/rpwexp_naive.R

@@ -0,0 +1,120 @@
+#  Copyright (c) 2023 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the simtrial program.
+#
+#  simtrial is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' The piecewise exponential distribution
+#'
+#' The piecewise exponential distribution allows a simple method to specify
+#' a distribution where the hazard rate changes over time.
+#' It is likely to be useful for conditions where failure rates change,
+#' but also for simulations where there may be a delayed treatment effect
+#' or a treatment effect that that is otherwise changing
+#' (for example, decreasing) over time.
+#' `rpwexp_naive()` is to support simulation of both the Lachin and Foulkes (1986)
+#' sample size method for (fixed trial duration) as well as the
+#' Kim and Tsiatis (1990) method (fixed enrollment rates and either
+#' fixed enrollment duration or fixed minimum follow-up);
+#' see [gsDesign::nSurv()].
+#'
+#' Using the `cumulative = TRUE` option, enrollment times that piecewise
+#' constant over time can be generated.
+#'
+#' @param n Number of observations to be generated.
+#' @param fail_rate A data frame containing `duration` and `rate` variables.
+#'   `rate` specifies failure rates during the corresponding interval duration
+#'   specified in `duration`. The final interval is extended to be infinite
+#'   to ensure all observations are generated.
+#'
+#' @noRd
+#'
+#' @examples
+#' # Example 1
+#' # Exponential failure times
+#' x <- rpwexp_naive(
+#'   n = 10000,
+#'   fail_rate = data.frame(rate = 5, duration = 1)
+#' )
+#' plot(sort(x), (10000:1) / 10001,
+#'   log = "y", main = "Exponential simulated survival curve",
+#'   xlab = "Time", ylab = "P{Survival}"
+#' )
+#'
+#' # Example 2
+#'
+#' # Get 10k piecewise exponential failure times.
+#' # Failure rates are 1 for time 0 to 0.5, 3 for time 0.5 to 1, and 10 for > 1.
+#' # Intervals specifies duration of each failure rate interval
+#' # with the final interval running to infinity.
+#' x <- rpwexp_naive(
+#'   n = 1e4,
+#'   fail_rate = data.frame(rate = c(1, 3, 10), duration = c(.5, .5, 1))
+#' )
+#' plot(sort(x), (1e4:1) / 10001,
+#'   log = "y", main = "PW Exponential simulated survival curve",
+#'   xlab = "Time", ylab = "P{Survival}"
+#' )
+rpwexp_naive <- function(
+    n = 100,
+    fail_rate = data.frame(duration = c(1, 1), rate = c(10, 20))) {
+  n_rate <- nrow(fail_rate)
+
+  if (n_rate == 1) {
+    # Set failure time to Inf if 0 failure rate
+    if (fail_rate$rate == 0) {
+      ans <- rep(Inf, n)
+    } else {
+      # Generate exponential failure time if non-0 failure rate
+      ans <- stats::rexp(n, fail_rate$rate)
+    }
+  } else {
+    # Start of first failure rate interval
+    start_time <- 0
+    # Ends of failure rate interval
+    end_time <- cumsum(fail_rate$duration)
+    # Initiate vector for failure times
+    ans <- rep(0, n)
+    # Index for event times not yet reached
+    indx <- rep(TRUE, n)
+
+    for (i in 1:n_rate) {
+      # Number of event times left to generate
+      n_event_left <- sum(indx)
+
+      # Stop if event is arrived
+      if (n_event_left == 0) {
+        break
+      }
+
+      # Set failure time to Inf for interval i if 0 fail rate
+      if (fail_rate$rate[i] == 0) {
+        ans[indx] <- start_time + rep(Inf, n_event_left)
+      } else {
+        # Generate exponential failure time for interval i if non-0 failure rate
+        ans[indx] <- start_time + stats::rexp(n_event_left, fail_rate$rate[i])
+      }
+
+      # Skip this for last interval as all remaining times are generated there
+      if (i < n_rate) {
+        start_time <- end_time[i]
+        # Update index of event times not yet reached
+        indx <- (ans > end_time[i])
+      }
+    }
+  }
+
+  ans
+}

R/sim_pw_surv.R

@@ -174,11 +174,11 @@ sim_pw_surv <- function(
   for (sr in unique_stratum) {
     for (tr in unique_treatment) {
       indx <- x$stratum == sr & x$treatment == tr
-      x$fail_time[indx] <- rpwexpinvRcpp(
+      x$fail_time[indx] <- rpwexp(
         n = sum(indx),
         fail_rate = fail_rate[fail_rate$stratum == sr & fail_rate$treatment == tr, , drop = FALSE]
       )
-      x$dropout_time[indx] <- rpwexpinvRcpp(
+      x$dropout_time[indx] <- rpwexp(
         n = sum(indx),
         fail_rate = dropout_rate[dropout_rate$stratum == sr & dropout_rate$treatment == tr, , drop = FALSE]
       )

_pkgdown.yml

@@ -27,8 +27,6 @@ reference:
 - title: "Distributions"
   contents:
     - rpwexp
-    - rpwexpRcpp
-    - rpwexpinvRcpp
     - rpwexp_enroll
     - fit_pwexp
 

src/RcppExports.cpp

@@ -10,34 +10,34 @@ Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
 Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
 #endif
 
-// rpwexpRcpp
-NumericVector rpwexpRcpp(int n, DataFrame fail_rate);
-RcppExport SEXP _simtrial_rpwexpRcpp(SEXP nSEXP, SEXP fail_rateSEXP) {
+// rpwexp_inverse_cdf_cpp
+NumericVector rpwexp_inverse_cdf_cpp(int n, DataFrame fail_rate);
+RcppExport SEXP _simtrial_rpwexp_inverse_cdf_cpp(SEXP nSEXP, SEXP fail_rateSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< int >::type n(nSEXP);
     Rcpp::traits::input_parameter< DataFrame >::type fail_rate(fail_rateSEXP);
-    rcpp_result_gen = Rcpp::wrap(rpwexpRcpp(n, fail_rate));
+    rcpp_result_gen = Rcpp::wrap(rpwexp_inverse_cdf_cpp(n, fail_rate));
     return rcpp_result_gen;
 END_RCPP
 }
-// rpwexpinvRcpp
-NumericVector rpwexpinvRcpp(int n, DataFrame fail_rate);
-RcppExport SEXP _simtrial_rpwexpinvRcpp(SEXP nSEXP, SEXP fail_rateSEXP) {
+// rpwexp_naive_cpp
+NumericVector rpwexp_naive_cpp(int n, DataFrame fail_rate);
+RcppExport SEXP _simtrial_rpwexp_naive_cpp(SEXP nSEXP, SEXP fail_rateSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
     Rcpp::traits::input_parameter< int >::type n(nSEXP);
     Rcpp::traits::input_parameter< DataFrame >::type fail_rate(fail_rateSEXP);
-    rcpp_result_gen = Rcpp::wrap(rpwexpinvRcpp(n, fail_rate));
+    rcpp_result_gen = Rcpp::wrap(rpwexp_naive_cpp(n, fail_rate));
     return rcpp_result_gen;
 END_RCPP
 }
 
 static const R_CallMethodDef CallEntries[] = {
-    {"_simtrial_rpwexpRcpp", (DL_FUNC) &_simtrial_rpwexpRcpp, 2},
-    {"_simtrial_rpwexpinvRcpp", (DL_FUNC) &_simtrial_rpwexpinvRcpp, 2},
+    {"_simtrial_rpwexp_inverse_cdf_cpp", (DL_FUNC) &_simtrial_rpwexp_inverse_cdf_cpp, 2},
+    {"_simtrial_rpwexp_naive_cpp", (DL_FUNC) &_simtrial_rpwexp_naive_cpp, 2},
     {NULL, NULL, 0}
 };