Solved – Simulation Analysis of a Cox Survival Model with Change Point.

Question

I wish to simulate survival data for the following probabilistic model which will be analyzed using a Cox model:

An exposure $X$ is modeled as binary having $X \sim_{iid} \mbox{Bernoulli}(p)$ where $p = \frac{1}{3}$. The hazard is a function of $X$ with:

\begin{eqnarray}
\log \left( \lambda(t|X=0) \right) &=& \beta_0\\
\log \left( \lambda(t|X=1) \right) &=& \beta_0 + \Delta\cdot\mathcal{I}(t \geq \tau)
\end{eqnarray}

I also wish to generate a censoring variable having an independent exponential distribution.

Answer 1

Constant hazards mean that survival times are exponentially distributed. It's possible, in the case of $X=1$ to simulate first whether your exponential time is greater than $\tau$ and then generate what that time will be according to the correct distribution (and, of course, memorylessness).

For an exponential RV $T \sim \mbox{Exp}(\beta_0)$ the probability that $T$ is greater than $\tau$ is just $\exp(-\theta\tau)$. If you compute the indicator of whether this exceeds that value, you can generate the actual value of $T$ using that $Pr(t \leq T | t \geq \tau) = 1-\exp( \left(t-\tau \right)\left(\theta + \Delta \right) )$

In R it looks like this:

n <- 100 ## set sample size
x <- rbinom(n, 1, 1/3) ## generate exposure
tau <- 1 ## set cut off for rate change point
theta <- 0.8 ## baseline rate
delta <- 0.3 ## rate difference after change point
eta <- 2 ## censoring rate
t <- numeric(n)  ## init failure times
t[x==0] <- rexp(sum(x==0), theta) ## generate unconditional failure times for unexposed
i <- rbinom(sum(x==1), 1, exp(-tau*theta)) ## generate indicator for change point
t[x==1] <- rexp(sum(x==1), theta + delta*i) + i*tau ## generate conditional failure times
r <- rexp(n, eta)
eventtm <- ifelse(t < r, t, r)
outcome <- ifelse(t < r, 1, 0)

coxph(Surv(eventtm, outcome) ~ x)