6
$\begingroup$

This question comes from:

How to simulate a Cox proportional hazards model with change point and code it in R (See answer)

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?

$\endgroup$
5
$\begingroup$

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)
$\endgroup$
  • 1
    $\begingroup$ Thank you for your answer. I'll adapt the code to my simulation study in order to get some results. $\endgroup$ – anxoestevez Jul 21 '13 at 14:51

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.