How to simulate a Cox proportional hazards model with change point and code it in R

Question

I have a model that has the following characteristics:

• The covariate $X$ follows a $Be(1/3)$.
• If $X=0$, survival time $Y$ follows an $E=Exponential (1)$.
• If $X=1$, survival time $Y$ is generated as $E$ if $E\le\Psi$, and as $\Psi+E_\lambda$ if $E>\Psi$, where $\Psi$ is the 'change point' (let $\Psi=1$) and $E_\lambda=Exponential(\lambda)$ independent from $E$.

Does anyone know how to simulate this model and translate it in R code?

Answer 1

Well, after a little research I found the answer to my question, so here is the code that allows the simulation:

simulation <- function(n,lambda,changepoint, surv.df=TRUE) {

# Define the covariate, restrictions and parameters.

X <- rbinom(n,prob=1/3,size=1)
E <- rexp(n,1)
EL <- rexp(n,lambda)

# Define the piecewise function.

Y <- ifelse(X==0,E,
          ifelse(E<=changepoint,E,changepoint+EL))

## Construction of the data frame.

if (surv.df) data.frame(Y,X) else cbind(Y,X)
}

Now, if we would like to do m Monte Carlo replicas of this dataframe:

m = "number of replicas"

survive.df <- replicate(m, simulation(n = 100 ,lambda = 1,changepoint = 1), simplify=FALSE)

To view an i dataset:

survive.df[[i]]

Shape a Cox Proportional Hazards Model with the i dataset:

coxph(Surv(Y)~X , data=survive.df[[i]]) # Model.

cox.zph(coxph(Surv(Y)~X , data=survive.df[[i]])) # Schoenfield residuals.

Any suggestions?