weib_cdfsim {cpsurvsim}

Inverse CDF simulation for the Weibull change-point hazard distribution

Description

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

Usage

weib_cdfsim(n, endtime, gamma, theta, tau = NA)

Arguments

n

Sample size

endtime

Maximum study time, point at which all participants are censored

gamma

Shape parameter \gamma

theta

Scale parameter \theta

tau

Change-point(s) \tau

Details

This function simulates data from the Weibull 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 and \gamma is held constant.

Value

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

Examples

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

#' # Pay attention to how you parameterize your model!
# This simulates an increasing hazard
set.seed(9945)
increasingHazard <- weib_cdfsim(n = 100, endtime = 20, gamma = 2,
  theta = c(0.001, 0.005, 0.02), tau = c(8, 12))
# This tries to fit a decreasing hazard, resulting in biased estimates
cp2.nll <- function(par, tau = tau, gamma = gamma, dta = dta){
  theta1 <- par[1]
  theta2 <- par[2]
  theta3 <- par[3]
  ll <- (gamma - 1) * sum(dta$censor * log(dta$time)) +
    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/gamma) * sum((dta$time^gamma) * (dta$time < tau[1]) +
      (tau[1]^gamma) * (dta$time >= tau[1])) -
    (theta2/gamma) * sum((dta$time^gamma - tau[1]^gamma) *
      (dta$time >= tau[1]) * (dta$time<tau[2]) +
      (tau[2]^gamma - tau[1]^gamma) * (dta$time >= tau[2])) -
    (theta3/gamma) * sum((dta$time^gamma - tau[2]^gamma) *
      (dta$time >= tau[2]))
  return(-ll)
}
optim(par = c(0.2, 0.02, 0.01), fn = cp2.nll,
  tau = c(8, 12), gamma = 2,
  dta = increasingHazard)