Sie sind hier: Startseite Competing Risks and … R Code R code for chapter 3

R code for chapter 3

Section 3.2

First simulation example: with constant cause-specifc hazards.

set.seed(1383)
event.times <- rexp(100, 0.9)
f.cause <- rbinom(100, size = 1, prob = 1/3)
f.cause <- ifelse(f.cause == 0, 2, 1)
cens.times <- runif(100,0,5)
## truncation times
lt.times <- rgamma(100, shape= 0.5, rate = 2) 

obs.times <- pmin(event.times, cens.times)
obs.cause <- c(event.times <= cens.times) * f.cause

table(obs.cause)

Second example: α02(t) Weibull. Here, we use the inversion method to generate the event times.

## Drawing of uniformly distributed random variables
my.times <- runif(100)
## inversion method
my.times <- -0.5 + (0.25 - log(1 - my.times))^0.5

## function that computes the cause of failure at time x
cause1 <- function(x) {
    out <- rbinom(length(x), 1, prob = 1 / (1 + 2 * x))
    ifelse(out == 0, 2, 1)
}

my.cause <- cause1(my.times)

For illustration, we create a function to generate event times using uniroot(). It takes as arguments n, the sample size and max.int which means that the solution will be searched between 0 and max.int. As output you will get the event times.

generate.my.times.v2 <- function(n, max.int) {
    temp <- function(x,y) { return(x + x^2 + y) }
    stime <- NULL
    i <- 1
    
    while(length(stime) < n) {                           
        u <- runif(1)
        ## If endpoints are of opposite sign, call uniroot:
        if (temp(0, log(1 - u)) * temp(max.int, log(1 - u)) 
            < 0) {
            res <- uniroot(temp, c(0, max.int),
                           tol = 0.0001, y = log(1 - u))                        
            stime[i] <- res$root
            i <- i + 1
        }
        else cat("Values at endpoints not 
                  of opposite signs. \n")
    }
    return(stime)
}

Example:

my.times <- generate.my.times.v2(n=100, max.int=99)