[R] Kaplan-Meier for left truncated data?

[gb]

On Sun, 21 May 2000, gb wrote:

> 
> Is there a function somewhere for estimating the survivor function
> (plus confidence limits) when data are both left truncated and right
> censored?  survfit in  survival5  doesn't like left truncation.

I got one answer, but a function for left censoring (for which I am 
grateful). So I wrote my own function. Here it is. Be careful, there
are a few nasty loops... and no serious testing yet.

Tips on improvements are welcome! (Note, it just plots, no return
values. Easy to change, though.)

Göran
------------------------------------------------------------------------
plot.surv <- function(enter = rep( 0, length(exit) ),
                      exit,
                      event = rep( 1, length(exit) ),
                      group = rep( 1, length(exit) ),
                      limits = F,
                      conf = 0.95)
  {
    ## Input data:
    ##
    ## enter : left truncation point
    ## exit  : exit time point
    ## event : if zero, a censored observation; otherwise an event.
    ## group : one curve for each distinct value in group.
    ## limits: if TRUE, and only one group, pointwise confidence
    ##         limits (Greenwoods formula, log(-log) type).
    ## conf  : confidence level. Can be given as a percentage.
    
    ## Check input data:
    n <- length(exit)
    if (length(enter) != n)stop("enter and exit must have equal length.")
    if (length(event) != n)
      stop("event and exit must have equal length.")
    if (length(group) != n)
      stop("group and exit must have equal length.")
    if (min(exit - enter) <= 0) stop("Interval length must be positive.")
    if (conf > 1) conf <- conf / 100 ## conf given as a percentage(?)
    if ( (conf < 0.5) | (conf >=1) ) stop("Bad conf value")

    grupp <- as.character(group)
   
    strata <- sort(unique(grupp))
    if (length(strata) > 9)
      stop("Too many groups. No more than 9 are allowed.")
    
    ##
    if (length(strata) > 1) limits <- F # No limits if multiple curves.
    gang <- 0
    for (stratum in strata)
      {
        atom <- table.events(enter[grupp == stratum],
                             exit[grupp == stratum],
                             event[grupp == stratum])
        
        gang <- gang + 1
        surv <- c( 1, cumprod(1 - atom$events / atom$riskset.sizes) )
        if (gang == 1)
          {
            plot(c(0, atom$times), surv, type = "s",
                 xlab = "duration", ylab = "Surviving fraction",
                 main = "Survivor function(s)", ylim = c(0, 1),
                 col = gang%%5)
            if (limits) ## Greenwood's formula,
                        ## Kalbfleisch & Prentice, p. 15 (note error!).
              {
                q.alpha <- abs(qnorm((1 - conf) / 2))
                survived <- (atom$riskset.size - atom$events)
                se <- sqrt(cumsum(atom$events /
                                   ( atom$riskset.sizes * survived )
                                   )
                           )/
                            cumsum(-log(survived / atom$riskset.sizes))
                upper <- surv ^ exp(q.alpha * c(0, se))
                lower <- surv ^ exp(-q.alpha * c(0, se))
                lines(c(0, atom$times), upper, type = "s",
                      col = gang%%5 + 1)
                lines(c(0, atom$times), lower, type = "s",
                      col = gang%%5 + 1)
              }
          }
        else
          {
            lines(c(0, atom$times), surv, type = "s", 
                  col = gang%%5)
          }
      }
    abline(h=0)
    if (length(strata) > 1)
      {
        ## legend.txt <- as.character(strata)
        colors <- (1:length(strata))%%5
        legend(0, 0.6, lty = 1, legend = strata, col = colors)
      }
  }

table.events <- function(enter = rep(0, length(exit)),
                         exit,
                         event)
{
  n <- length(exit)

  ## Check input data:
  if ( length(enter) != n ) stop("enter and exit must have equal length.")
  if ( length(event) != n ) stop("event and exit must have equal length.")
  ##
  
  event <- (event != 0) ## 0 (F) = censoring, else (T) = event

  times <- c(unique(sort(exit[event])))
  nn <- length(times)

  rs.size <- double(nn)
  n.events <- double(nn)

  for (i in 1:nn) ## Try to avoid this loop!
    {
      rs.size[i] <- sum( (enter < times[i]) &
                        (exit >= times[i]) )
      n.events[i] <- sum( (exit == times[i]) & event )
    }

  stop.at <- which(rs.size == n.events)
  if (length(stop.at))
    {
      stop.at <- min(stop.at) - 1
      if (stop.at <= 0) stop("Bad data. All died immediately!")
      times <- times[1:stop.at]
      n.events <- n.events[1:stop.at]
      rs.size <- rs.size[1:stop.at]
    }
      
  return ( list(times         = times,
                events        = n.events,
                riskset.sizes = rs.size)
          )
}

Have fun! And send bug reports to
-- 
 Göran Broström                    e-mail: gb at stat.umu.se
 Professor                           tel: +46 90 786 5223
 Department of Statistics            fax: +46 90 786 6614
 Umeå University
 SE-90187 Umeå, Sweden             http://www.stat.umu.se/egna/gb

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._