# Simulating survival data for a Cox model with a dichotomous time-varying covariate in R

I am interested in simulating non-proportional data following the formulation in Section 3.1 of the article: [Generating survival times to simulate Cox proportional hazards models with time‐varying covariates][1]. The problem is that the estimates for bt are not consistent with the true value of 1 and thus resulting in very low coverage probability. Can someone please help me, I've been working on this for a while now and I don't know what I'm doing wrong.

sim.50=function(n,b,bt,lambda,lambdac,alpha,theta){

#Generate covariates: x (time invariant), z (time-varying)
x=rnorm(n,0,1)
z=rbinom(n,1,0.5)

#Uniform random variable
u=runif(n,0,1)

t=(-log(u))/(lambda*exp(b*x))
t0=lambdac

t.e=ifelse(-log(u)<lambda*exp(b*x)*t0,t,
(-log(u)-lambda*exp(b*x)*t0+lambda*exp(b*x+bt)*t0)/(lambda*exp(b*x+bt)))

c=runif(n,1,theta )

time.e=pmin(t.e,c)
status.e=as.numeric(t.e<=c)

data.frame(x=x,z=z,u=u,t=t,t0=t0,t.e=t.e,c=c,time.e=time.e,status.e=status.e)

}

library(survival)
runs=1000

#Data storage for Cox model
cox=data.frame(matrix(nrow = runs, ncol = 2+2+2+2+2))
colnames(cox)=c("b", "bt","se_b","se_bt","cil_b","cil_bt", "ciu_b","ciu_bt", "coverage_b","coverage_bt" )

set.seed(11)
for(i in 1:runs){

#Data
data.e=sim.50(50,1,1,0.5,10,1.5,20.5)

#1.Fit the Cox proportional hazards model
fit.cox=coxph(Surv(time.e, status.e) ~ x+z,data =data.e)
#Obtain estimates, standard errors, 95% CI limits and coverage probability
cox[i,1:2]=fit.cox$$coef cox[i,3:4]=sqrt(diag(fit.cox$$var))
cox[i,5:6]=cox[i,1:2]+ qnorm(.025) * cox[i,3:4]
cox[i,7:8]=cox[i,1:2] + qnorm(.975) * cox[i,3:4]
cox[i,9:10]=as.integer ((cox[i,5:6]<=1) & (cox[i,7:8]>=1))

}

#Estimates
Estimates=colMeans(cox[,1:2])
Estimates
#Bias
true_values=c(1,1)
Bias=Estimates-true_values
Bias
#Estimated standard error
Est_Se=colMeans(cox[,3:4])
Est_Se
#Empirical standard error
Emp_Se=sqrt(diag(var(cox[,1:2])))
Emp_Se
#MSE
s=sum((cox[,1:2]-true_values)^2)
s
MSE=s/runs
MSE
#Coverage probability
cp=colMeans(cox[,9:10])
cp

[1]: https://onlinelibrary.wiley.com/doi/full/10.1002/sim.5452


To simulate from a time-varying Cox model you will need to simulate event times from the survival function: $$S_i(t) = \exp \biggl ( -\int_0^t h_0(s) \exp\{\gamma_1 x_i1 + \gamma_2 z_i(s)\} \, ds \biggr),$$ for a specific choice of the baseline hazard function $$h_0(t)$$.
You could, for example, use the inversion method in which first you simulate a value $$u^*$$ from a uniform distribution $$U(0, 1)$$, and then you need to find the event time $$T^*$$ for which $$S_i(T^*) = u^*$$. This last step in R could be done with function uniroot().