0
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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().

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I'm still unclear on your suggestion, how does it tie in to my existing program? $\endgroup$ – ike Nov 27 '19 at 20:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.