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

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?

• Before proceeding you would need to have justification for there being a discontinuity in the system. I.e., does a discrete external event happen at the supposed change point? If this is not the case, then it may be more reasonable to estimate the smooth relationship between time and the effect of $X$. This can be done using time-dependent covariables (e.g. a spline function in time interacting with $X$). You could start with a smooth scaled Schoenfeld residual plot in the context of the Cox proportional hazards model. – Frank Harrell Jul 16 '13 at 11:41
• @FrankHarrell I know your point, and in fact I want to use cox.zph function in order to compute Schoenfield residuals, but my main question resides on how to express this statement more formally in a mathematical way (although I'm not sure if it's already enough), and additionally, how to translate it in R code, because I don't know which object / function I can apply cox.zph. Many thanks for your consideration. – anxoestevez Jul 16 '13 at 12:07
• What are you trying to do here? Estimating $\lambda$? – boscovich Jul 16 '13 at 13:29
• No, I don't. The only two things I'm trying to do are well explained in the title: Formalize the statement (in terms of $S_x(t)$ and $S_0(t)$ in a Cox Regression Model) and code the expression in R in order to use it with coxph and cox.zph functions. – anxoestevez Jul 16 '13 at 13:47
• R survival and rms package cox.zph function when using "identity" as the transformation in plotting directly estimates $\beta(t)$. To put such an estimate into action would require a model to be fitted, such as a Cox model with time-dependent covariates, or a parametric version of same. – Frank Harrell Jul 16 '13 at 16:02

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?