Quantile Function of Piecewise Exponential Distribution

To generate endpoint that follows a piecewise exponential distribution and is correlated with other endpoints, the copula method is commonly used. The quantile function needs to be specified (e.g., in the simdata package). If a piecewise exponential distributed endpoint is independent to other endpoints, one can simply use TrialSimulator::PiecewiseConstantExponentialRNG() to specify the generator argument in endpoint().

There are many R packages implementing the quantile function of the piecewise exponential distributed random variable. Why do I implement it again? The reason is that this function is extremely important for simulating time-to-event endpoint in clinical trial simulation, thus the speed matters. qPiecewiseExponential() is implemented purely in R for code transparency, and is much faster than other packages.

Usage

qPiecewiseExponential(p, times, piecewise_risk)

Arguments

p

numeric. A vector of probabilities.

times

numeric. A vector of time points where risk (event rates) change. 0 shouldn't be in times.

piecewise_risk

numeric. A vector of constant risk (event rates) in a time window. length(piecewise_risk) = length(times) + 1. To fully specify a piecewise exponential distribution, the number of risk parameters is one greater than the number of changepoints in times.

Value

A vector of quantiles.

Examples

## This code snippet can take > 10s to execute
if(interactive()){
library(TrialSimulator)

run <- TRUE
if(!requireNamespace("rpact", quietly = TRUE)){
  run <- FALSE
  message("Please install 'rpact' to run this example.")
}

if(!requireNamespace("PWEXP", quietly = TRUE)){
  run <- FALSE
  message("Please install 'PWEXP' to run this example.")
}

if(!requireNamespace("PWEALL", quietly = TRUE)){
  run <- FALSE
  message("Please install 'PWEALL' to run this example.")
}

if(run){

x <- solvePiecewiseConstantExponentialDistribution(
  surv_prob = c(.9, .75, .64, .42, .28),
  times = c(.4, 1.2, 4, 5.5, 9)
)

p <- stats::runif(1e5)

## fast, and only five lines of R codes
message('TrialSimulator::qPiecewiseExponential(): ')
system.time(
  a <- qPiecewiseExponential(
    p, times = x$end_time, piecewise_risk = c(x$piecewise_risk, .1)
  )
) |> print()

## > 10s
message('rpact::getPiecewiseExponentialQuantile(): ')
system.time(
  b <- rpact::getPiecewiseExponentialQuantile(
    p, piecewiseSurvivalTime = c(0, x$end_time),
    piecewiseLambda=c(x$piecewise_risk, .1)
  )
) |> print()

## equally fast, but implemented in Fortran
message('PWEALL::qpwe(): ')
system.time(
  c <- PWEALL::qpwe(p, rate = c(x$piecewise_risk, .1), tchange = c(0, x$end_time))$q
) |> print()

## equally fast, long codes in R (maybe more versatile?)
message('PWEXP::qpwexp(): ')
system.time(
  d <- PWEXP::qpwexp(p, rate = c(x$piecewise_risk, .1), breakpoint = x$end_time)
) |> print()

message('a == b: ')
all.equal(a, b) |> print()
message('a == c: ')
all.equal(a, c) |> print()
message('a == d: ')
all.equal(a, d) |> print()

}
}