| The prodint() function used throughout this chapter for product integration requires two arguments: a vector of time points and a function. |
- prodint <- function(time.points, A) { times <- c(0, sort(unique(time.points))) S <- prod(1 - diff(apply(X = matrix(times), MARGIN = 1, FUN = A))) return(S) }
|
| A.exp is a function which computes the cumulative hazard with hazard(t) = 0.9 using a vector of time points that you have to enter. |
- A.exp <- function(time.point) { return(0.9 * time.point) }
|
| The first example (p. 12) uses a sequence of timepoints between 0 and 1. |
- times <- seq(0, 1, 0.001) prodint(times, A.exp) exp(-0.9 * max(times))
|
| The next example uses a vector of uniformly distributed random numbers. |
- prodint(runif(n = 1000, min = 0, max = 1), A.exp)
|
| Creation of a similar function for Weibull distributed survival times. |
- A.weibull <- function(time.point){ return(2 * time.point^0.25) prodint(times, A.weibull) ## True result exp(-2 * max(times)^0.25) } ## finer grid of time points prodint(seq(0, 1, 0.000001), A.weibull)
|
| Simulate some event times with constant hazard equal to 0.9. |
- event.times <- rexp(100,0.9)
|
| Computing survival function S(t) with survfit function, which is part of the survival package. |
- library(survival) fit.surv <- survfit(Surv(event.times) ~ 1)
|
| Compute the Nelson-Aalen estimator based on what survfit returns. |
- A <- function(time.point) { sum(fit.surv$n.event[fit.surv$time <= time.point]/ fit.surv$n.risk[fit.surv$time <= time.point]) }
|
| One can see that using product integration leads to the same result. |
- prodint(event.times[event.times <= 1], A)
|
| Computation of the empirical survival function: |
- sum(event.times > 1) / length(event.times)
|
| Now, we simulate censoring times following a uniform distribution. |
- cens.times <- runif(100,0,5)
|
| Next, we take the minimum of the event times and the censoring times. |
- obs.times <- pmin(event.times, cens.times)
|
| and with their help we create an event indicator |
- event.times <= cens.times
|
| Survival probability estimate: |
- fit.surv <- survfit(Surv(obs.times, event.times <= cens.times) ~ 1)
|
| with product integration |
- prodint(obs.times[obs.times<=1], A)
|
| estimate at time 1 using survfit |
- S <- fit.surv$surv S[fit.surv$time <= 1][length(S[fit.surv$time <= 1])]
|
| Figure 2.3: nt counting process |
- nt <- cumsum(fit.surv$n.event)
|
| Function that computes the compensator |
- integral2 <- function(x) { tmp <- rep(0,length(x)) for(i in 1:length(x)){ if(x[i]<=min(fit.surv$time)){ tmp[i] <- x[i] * fit.surv$n.risk[1] } else{ if(x[i]>=max(fit.surv$time)){ tmp[i] <- sum(diff(c(0,fit.surv$time)) * fit.surv$n.risk) } else{ bla <- diff(c(0,fit.surv$time)[c(0,fit.surv$time) < x[i]]) tmp[i] <- sum(bla * fit.surv$n.risk[1:length(bla)]) + (x[i]-max(fit.surv$time[fit.surv$time<x[i]])) * fit.surv$n.risk[length(bla)+1] } } } return(0.9 * tmp) }
|
| Plot the compensator and counting process. |
- plot(x = c(0, 4), y = c(0, 80), xlab = expression(paste(Time, " ", italic(t))), ylab = "", type="n", axes = F, main = "", cex.main = 1.7, cex.lab = 1.5) axis(1, at = seq(0, 4, 1), cex.axis = 1.5) axis(2, at = seq(0, 80, 10), cex.axis = 1.25) box() lines(fit.surv$time, nt, type="s", lwd=2) curve(integral2, from = 0, to = 5, add = TRUE, lwd=2, n = 1001)
|
