# Generating censoring times for the cox proportional hazards model

I am trying to understand the different ways of simulating survival times in a cox proportional hazards model. A simple example consists in simulating the event times following a Weibull distribution. I refer to this answer (https://stats.stackexchange.com/q/135129) for more deatils on the simulation procedure:

I paste part of the answer that details the simulation of event times in the Weibull case below :

Example [Weibull baseline hazard]

Let $$h_0(t) = \lambda \rho t^{\rho - 1}$$ with shape $$\rho > 0$$ and scale $$\lambda > 0$$. Then $$H_0(t) = \lambda t^\rho$$ and $$H^{-1}_0(t) = (\frac{t}{\lambda})^{\frac{1}{\rho}}$$. Following the inverse probability method, a realisation of $$T \sim S(\cdot \,|\, \mathbf{x})$$ is obtained by computing $$t = \left( - \frac{\log(v)}{\lambda \exp(\mathbf{x}^\prime \mathbf{\beta})} \right)^{\frac{1}{\rho}}$$ with $$v$$ a uniform variate on $$(0, 1)$$. Using results on transformations of random variables, one may notice that $$T$$ has a conditional Weibull distribution (given $$\mathbf{x}$$) with shape $$\rho$$ and scale $$\lambda \exp(\mathbf{x}^\prime \mathbf{\beta})$$.

Right censoring of the event times implies that we observe min($$t$$,$$C$$) where $$C$$ is the censoring variable. I try out two ways of generating the censoring times below (20 percent censoring).

I also start from the same code in the answer with a different way of dealing with the censoring mechanism. I only changed the line dealing with generation of the censoring variable C in the code. Instead of simulating following an exponential distribution (as in the answer cited), I randomly generate 0's and 1's to censor the event time variable. :

simulWeib <- function(N, lambda, rho, beta, rateC, method) {
# covariate --> N Bernoulli trials
x <- sample(x=c(0, 1), size=N, replace=TRUE, prob=c(0.5, 0.5))

# Weibull latent event times
v <- runif(n=N)
Tlat <- (- log(v) / (lambda * exp(x * beta)))^(1 / rho)

# censoring times
if (method == "method 1") {
C = c(rep(0,N*rateC),rep(1,N*(1-rateC)))
C <- sample(C)

# follow-up times and event indicators
time <- Tlat
status <- C
} else if (method == "method 2") {
# censoring times
C <- rexp(n=N, rate=rateC)

# follow-up times and event indicators
time <- pmin(Tlat, C)
status <- as.numeric(Tlat <= C)
}

# data set
data.frame(id=1:N,
time=time,
status=status,
x=x)
}

set.seed(1234)
betaHat_1 <- rep(NA, 1e3)
betaHat_2 <- rep(NA, 1e3)
for(k in 1:1e3)
{
dat1 <- simulWeib(N=100, lambda=0.01, rho=1, beta=-0.6, rateC=0.2, method = "method 1")
fit1 <- coxph(Surv(time, status) ~ x , data=dat1)
betaHat_1[k] <- fit1$coef dat2 <- simulWeib(N=100, lambda=0.01, rho=1, beta=-0.6, rateC=0.001, method = "method 2") fit2 <- coxph(Surv(time, status) ~ x , data=dat2) betaHat_2[k] <- fit2$coef
}

mean(betaHat_1)
> -0.6169394
mean(betaHat_2)
> -0.613186


With method 1, rateC gives directly the censoring rate but I am not sure this could still be considered right censoring. The common way of generating censored times is by method 2. In this example, the two estimates are similar though. Method 2 can be considered right censoring but what about method 1 ?

If the question or part of the code is not clear, feel free to ask for edits.

Thank you !

EDIT : I am mainly interested in random right censoring. If we look at "method 2" in the code, this is the classic way of simulating right censoring times. That is, you simulate event times using the inverse probability method as shown in the answer. Then the censoring times are simulated following an exponential distribtuion with a rate parameter that is adjusted so that we get the desired censoring level when we compute the right censored time : right-censored time = min(event_time,censoring_time). This part corresponds to "method 2" in the code and is clear.

When I was less well versed in the subject of survival analysis (3 years ago), I naively censored the event times by generating a censor variable with 0's (censored) and 1's (not censored) and randomly allocated the 0's and 1's to the event times. This corresponds to "method 1" in the code. What type of censoring does this correspond to ? Is this even considered censoring ? These are my questions.

## 2 Answers

Ok, so your aim is to create/generate censoring times for the cox proportional hazards model, and you have already tried the solutions proposed here and here and they don´t answer your question.

I would have recommended to look closer into properties of the data you are trying to model, but it became clear from your comments that you have no interest in real data or doing any predictive modeling, but rather trying to understand the different ways of censoring data.

About censoring.

Censoring has been often categorized into type I and type II censoring. Also, in real life, you come across censored data quite a lot as you get the data from measurements, which often have some cut-off value and there is limit-of-detection censoring. All these types are different from random cencoring.

Type I. Imagine a trial (or a study), where participants are followed for a prespecified, fixed duration of time. Than the event can get any value smaller than this fixed duration - thus time-to-event is right-censored, and this is type I censoring.

Type II. In type II censoring, the process/study/trial is stopped when a prespecified number of participants have experienced the event of interest.

Limit-of-detection censoring. This can be both - type I or type II, depending on the parameters.

Randomly-censored data is often considered as randomly-missing data, and is treated as such.

What type of censoring you are interested in?

• not really. I refer to that answer in my question as the code I use as example is taken from that answer. My question is about the censoring time in particular. I illustrate this using two methods and from the estimates obtained, it seems the methods are equivalent. However in a more slightly more complex setting I am working on, I get biased estimates when using method 1. I am struggling to understand why though. Apr 23, 2020 at 13:40
• and what about this: stats.stackexchange.com/questions/267481/…? ?
– Oka
Apr 23, 2020 at 15:40
• this doesn't answer my question. I have no interest in doing any predictive modeling. I am just trying to understand the different ways of censoring data and what the method 1 I implemented implies Apr 23, 2020 at 20:16
• I added an EDIT part to the question hoping that my question is clearer now. Apr 24, 2020 at 9:44

In theory you should have an independent random variable (any) for the time you "observe", and event status is then 1 if you observed long enough till the event time, or 0 if not. (So method 2, but you may use any other distribution, not just exponential, just independent from event times; can be a constant).

In the 1st method you generated event times with a given parameters and distribution, and then said that it is actually not the event time for some observations (by assigning censored = 0), but event actually happened sometime later (if I understood the code correctly). So this creates a non-random censoring as you artificially made some event times look longer than they were. For Cox model and for the estimate of the HR this may not be a problem a those non-censored would still exhibit dependencies as modelled. However I guess if you use survreg() this may give different estimate, and baseline hazard in both case may be skewed to the one with lower hazard.