Zu den Inhalten springen

Institute of Medical Biometry and Statistics (IMBI)

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)