0
$\begingroup$

I'm trying to simulate a survival data with very large sample size and show that the KM estimators approach the true survival function, however it ends up that the two are visually quite different. I think I wasn't comparing the right thing but I couldn't find where I did wrong.

First I simulate a survival data with 100K observations, max follow-up time 50, censoring rate 0.8, and specify the baseline hazard function. I use the Weibull distribution but go for the easy way that lambda=0.02 and rho=1 so that baseline hazard is a constant:

library(coxed)
library(KMsurv)

set.seed(20210330)
nsample=1e5
lambda=0.02
rho=1
true.beta=0.4

myhaz=function(t){
  # lambda*rho*(t^(rho-1))
  0.02
}
simdata <- sim.survdata(N=nsample, T=50, censor=.8, 
                        num.data.frames = 1, beta = 0.4, 
                        hazard.fun=myhaz,
                        X=data.frame("X"=r_sample_binary(nsample, 
                        x = 0:1, prob = c(0.8,0.2), name = "Binary")))
simdata$data$event=as.numeric(simdata$data$failed)

If I look at the baseline survival probabilities from the simulated data:

> simdata$baseline$survivor
 [1] 0.9801987 0.9607894 0.9417645 0.9231163 0.9048374 0.8869204 0.8693582 0.8521438 0.8352702 0.8187308
[11] 0.8025188 0.7866279 0.7710516 0.7557837 0.7408182 0.7261490 0.7117703 0.6976763 0.6838614 0.6703200
[21] 0.6570468 0.6440364 0.6312836 0.6187834 0.6065307 0.5945205 0.5827483 0.5712091 0.5598984 0.5488116
[31] 0.5379444 0.5272924 0.5168513 0.5066170 0.4965853 0.4867523 0.4771139 0.4676664 0.4584060 0.4493290
[41] 0.4404317 0.4317105 0.4231621 0.4147829 0.4065697 0.3985190 0.3906278 0.3828929 0.3753111 0.3678794

and compare to my true survival probability, which can be calculated from my true baseline hazard:

St=exp(-lambda*(c(1:50)^rho))
> St
 [1] 0.9801987 0.9607894 0.9417645 0.9231163 0.9048374 0.8869204 0.8693582 0.8521438 0.8352702 0.8187308
[11] 0.8025188 0.7866279 0.7710516 0.7557837 0.7408182 0.7261490 0.7117703 0.6976763 0.6838614 0.6703200
[21] 0.6570468 0.6440364 0.6312836 0.6187834 0.6065307 0.5945205 0.5827483 0.5712091 0.5598984 0.5488116
[31] 0.5379444 0.5272924 0.5168513 0.5066170 0.4965853 0.4867523 0.4771139 0.4676664 0.4584060 0.4493290
[41] 0.4404317 0.4317105 0.4231621 0.4147829 0.4065697 0.3985190 0.3906278 0.3828929 0.3753111 0.3678794

The two are identical because I have large enough sample size. However if I calculate the KM estimators, they are quite different from St:

km_fit <- survfit(Surv(y, event) ~ X, data=simdata$data)
km=summary(km_fit, times = c(1:50))$surv[1:50]
> km
 [1] 0.9962515 0.9919738 0.9882571 0.9844797 0.9806792 0.9772240 0.9730233 0.9694796 0.9654927 0.9618635
[11] 0.9578942 0.9544436 0.9505571 0.9467933 0.9431748 0.9395630 0.9361808 0.9324463 0.9291364 0.9256540
[21] 0.9217372 0.9181445 0.9147239 0.9108842 0.9078288 0.9048339 0.9014076 0.8976344 0.8943120 0.8910748
[31] 0.8879073 0.8843805 0.8806476 0.8772176 0.8740362 0.8708713 0.8671436 0.8635299 0.8604316 0.8565315
[41] 0.8529732 0.8495042 0.8459996 0.8426619 0.8391398 0.8359729 0.8327033 0.8292976 0.8264344 0.0000000

Here's also a plot comparing the two:

Survival plot

Why are the two so different?

$\endgroup$
2
  • 1
    $\begingroup$ What are you putting in for the argument event to the Surv function when you haven't simulated any censoring process? Also on what value of X are you conditioning for the output of survfit? One may be using prediction at the means, while the other sets X to 0. $\endgroup$
    – AdamO
    Apr 2 at 15:47
  • 1
    $\begingroup$ If you do with (simdata$data,table(event, y==50)) you will see that all the times that reach 50 (ie, censored) have event==1, which suggests there's an issue with the coding of the event indicator. $\endgroup$ Apr 2 at 22:14
4
$\begingroup$

I am reasonably sure there's a problem with the simulation software. The documentation says

if censor.cond is FALSE then a proportion of the observations specified by censor is randomly and uniformly selected to be right-censored

That matches the code, which has

            ifelse(censor.cond, data$failed <- !censor.x(xdata, 
                censor = censor), data$failed <- !(runif(N) < 
                censor))

But that isn't what censoring should do. That sets the censoring indicator independent of the survival time, when actually it's the censoring hazard that should be independent of the survival time. Setting the censoring indicator independent of the survival time implies a lower censoring rate for long survival times and a higher one for short survival times.

On top of that, the censoring indicator appears to be coded the opposite way to how I would expect and how you used it

With censor=0.8

> with(simdata$data, table(event, y==50))
     
event FALSE  TRUE
    0 52629     0
    1 12818 34553

With censor=0.99

>  with(simdata$data, table(event, y==50))
     
event FALSE  TRUE
    0 64798     0
    1   649 34553

And if you look at the results of survfit as the censor parameter is increased, the fitted survival probabilities get higher and higher -- and further and further from your known baseline curve. That's what would happen with the documented censoring behaviour, but it's not what you'd usually want. The published paper about the method (Harden & Kropko) doesn't give much information about how censoring would be generated; it goes into detail on the survival time generation.

The code quoted above also makes me nervous. I looked up the R documentation and I think it's correct, since censor.cond is a single value and is always either TRUE or FALSE, so exactly one of the two assignments will be evaluated and the rules for lazy evaluation mean it will be evaluated in the calling frame of the ifelse() function rather than inside the function. But I think the authors wanted if rather than ifelse.

$\endgroup$
1
  • $\begingroup$ Thank you very much, Thomas, what you wrote makes a lot of senses to me. I actually did simulation again using my own written function, and the KM curve and true survival probability curve are quite consistent with each other. I'm going to write to the author reporting this issue. I'll post updates here if any. However I just don't get the part "the censoring indicator appears to be coded the opposite way" - I think it's alright as the failed variable is the failure indicator which is opposite to the censoring indicator. Please advise, thank you. $\endgroup$
    – Weiyu Qiu
    Apr 8 at 2:15
0
$\begingroup$

Thanks to everyone for the question and feedback. We are in the process of making some changes to address these issues. See the Github page for updates: https://github.com/jkropko/coxed.

On the censoring problem, we followed some previously published papers' simulations with that decision primarily because it was simple and censoring was not the main focus of our work. But we agree with Prof. Lumley that we need to generate a censoring hazard independent of survival time. We are working on that.

Thanks again!

$\endgroup$

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.