How to generate survival data with time dependent covariates using R

I want to generate survival time from a Cox proportional hazards model that contains time dependent covariate. The model is

$h(t|X_i) =h_0(t) \exp(\gamma X_i + \alpha m_{i}(t))$

where $X_i$ is generated from Binomial(1,0.5) and $m_{i}(t)=\beta_0 + \beta_1 X_{i} + \beta_2 X_{i} t$.

The true parameter values are used as $\gamma = 1.5, \beta_0 = 0, \beta_1 = -1, \beta_2 = -1.5, h_0(t) = 1$

For time-independent covariate (i.e. $h(t|X_i) =h_0(t) \exp(\gamma X_i)$ I generated as follows

#For time independent case
# h_0(t) = 1
gamma <- -1
u <- runif(n=100,min=0,max=1)
Xi <- rbinom(n=100,size=1,prob=0.5)
T <- -log(u)/exp(gamma*Xi)

• What sort of function is $m_i(t)$? Is it continuous? Piecewise constant? A different algorithm will probably be needed accordingly. Jul 24 '14 at 20:42
• $m_{i}(t)$ is a time dependent covariate, for simplicity you can consider a proportional relationship with time. Jul 25 '14 at 13:42
• I've edited my question, considering a function of $m_i(t)$ Jul 29 '14 at 5:52
• how did you perform the R code from the above equation? means that At each death time within the same id the program needs to figure out what the covariates are for everyone which is either x is equal to 1 or 0. if all equal to 1 cumsum the hazard. after that calculate the survival function. lets it pick the right line for each subject. Jun 14 '15 at 4:47
• As Z. Zhang points out then have a look at this article. Further, you can see my answer to his question where I show how to simulate for those in the $X_i = 1$ group in R. Nov 12 '17 at 20:24

OK from your R code you are assuming an exponential distribution (constant hazard) for your baseline hazard. Your hazard functions are therefore:

$$h\left(t \mid X_i\right) = \begin{cases} \exp{\left(\alpha \beta_0\right)} & \text{if X_i = 0,} \\ \exp{\left(\gamma + \alpha\left(\beta_0+\beta_1+\beta_2 t\right)\right)} & \text{if X_i = 1.} \end{cases}$$

We then integrate these with respect to $t$ to get the cumulative hazard function:

\begin{align} \Lambda\left(t\mid X_i\right) &= \begin{cases} t \exp{\left(\alpha \beta_0\right)} & \text{if X_i=0,} \\ \int_0^t{\exp{\left(\gamma + \alpha\left(\beta_0+\beta_1+\beta_2 \tau\right)\right)} \,d\tau} & \text{if X_i=1.} \end{cases} \\ &= \begin{cases} t \exp{\left(\alpha \beta_0\right)} & \text{if X_i=0,} \\ \exp{\left(\gamma + \alpha\left(\beta_0+\beta_1\right)\right)} \frac{1}{\alpha\beta_2} \left(\exp\left(\alpha\beta_2 t\right)-1\right) & \text{if X_i=1.} \end{cases} \end{align}

These then give us the survival functions:

\begin{align} S\left(t\right) &= \exp{\left(-\Lambda\left(t\right)\right)} \\ &= \begin{cases} \exp{\left(-t \exp{\left(\alpha \beta_0\right)}\right)} & \text{if X_i=0,} \\ \exp{\left(-\exp{\left(\gamma + \alpha\left(\beta_0+\beta_1\right)\right)} \frac{1}{\alpha\beta_2} \left(\exp\left(\alpha\beta_2 t\right)-1\right)\right)} & \text{if X_i=1.} \end{cases} \end{align}

You then generate by sampling $X_i$ and $U\sim\mathrm{Uniform\left(0,1\right)}$, substituting $U$ for $S\left(t\right)$ and rearranging the appropriate formula (based on $X_i$) to simulate $t$. This should be straightforward algebra you can then code up in R but please let me know by comment if you need any further help.

• Thank you very much for the algebra. I will code up in R and will contact you for further help. Jul 29 '14 at 17:07