tests/ekm.R

In survival: Survival Analysis

library(survival)
aeq <- function(x, y, ...) all.equal(as.vector(x), as.vector(y), ...)

# test for the "extended KM", where subjects change arms midstream
# (I don't like it statistically, but some use it).  

tdata <- aml
tdata$id <- 1:nrow(tdata)
tdata <- survSplit(Surv(time, status) ~ ., tdata, cut= c(9, 17, 30))
tdata$trt <- rep(c(1,1,2,2,2), length=nrow(tdata))
# different weights for different rows of the same subject = hardest case
tdata$wt <-  rep(1:6, length= nrow(tdata))
tdata$status[tdata$time==13] <- 1   # force at least 1 tied event

# not exported, but used in byhand
if (!exists("survflag")) survflag <- survival:::survflag

byhand <- function(t1, t2, status, grp, id, wt, debug=FALSE) {
    if (missing(wt)) wt <- rep(1, length(t1))
    ugrp <- unique(grp)
    ngrp <- length(ugrp)
    out <- vector("list", ngrp)
    names(out) <- ugrp
    pos <- survflag(Surv(t1, t2, status), id, grp)
    for (i in ugrp) { # create this curve
        keep <- (grp ==i)
        n <- sum(keep)
        utime <- sort(unique(c(t1[keep], t2[keep])))
        ntime <- length(utime)
        nrisk <- ncensor <- nevent <- entry <- double(ntime)
        surv <- cumhaz <- double(ntime)
        U <- C <- matrix(0, n, ntime) # influence
        U2 <- dN <- U  # portions of U, useful for debugging the AUC
        utemp <- ctemp <- double(n)
        km <- 1.0; chaz <- 0
        for (j in 1:ntime) {
            atrisk <- (keep & t1 < utime[j] &  t2 >= utime[j])
            nrisk[j] <- sum(wt[atrisk])
            nevent[j] <- sum(wt[keep & t2== utime[j]  & status ==1])
            ncensor[j] <- sum(wt[keep & t2==utime[j]  & status==0 & pos >1])
            entry[j] <- sum(wt[keep & t1== utime[j] & pos%%2 ==1])
            if (nrisk[j] >0) {
                km <- km * (nrisk[j]- nevent[j])/ nrisk[j]
                chaz <- chaz + nevent[j]/nrisk[j]
            }
            surv[j] =km
            cumhaz[j] =chaz

            # influence
            if (nrisk[j] > 0) {
                haz <- nevent[j]/nrisk[j]
                death <- (t2[keep]== utime[j]  & status[keep] ==1)
                temp <- double(n)
                temp[death] <- 1/nrisk[j]
                temp[atrisk[keep]] <- temp[atrisk[keep]] - haz/nrisk[j]
                ctemp <- ctemp + temp
                if (haz <1) utemp <- utemp - temp/(1-haz) else utemp <- 0
                dN[death,j] <- 1/(nrisk[j]* (1-haz))
                U2[atrisk[keep],j] <- haz/(nrisk[j] * (1-haz))
            }
            U[,j] <- utemp*km
            C[,j] <- ctemp
        }   
        out[[i]] <- list(n.id = length(unique(id[keep])), n= n,
                         time= utime, n.enter=entry, n.risk=nrisk, 
                         n.event=nevent, n.censor=  ncensor, surv= surv,
                         cumhaz = cumhaz, U=U, C=C, U2=U2, dN=dN)
        }
    out
}

true <- with(tdata, byhand(tstart, time, status, trt, id, wt))
ekm <- survfit(Surv(tstart, time, status) ~ trt, tdata, id=id, entry=TRUE,
               influence= TRUE, weights=wt)

aeq(ekm$n.id,    unlist(lapply(true, function(x) x$n.id)))
aeq(ekm$n,       unlist(lapply(true, function(x) x$n)))
aeq(ekm$time,    unlist(lapply(true, function(x) x$time)))
aeq(ekm$n.risk,  unlist(lapply(true, function(x) x$n.risk)))
aeq(ekm$n.enter, unlist(lapply(true, function(x) x$n.enter)))
aeq(ekm$n.event, unlist(lapply(true, function(x) x$n.event)))
aeq(ekm$n.censor,unlist(lapply(true, function(x) x$n.censor)))
aeq(ekm$surv    ,unlist(lapply(true, function(x) x$surv)))

# The byhand function gives per-observation influence, ekm has per-subject,
#  residuals can do either, but will fail with an error for this data when
#  collapse=TRUE with a message "same id appears in multiple curves "
rr <- residuals(ekm, times= c(9, 17, 30, 45), type="cumhaz")
aeq(rr[tdata$trt==1,], true[[1]]$C[, match(c(9,17,30,45), true[[1]]$time)])
aeq(rr[tdata$trt==2,], true[[2]]$C[, match(c(9,17,30,45), true[[2]]$time)])

rr <- residuals(ekm, times= c(9, 17, 30, 45), type= "pstate")
aeq(rr[tdata$trt==1,], true[[1]]$U[, match(c(9,17,30,45), true[[1]]$time)])
aeq(rr[tdata$trt==2,], true[[2]]$U[, match(c(9,17,30,45), true[[2]]$time)])

# Check influence returned by survfit
tdata1 <- subset(tdata, trt==1)
tdata2 <- subset(tdata, trt==2)
inf1 <- rowsum(tdata1$wt* true[[1]]$U, tdata1$id, reorder=FALSE)
inf2 <- rowsum(tdata2$wt* true[[2]]$U, tdata2$id, reorder=FALSE)
aeq(inf1, ekm$influence.surv[[1]])
aeq(inf2, ekm$influence.surv[[2]])

c1 <- rowsum(tdata1$wt* true[[1]]$C, tdata$id[tdata$trt==1], reorder=FALSE)
c2 <- rowsum(tdata2$wt* true[[2]]$C, tdata$id[tdata$trt==2], reorder=FALSE)
aeq(c1, ekm$influence.chaz[[1]])
aeq(c2, ekm$influence.chaz[[2]])

# Look at the AUC
t1 <- true[[1]]
width <- diff(c(t1$time, 50)) 
aucr <- width* t1$surv  # the rectangles that make up auc(50)
ainf <- t1$U %*% diag(width) # row i= influence of obs i on each rectangle
aucinf <- t(apply(ainf,1,cumsum))

rr <- residuals(ekm, times= c(t1$time[-1], 50), type= "auc")
aeq(rr[tdata$trt==1,], aucinf)

# use factor(status) to force the multi-state code
ekm2 <- survfit(Surv(tstart, time, factor(status)) ~ trt, tdata, id=id, 
                entry=TRUE, influence= TRUE, weights=wt)

aeq(ekm$n,        ekm2$n)
aeq(ekm$time,     ekm2$time)
aeq(ekm$n.risk,   ekm2$n.risk[,1])
aeq(ekm$n.event,  ekm2$n.event[,2])
aeq(ekm$n.censor, ekm2$n.censor[,1])
aeq(ekm$n.enter,  ekm2$n.enter[,1])
aeq(ekm$surv,     ekm2$pstate[,1])
aeq(ekm$std.err,  ekm2$std.err[,1])
aeq(ekm$cumhaz,   ekm2$cumhaz[,1])
aeq(ekm$std.chaz, ekm2$std.chaz[,1])
aeq(ekm$strata,   ekm2$strata)
aeq(ekm$n.id,     ekm2$n.id)
aeq(ekm$counts, 
    ekm2$counts[,c("nrisk:1", "ntrans:1.2", "ncensor:1", "nenter:1")])
aeq(ekm$influence.surv[[1]], ekm2$influence.pstate[[1]][,,1])
aeq(ekm$influence.surv[[2]], ekm2$influence.pstate[[2]][,,1])