Sie sind hier: R code for chapter 2

# R code for chapter 2

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
}
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)
``` 