Kaplan-Meier and Nelson-Aalen estimates with and without survfit function from survival package

Data reading:


rm(list=ls(all=TRUE)) # wipe out everything else
require(survival)
data = aml[aml$x == "Maintained",]
data



   time status          x
1     9      1 Maintained
2    13      1 Maintained
3    13      0 Maintained
4    18      1 Maintained
5    23      1 Maintained
6    28      0 Maintained
7    31      1 Maintained
8    34      1 Maintained
9    45      0 Maintained
10   48      1 Maintained
11  161      0 Maintained

time <- data[,1]
status <- data[,2]

Function required to estimate Kaplan-Meier and Nelson-Aalen estimates without survfit function (almost as if the calculations were done by hand):

KM.tab <- function(time.var, status.var){

 data.temp <- cbind(time.var, status.var)

 ordered.time <- order(time.var)

 ordered.data <- data.temp[ordered.time,c(1:2)]

 unique.times <- !duplicated(ordered.data[,1])

 n.at.risk <- nrow(data.temp):1

 n.events <- rev(cumsum(rev(ordered.data[,2])))

 times <- ordered.data[unique.times,1]

 n.events <- rev(diff(rev(c(n.events[unique.times],0))))

 n.at.risk <- n.at.risk[unique.times]

 result.table <- as.data.frame(cbind(times, n.at.risk, n.events))

 result.table$N.A. <- cumsum(result.table$n.events/result.table$n.at.risk)

 result.table$N.A.Var <- cumsum(result.table$n.events/result.table$n.at.risk^2)

 result.table$Tsiatis.SE <- sqrt(result.table$N.A.Var)

 result.table$K.M. <- cumprod(1-result.table$n.events/result.table$n.at.risk)

 result.table$KM.Var <- cumsum(result.table$n.events/(result.table$n.at.risk*(result.table$n.at.risk+1.e-5-result.table$n.events)))*result.table$K.M.^2

 result.table$Grnwd.SE <- result.table$K.M.* sqrt(cumsum(result.table$n.events/(result.table$n.at.risk*(result.table$n.at.risk-result.table$n.events))))
 return(round(data.matrix(result.table[result.table$n.events>0,]),3))

}

Use the above function to get the estimates:

KM.tab(time, status)

Output

  times n.at.risk n.events  N.A. N.A.Var Tsiatis.SE  K.M. KM.Var Grnwd.SE
1     9        11        1 0.091   0.008      0.091 0.909  0.008    0.087
2    13        10        1 0.191   0.018      0.135 0.818  0.014    0.116
3    18         8        1 0.316   0.034      0.184 0.716  0.020    0.140
4    23         7        1 0.459   0.054      0.233 0.614  0.023    0.153
6    31         5        1 0.659   0.094      0.307 0.491  0.027    0.164
7    34         4        1 0.909   0.157      0.396 0.368  0.026    0.163
9    48         2        1 1.409   0.407      0.638 0.184  0.024    0.153

Use of survfit function from survival package

survfit.obj <- summary(survfit(Surv(time,status)~1))

survfit.obj

Call: survfit(formula = Surv(time, status) ~ 1)

 time n.risk n.event survival std.err lower 95% CI upper 95% CI

    9     11       1    0.909  0.0867       0.7541        1.000

   13     10       1    0.818  0.1163       0.6192        1.000

   18      8       1    0.716  0.1397       0.4884        1.000

   23      7       1    0.614  0.1526       0.3769        0.999

   31      5       1    0.491  0.1642       0.2549        0.946

   34      4       1    0.368  0.1627       0.1549        0.875

   48      2       1    0.184  0.1535       0.0359        0.944

round(survfit.obj$std.err[survfit.obj$n.ev>0],3)

0.087 0.116 0.140 0.153 0.164 0.163 0.153

Acknowledgement: original package, teaching link