# Get a desired percentage of censored observations in a simulation of Cox PH Model

This question comes from:

I want to generate a censoring variable $C=Exponential(\theta)$ that create a desired percentage of censorship on $Y$, like 33%, 55%, etc, by choosing an appropriate value for $\theta$.

In order to include this censoring variable in the code I did for the previous question, I added this to the dataframe generator:

cen <- rexp(n,theta) # Censoring variable.

ycen <- pmin(Y, cen) # Censoring effect.
di <- as.numeric(Y <= cen)


...

if (surv.df) data.frame(Y,X,ycen,di) else cbind(Y,X,ycen,di)


To do something like this (in the context of the aforementioned simulation):

coxph(Surv(ycen, di)~X)


Which $\theta$ value should I choose?

I am not sure if this answers your question but please find below some R code which follows Bender et al. (2005). They describe an approach to simulate a Cox PH regression model with given properties like the proportion of censored events (see line "dat <- data.frame(T = T, X, event = rbinom(n, 1, 0.30))", i.e. 70% of all events are censored).

Reference

Bender, Ralf, Thomas Augustin, und Maria Blettner. 2005. Generating survival times to simulate Cox proportional hazards models. Statistics in Medicine 24: 1713–1723.

##' Generate survival data with $p$ (correlated) predictors
##'
##'
##'
##' @title Generate survival data
##' @param n Sample size
##' @param beta Vector of coefficients
##' @param r Correlation between predictors
##' @param id.iter
##' @param id.study
##' @return matrix with identification variables id.iter and id.study,
##'         T (survival time), event (0: censored),
##'         predictors X1 to X$p$
##' @author Bernd Weiss
##' @references Bender et al. (2005)
genSurvData <- function(n = 100000,
beta = c(0.8, 2.2, -0.5, 1.1, -1.4),
r = 0.1,
id.iter = NA,
id.study = NA){

## Scale parameter (the smaller lambda, the greater becomes T)
lambda <- 0.000001#1.7

## Shape parameter
nue <- 8.9#9.4

## Sample size
n <- n

## Number of predictors
p <- length(beta)

## Generate column vector of coefficients
beta <- matrix(beta, ncol = 1)

## Generate correlated covariate vectors using a multivariate normal
## distribution with X ~ N(mu, S) and a given correlation matrix R, with:
## R: A p x p correlation matrix
## mu: Vector of means
## SD: Vector of standard deviations
## S: Variance-covariance matrix
R <- matrix(c(rep(r, p^2)), ncol = p)
diag(R) <- 1
R
mu <- rep(0, p)
SD <- rep(1, p)
S <- R * (SD %*% t(SD))
X <- mvrnorm(n, mu, S)
cov(X)
cor(X)
sqrt(diag(cov(X)))

## Calculate survival times
T <- (-log(runif(n)) / (lambda * exp(X %*% beta)))^(1/nue)

## 30% (0.30) of all marriages are getting divorced, i.e. 70% of all
## observations are censored ("event = rbinom(n, 1, 0.30)")
dat <- data.frame(T = T, X, event = rbinom(n, 1, 0.30))
## Also, all T's > 30 yrs are by definition censored and T is set to 30 yrs
dat$event <- ifelse(dat$T >= 30, 0, dat$event) dat$T <- ifelse(dat$T >= 30, 30, dat$T)

dat$id.iter <- id.iter dat$id.study <- id.study

## Reorder data frame: T, event, covariates
tmp.names <- names(dat)
dat <- dat[, c("id.iter", "id.study", "T", "event", tmp.names[grep("X", tmp.names)])]

## Returning a matrix speeds-up things a lot... lesson learned.
dat <- as.matrix(dat)
return(dat)
}

library(survival)
library(MASS)

dat <- genSurvData(n = 1000)
dat <- as.data.frame(dat)

survfit(Surv(time = T, event = event) ~ 1, data = dat)

coxph(Surv(time = T, event = event) ~ X1 + X2 + X3 +X4 +X5, data = dat)

• Thank you for your answer. I'll adapt the code to my simulation study in order to get some results. – anxoestevez Jul 21 '13 at 14:51