exp_cdfsim {cpsurvsim}

Inverse CDF simulation for the exponential change-point hazard distribution

Description

exp_cdfsim simulates time-to-event data from the exponential change-point hazard distribution by implementing the inverse CDF method.

Usage

exp_cdfsim(n, endtime, theta, tau = NA)

Arguments

n

Sample size

endtime

Maximum study time, point at which all participants are censored

theta

Scale parameter \theta

tau

Change-point(s) \tau

Details

This function simulates data for the exponential change-point hazard distribution with K change-points by simulating values of the exponential distribution and substituting them into the inverse hazard function. This method applies Type I right censoring at the endtime specified by the user. This function allows for up to four change-points.

Value

Dataset with n participants including a survival time and censoring indicator (0 = censored, 1 = event).

Examples

nochangepoint <- exp_cdfsim(n = 10, endtime = 20, theta = 0.05)
onechangepoint <- exp_cdfsim(n = 10, endtime = 20,
  theta = c(0.05, 0.01), tau = 10)
twochangepoints <- exp_cdfsim(n = 10, endtime = 20,
  theta = c(0.05, 0.01, 0.05), tau = c(8, 12))

# Pay attention to how you parameterize your model!
# This simulates a decreasing hazard
set.seed(7830)
decreasingHazard <- exp_cdfsim(n = 10, endtime = 20,
  theta = c(0.5, 0.2, 0.01), tau = c(8, 12))
# This tries to fit an increasing hazard, resulting in biased estimates
cp2.nll <- function(par, tau = tau, dta = dta){
  theta1 <- par[1]
  theta2 <- par[2]
  theta3 <- par[3]
  ll <- log(theta1) * sum(dta$time < tau[1])+
        log(theta2) * sum((tau[1] <= dta$time) * (dta$time < tau[2])) +
        log(theta3) * sum((dta$time >= tau[2]) * dta$censor) -
        theta1 * sum(dta$time * (dta$time < tau[1]) +
          tau[1] * (dta$time >= tau[1])) -
        theta2 * sum((dta$time - tau[1]) * (dta$time >= tau[1]) *
          (dta$time < tau[2]) + (tau[2] - tau[1]) * (dta$time >= tau[2])) -
        theta3 * sum((dta$time - tau[2]) * (dta$time >= tau[2]))
  return(-ll)
}
optim(par = c(0.001, 0.1, 0.5), fn = cp2.nll,
      tau = c(8, 12), dta = decreasingHazard)