R/lines.survfit.R

In survival: Survival Analysis

# This has a lot of common code with plot.survfit, make sure to keep them
#  in sync!
lines.survfit <- function(x, type='s', 
                          pch=3, col=1, lty=1, lwd=1,
                          cex=1,
                          mark.time=FALSE, xmax,
                          fun,  conf.int=FALSE,  
                          conf.times, conf.cap=.005, conf.offset=.012,
                          conf.type=c('log',  'log-log',  'plain', 
                                  'logit', "arcsin"),
                          mark, noplot="(s0)", cumhaz=FALSE, cumprob=FALSE, 
                          ...) {
    x <- survfit0(x, x$start.time)

    xlog <- par("xlog")
    # The default for plot and lines is to add confidence limits
    #  if there is only one curve
    if (!missing(conf.type) || is.null(x$conf.type)) 
        conf.type <- match.arg(conf.type) # check for legal
    else conf.type <- x$conf.type  # use the default in the curve
    if (conf.type=="none") conf.int <- FALSE

    if (missing(conf.int) && missing(conf.times))  
        conf.int <- (!is.null(x$std.err) && prod(dim(x) ==1))

    if (missing(conf.times)) conf.times <- NULL   
    else {
        if (!is.numeric(conf.times)) stop('conf.times must be numeric')
        if (missing(conf.int)) conf.int <- TRUE
    }
    if (!missing(conf.int)) {
        if (is.numeric(conf.int)) {
            conf.level <- conf.int
            if (conf.level<0 || conf.level > 1)
                stop("invalid value for conf.int")
            if (conf.level ==0) conf.int <- FALSE
            else if (conf.level != x$conf.int) {
                x$upper <- x$lower <- NULL  # force recomputation
            }
            conf.int <- TRUE
        }
        else conf.level = 0.95
    }

    # Organize data into stime, ssurv, supper, slower
    stime <- x$time
    std   <- NULL
    yzero <- FALSE   # a marker that we have an "ordinary survival curve" with min 0
    smat <- function(x) {
        # the rest of the routine is simpler if everything is a matrix
        dd <- dim(x)
        if (is.null(dd)) as.matrix(x)
        else if (length(dd) ==2) x
        else matrix(x, nrow=dd[1])
    }

    if (is.numeric(cumhaz)) { # plot the cumulative hazard
        if (!inherits(x, "survfitms") && any(cumhaz != 1))
            stop("numeric cumhaz argument only applies to multi-state")
        dd <- dim(x$cumhaz)
        if (is.null(dd)) nhazard <- 1
        else nhazard <- prod(dd[-1])

        if (!all(cumhaz == floor(cumhaz))) stop("cumhaz argument is not integer")
        if (any(cumhaz < 1 | cumhaz > nhazard)) stop("subscript out of range")
        ssurv <- smat(x$cumhaz)[,cumhaz, drop=FALSE]
        if (!is.null(x$std.chaz)) std <- smat(x$std.chaz)[,cumhaz, drop=FALSE]
        cumhaz <- TRUE # for the rest of the code
    } else if (cumhaz) {
        if (is.null(x$cumhaz)) 
            stop("survfit object does not contain a cumulative hazard")
        ssurv <- smat(x$cumhaz)
        if (!is.null(x$std.chaz)) std <- smat(x$std.chaz)
    }
    else if (inherits(x, "survfitms")) {
        if (!missing(cumprob) && !(is.logical(cumprob) && !cumprob)) {
            if (conf.int) 
                stop("confidence intervals not available when cumprob=TRUE")
            dd <- dim(x)
            j <- match("states", names(dd), nomatch=0)
            if (j==0) stop("survfitms object with no states dimension")
            
            #  cumprob is T/F or a vector of integers
            if (is.logical(cumprob)) cumprob <- 1:dd[j]
            else if (!is.numeric(cumprob) || any(cumprob <1 | cumprob > dd[j])
                     || any(cumprob != floor(cumprob)))
                stop("cumprob contains an invalid numeric")
                
            if (dd[j] ==1) {
                # nothing to do, user subscripted to only 1 state
                ssurv <- x$pstate
            } else {
                # reorder the states, pstate has dimension 2 or 3,
                #  time/strata is first, data (if present), then states
                #  (dd is the dimension from the user's point of view, of
                #  strata, data, state)
                if (length(dim(x$pstate))==2) {
                    # drop = FALSE for the rare case of a single time point
                    ssurv <- t(apply(x$pstate[,cumprob, drop=FALSE],1,cumsum))
                } else {
                    temp <- apply(x$pstate[,,cumprob, drop=FALSE],1:2, cumsum)
                    ssurv <- smat(aperm(temp, c(2,3,1)))
                }
                cumprob <- TRUE  # for the lastx line
            }
        } else {
            i <- !(x$states %in% noplot)
            if (all(i) || !any(i)) {
                # the !any is a failsafe, in case none are kept we ignore noplot
                ssurv <- smat(x$pstate)
                if (!is.null(x$std.err)) std <- smat(x$std.err)
                if (!is.null(x$lower)) {
                    slower <- smat(x$lower)
                    supper <- smat(x$upper)
                }
            }
            else {
                i <- which(i)  # the states to keep
                # we have to be careful about subscripting
                if (length(dim(x$pstate)) ==3) {
                    ssurv <- smat(x$pstate[,,i, drop=FALSE])
                    if (!is.null(x$std.err))
                        std <- smat(x$std.err[,,i, drop=FALSE])
                    if (!is.null(x$lower)) {
                        slower <- smat(x$lower[,,i, drop=FALSE])
                        supper <- smat(x$upper[,,i, drop=FALSE])
                    }
                }
                else {
                    ssurv <- x$pstate[,i, drop=FALSE]
                    if (!is.null(x$std.err)) std <- x$std.err[,i, drop=FALSE]
                    if (!is.null(x$lower)) {
                        slower <- smat(x$lower[,i, drop=FALSE])
                        supper <- smat(x$upper[,i, drop=FALSE])
                    }
                }
            }
        }
    }
    else {
        yzero <- TRUE
        ssurv <- as.matrix(x$surv)   # x$surv will have one column
        if (!is.null(x$std.err)) std <- as.matrix(x$std.err)
        # The fun argument usually applies to single state survfit objects
        #  First deal with the special case of fun='cumhaz', which is here for
        #  backwards compatability; people should use the cumhaz argument
        if (!missing(fun) && is.character(fun) && fun=="cumhaz") {
            cumhaz <- TRUE
            if (!is.null(x$cumhaz)) {
                ssurv <- as.matrix(x$cumhaz)
                if (!is.null(x$std.chaz)) std <- as.matrix(x$std.chaz)
            } 
            else {
                ssurv <- as.matrix(-log(x$surv))
                if (!is.null(x$std.err)) {
                    if (x$logse) std <- as.matrix(x$std.err)
                    else std <- as.matrix(x$std.err/x$surv)
                }
             }
        }
    }

    # set up strata
    if (is.null(x$strata)) {
        nstrat <- 1
        stemp <- rep(1, length(x$time)) # same length as stime
    }
    else {
        nstrat <- length(x$strata)
        stemp <- rep(1:nstrat, x$strata) # same length as stime
    }
    ncurve <- nstrat * ncol(ssurv)
    if (missing(conf.type)) {
        missingtype <- TRUE
        conf.type <- match.arg(conf.type)
    } else missingtype <- FALSE  # used below for cumhaz
    if (conf.type=="none") conf.int <- FALSE
    if (conf.int== "none") conf.int <- FALSE
    if (conf.int=="only") {
        plot.surv <- FALSE
        conf.int <- TRUE
        }
    else plot.surv <- TRUE

    if (conf.int) {
        if (is.null(std)) stop("object does not have standard errors, CI not possible")
        if (cumhaz) {
            if (missingtype) conf.type="plain"
            temp <- survfit_confint(ssurv, std, logse=FALSE,
                                    conf.type, conf.level, ulimit=FALSE)
            supper <- as.matrix(temp$upper)
            slower <- as.matrix(temp$lower)
        }
        else if (is.null(x$upper)) {
            if (missing(conf.type) && !is.null(x$conf.type))
                conf.type <- x$conf.type
            temp <- survfit_confint(ssurv, std, logse= x$logse,
                                    conf.type, conf.level, ulimit=FALSE)
            supper <- as.matrix(temp$upper)
            slower <- as.matrix(temp$lower)
        }
        else if (!inherits(x, "survfitms")) {
            supper <- as.matrix(x$upper)
            slower <- as.matrix(x$lower)
        }
    } else supper <- slower <- NULL
    if (!missing(fun)){
        if (is.character(fun)) {
            if (cumhaz) {
                tfun <- switch(tolower(fun),
                               'log' = function(x) x,
                               'cumhaz'=function(x) x,
                               'identity'= function(x) x,
                               stop("Invalid function argument")
                               )
            } else if (inherits(x, "survfitms")) {
                tfun <-switch(tolower(fun),
                              'log' = function(x) log(x),
                              'event'=function(x) x,
                              'cloglog'=function(x) log(-log(1-x)),
                              'cumhaz' = function(x) x,
                              'pct' = function(x) x*100,
                              'identity'= function(x) x,
                              stop("Invalid function argument")
                              )
            } else {
                yzero <- FALSE
                tfun <- switch(tolower(fun),
                           'log' = function(x) x,
                           'event'=function(x) 1-x,
                           'cumhaz'=function(x) x,
                           'cloglog'=function(x) log(-log(x)),
                           'pct' = function(x) x*100,
                           'logpct'= function(x) 100*x,  #special case further below
                           'identity'= function(x) x,
                           'f' = function(x) 1-x,
                           's' = function(x) x,
                           'surv' = function(x) x,
                           stop("Unrecognized function argument")
                           )
            }
        }
        else if (is.function(fun)) tfun <- fun
        else stop("Invalid 'fun' argument")
        
        ssurv <- tfun(ssurv )
        if (!is.null(supper)) {
            supper <- tfun(supper)
            slower <- tfun(slower)
        }
    }
    if (missing(mark.time) & !missing(mark)) mark.time <- TRUE
    if (missing(pch) && !missing(mark)) pch <- mark
    if (length(pch)==1 && is.character(pch)) pch <- strsplit(pch, "")[[1]]

    # Marks are not placed on confidence bands
    pch  <- rep(pch, length.out=ncurve)
    mcol <- rep(col, length.out=ncurve)
    if (is.numeric(mark.time)) mark.time <- sort(mark.time)

    # The actual number of curves is ncurve*3 if there are confidence bands,
    #  unless conf.times has been given.  Colors and line types in the latter
    #  match the curves
    # If the number of line types is 1 and lty is an integer, then use lty 
    #    for the curve and lty+1 for the CI
    # If the length(lty) <= length(ncurve), use the same color for curve and CI
    #   otherwise assume the user knows what they are about and has given a full
    #   vector of line types.
    # Colors and line widths work like line types, excluding the +1 rule.
    if (conf.int & is.null(conf.times)) {
        if (length(lty)==1 && is.numeric(lty))
            lty <- rep(c(lty, lty+1, lty+1), ncurve)
        else if (length(lty) <= ncurve)
            lty <- rep(rep(lty, each=3), length.out=(ncurve*3))
        else lty <- rep(lty, length.out= ncurve*3)
        
        if (length(col) <= ncurve) col <- rep(rep(col, each=3), length.out=3*ncurve)
        else col <- rep(col, length.out=3*ncurve)
        
        if (length(lwd) <= ncurve) lwd <- rep(rep(lwd, each=3), length.out=3*ncurve)
        else lwd <- rep(lwd, length.out=3*ncurve)
    }
    else {
        col  <- rep(col, length.out=ncurve)
        lty  <- rep(lty, length.out=ncurve)
        lwd  <- rep(lwd, length.out=ncurve)
    }

    # remember a prior xmax 
    if (missing(xmax)) xmax <- getOption("plot.survfit")$xmax 
    # Create a step function, removing redundancies that sometimes occur in
    #  curves with lots of censoring.
    dostep <- function(x,y) {
        keep <- is.finite(x) & is.finite(y) 
        if (!any(keep)) return()  #all points were infinite or NA
        if (!all(keep)) {
            # these won't plot anyway, so simplify (CI values are often NA)
            x <- x[keep]
            y <- y[keep]
        }
        n <- length(x)
        if (n==1)       list(x=x, y=y)
        else if (n==2)  list(x=x[c(1,2,2)], y=y[c(1,1,2)])
        else {
            # replace verbose horizonal sequences like
            # (1, .2), (1.4, .2), (1.8, .2), (2.3, .2), (2.9, .2), (3, .1)
            # with (1, .2), (.3, .2),(3, .1).  
            #  They are slow, and can smear the looks of the line type.
            temp <- rle(y)$lengths
            drops <- 1 + cumsum(temp[-length(temp)])  # points where the curve drops

            #create a step function
            if (n %in% drops) {  #the last point is a drop
                xrep <- c(x[1], rep(x[drops], each=2))
                yrep <- rep(y[c(1,drops)], c(rep(2, length(drops)), 1))
            }
            else {
                xrep <- c(x[1], rep(x[drops], each=2), x[n])
                yrep <- c(rep(y[c(1,drops)], each=2))
            }
            list(x=xrep, y=yrep)
        }
    }

    drawmark <- function(x, y, mark.time, censor, cex, ...) {
        if (!is.numeric(mark.time)) {
            xx <- x[censor>0]
            yy <- y[censor>0]
            if (any(censor >1)) {  # tied death and censor, put it on the midpoint
                j <- pmax(1, which(censor>1) -1)
                i <- censor[censor>0]
                yy[i>1] <- (yy[i>1] + y[j])/2
            }
        }
        else { #interpolate
            xx <- mark.time
            yy <- approx(x, y, xx, method="constant", f=0)$y
        }
        points(xx, yy, cex=cex, ...)
    }
    c1 <- 1  # keeps track of the curve number
    c2 <- 1  # keeps track of the lty, col, etc
    xend <- yend <- double(ncurve)
    if (length(conf.offset) ==1) 
        temp.offset <- (1:ncurve - (ncurve+1)/2)* conf.offset* diff(par("usr")[1:2])
    else temp.offset <- rep(conf.offset, length=ncurve) *  diff(par("usr")[1:2])
    temp.cap    <-  conf.cap    * diff(par("usr")[1:2])

    for (j in 1:ncol(ssurv)) {
        for (i in unique(stemp)) {  #for each strata
            who <- which(stemp==i)

            # if n.censor is missing, then assume any line that does not have an
            #   event would not be present but for censoring, so there must have
            #   been censoring then
            # otherwise categorize is 0= no censor, 1=censor, 2=censor and death
            if (is.null(x$n.censor)) censor <- ifelse(x$n.event[who]==0, 1, 0)
            else censor <- ifelse(x$n.censor[who]==0, 0, 1 + (x$n.event[who] > 0))
            xx <- stime[who]
            yy <- ssurv[who,j]
            if (conf.int) {
                ylower <- (slower[who,j])
                yupper <- (supper[who,j])
            }
            if (!is.null(xmax) && max(xx) > xmax) {  # truncate on the right
                xn <- min(which(xx > xmax))
                xx <- xx[1:xn]
                yy <- yy[1:xn]
                xx[xn] <- xmax
                yy[xn] <- yy[xn-1]
                if (conf.int) {
                    ylower <- ylower[1:xn]
                    yupper <- yupper[1:xn]
                    ylower[xn] <- ylower[xn-1]
                    yupper[xn] <- yupper[xn-1]
                }
            }
                

            if (plot.surv) {
                if (type=='s')
                    lines(dostep(xx, yy), lty=lty[c2], col=col[c2], lwd=lwd[c2]) 
                else lines(xx, yy, type=type, lty=lty[c2], col=col[c2], lwd=lwd[c2])
                if (is.numeric(mark.time) || mark.time) 
                    drawmark(xx, yy, mark.time, censor, pch=pch[c1], col=mcol[c1],
                             cex=cex)
            }
            xend[c1] <- max(xx)
            yend[c1] <- yy[length(yy)]

            if (conf.int && !is.null(conf.times)) {
                # add vertical bars at the specified times
                x2 <- conf.times + temp.offset[c1]
                templow <- approx(xx, ylower, x2,
                                  method='constant', f=1)$y
                temphigh<- approx(xx, yupper, x2,
                                  method='constant', f=1)$y
                segments(x2, templow, x2, temphigh,
                          lty=lty[c2], col=col[c2], lwd=lwd[c2])
                if (conf.cap>0) {
                    segments(x2-temp.cap, templow, x2+temp.cap, templow,
                             lty=lty[c2], col=col[c2], lwd=lwd[c2] )
                    segments(x2-temp.cap, temphigh, x2+temp.cap, temphigh,
                              lty=lty[c2], col=col[c2], lwd=lwd[c2])
                }
               
            }
            c1 <- c1 +1
            c2 <- c2 +1

            if (conf.int && is.null(conf.times)) {
                if (type == 's') {
                    lines(dostep(xx, ylower), lty=lty[c2], 
                          col=col[c2],lwd=lwd[c2])
                    c2 <- c2 +1
                    lines(dostep(xx, yupper), lty=lty[c2], 
                          col=col[c2], lwd= lwd[c2])
                    c2 <- c2 + 1
                }
                else {
                    lines(xx, ylower, lty=lty[c2], 
                          col=col[c2],lwd=lwd[c2], type=type) 
                    c2 <- c2 +1
                    lines(xx, yupper, lty=lty[c2], 
                          col=col[c2], lwd= lwd[c2], type= type)
                    c2 <- c2 + 1
                }
             }

        }
    }

    if (cumprob) {
        if (!is.null(xmax) && max(stime) > xmax) {  # truncate on the right
            keep <- (stime <= xmax)
            lastx <- list(x = stime[keep], y= ssurv[,keep])
        }
        else lastx <- list(x=stime, y=ssurv)
    }
    else lastx <- list(x=xend, y=yend)
    invisible(lastx)
}