I wish to simulate events from a Cox PH model where the censoring is informative, and to compare parameter estimator quality with estimates obtained from data generated by a Cox PH model with non-informative censoring. However I am struggling to do this, and I conclude rather than a technical flaw in my code I instead do not understand what is meant by informative and un-informative censoring.
My approach is to simulate an event time $T$ dependent on a vector of covariates $x$ having hazard function $h(t|x)=\lambda\exp(\beta'x)vt^{v-1}$. This corresponds to $T\sim Weibull(\lambda(x),v)$, where the scale parameter $\lambda(x)=\lambda\exp(\beta'x)$ depends on $x$ and the shape parameter $v$ is fixed. I have $N$ subjects where $T_{i}\sim Weibull(\lambda_{T}(x_{i}),v_{T})$, $\lambda_{T}(x_{i})=\lambda_{T}\exp(\beta_{T}'x_{i})$, for $i=1,...,N$. Here I assume the regression coefficients are $p$-dimensional.
I generate informative censoring times $C_{i}\sim Weibull(\lambda_{C}(x_{i}),v_{C})$, $\lambda_{C}(x_{i})=\lambda_{C}\exp(\beta_{C}'x_{i})$ and compute $Y^{I}_{i}=min(T_{i},C_{i})$ and a censored flag $\delta^{I}_{i}=1$ if $Y^{I}_{i}\leq C_{i}$ (an observed event), and $\delta^{I}_{i}=0$ if $Y^{I}_{i}> C_{i}$ (informatively censored: event not observed). I am convinced this is informative censoring because as long as $\beta_{T}\not=0$ and $\beta_{C}\not=0$ then for each subject the data generating process for $T$ and $C$ both depend on $x$.
In contrast I generate non-informative censoring times $D_{i}\sim Weibull(\lambda_{D}\exp(\beta_{D}),v_{D})$, and compute $Y^{NI}_{i}=min(T_{i},D_{i})$ and a censored flag $\delta^{NI}_{i}=1$ if $Y^{NI}_{i}\leq D_{i}$ (an observed event), and $\delta^{NI}_{i}=0$ if $Y^{NI}_{i}> D_{i}$ (non-informatively censored: event not observed). Here $\beta_{D}$ is a scalar.
I "scale" the simulation by choosing the $\lambda_{T}$, $\lambda_{C}$ and $\lambda_{D}$ parameters such that on average $T_{i}<C_{i}$ and $T_{i}<D_{i}$ to achieve $X\%$ of censored subjects for both $Y^{I}_{i}$ and $Y^{NI}_{i}$.
The problem is that even for say $30\%$ censoring (which I think is high), the Cox PH parameter estimates using both $Y^{I}_{i}$ and $Y^{NI}_{i}$ are unbiased when I expected the estimates using $Y^{I}_{i}$ to be biased, and I think I see why: however different $\beta_{C}$ is from $\beta_{T}$, a censored subject can presumably influence the estimation of $\beta_{T}$ only by affecting the set of subjects at risk at any time $t$, but this does not change the fact that every single $Y_{i}^{I}$ with $\delta_{i}^{I}=1$ will have been generated using $\beta_{T}$ only. Thus I do not see how my simulation can possibly produce biased estimates for $\beta_{T}$ using $Y_{i}^{I}$.
But then what is informative censoring if not based on this approach?
Edit: The code below I think demonstrates the point I refer to in my comments below this question, that is the bias in the estimation of the hazard function is not related to the bias in the estimation of the regression parameters. In my opinion (which may be wrong) the censoring in this code is not informative since the censoring does not depend on the covariates. If anybody can show me how to amend the code to produce informative censoring that produces biased regression estimates that would be great - then perhaps I will understand what informative censoring is.
##install.packages("muhaz")
library(muhaz)
censRand <- function(time, cens.t.5){ # cens.t.5 is the t 1/2 of censor process
ctime <- rexp(n = length(time), rate = 1/cens.t.5)
event <- (time <= ctime)
t_obs <- pmin(time, ctime)
return(data.frame(Times=t_obs, event=event))
}
N=10000
#simulate a hazard of the form;
# h(t|x) = lambda*exp[(beta_age*age)+(beta_treat*treat)]
# Do this by drawing from an exponential distribution
beta_age = 0.5
beta_treat = -1
age = rnorm(N,33,2)
treat = rbinom(N,1,0.5)
lp = (beta_age*age)+(beta_treat*treat)
lambda=0.000001
time=rexp(N, lambda*exp(lp))
##censoring times
ctime <- censRand( time, 0.1)
pc_cen = round((length(ctime[ctime[,2]==FALSE,2])/N)*100,2)
##plots
par(mfrow=c(2,2))
tit = paste("CTTE (censored = ",toString(pc_cen)," %)")
hist(time,breaks=20,main="TTE")
hist(ctime[,1],breaks=20,main=tit)
##fit the Cox model
cox_obj = coxph(Surv(ctime[ ,1],ctime[ ,2],type=c("right"))~age+treat )
summary(cox_obj)
plot( muhaz( ctime[ ,1], ctime[,2]) )
edit. $\endgroup$