\section{Plotting survival curves}
The plot, lines, and points routines use several common code blocks in order to 
maintain consistency.

The xmax argument has been a long term issue.  Using xmax on a plot call, we
would like that xmax to persist in a subsequent lines.survfit call.  
But, the problem with this is that lines might not be called after plot.survfit:
someone might have other data and then want to add a survfit line to it (rare
case I know).  If we save the xlimits in some global object, there is no way
to erase that object every time a high level call is made.  

<<plot.survfit>>=
plot.survfit<- function(x, conf.int,  mark.time=FALSE,
			pch=3,  col=1,lty=1, lwd=1, 
                        cex=1, log=FALSE,
			xscale=1, yscale=1, 
			xlim, ylim, xmax, 
			fun, xlab="", ylab="", xaxs='r', 
                        conf.times, conf.cap=.005, conf.offset=.012, 
                        conf.type=c("log",  "log-log",  "plain", 
                                  "logit", "arcsin", "none"),
                        mark.col, noplot="(s0)", cumhaz=FALSE,
                        firstx, ymin, cumprob=FALSE, ...) {

    dotnames <- names(list(...))
    if (any(dotnames =='type'))
        stop("The graphical argument 'type' is not allowed")
    x <- survfit0(x, x$start.time)   # align data at 0 for plotting

    <<plot-log>>
    <<plot-data>>
    <<plot-confint>>
    <<plot-transform>>
    <<plot-setup-marks>>
    <<plot-makebox>> 
    <<plot-functions>>
    type <- 's'
    <<plot-draw>> 
    invisible(lastx)
}

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", "none"),
                          mark, noplot="(s0)", cumhaz=FALSE, cumprob=FALSE, 
                          ...) {
    x <- survfit0(x, x$start.time)

    xlog <- par("xlog")
    <<plot-data>>
    <<plot-confint>>
    <<plot-transform>>
    <<plot-setup-marks>>

    # remember a prior xmax 
    if (missing(xmax)) xmax <- getOption("plot.survfit")$xmax 
    <<plot-functions>>
    <<plot-draw>>
    invisible(lastx)
}

points.survfit <- function(x, fun, censor=FALSE,
                           col=1, pch, noplot="(s0)", cumhaz=FALSE, ...) {

    conf.int <- conf.times <- FALSE  # never draw these with 'points'
    cumprob <- FALSE; conf.type <- "none" 
    x <- survfit0(x, x$start.time)

    <<plot-data>>
    <<plot-transform>>
    
    if (ncurve==1 || (length(col)==1 && missing(pch))) {
        if (censor) points(stime, ssurv, ...)
        else points(stime[x$n.event>0], ssurv[x$n.event>0], ...)
    }
    else {
        c2 <- 1  #cycles through the colors and characters
        col <- rep(col, length=ncurve)
        if (!missing(pch)) {
            if (length(pch)==1)
                pch2 <- rep(strsplit(pch, '')[[1]], length=ncurve)
            else pch2 <- rep(pch, length=ncurve)
        }
        for (j in 1:ncol(ssurv)) {
            for (i in unique(stemp)) {
                if (censor) who <- which(stemp==i)
                else who <- which(stemp==i & x$n.event >0)
                if (missing(pch))
                    points(stime[who], ssurv[who,j], col=col[c2], ...)
                else
                    points(stime[who], ssurv[who,j], col=col[c2], 
                           pch=pch2[c2], ...) 
                c2 <- c2+1
            }
        }
    }
}
@ 

<<plot-log>>=
# decide on logarithmic axes, yes or no
if (is.logical(log)) {
    ylog <- log
    xlog <- FALSE
    if (ylog) logax <- 'y'
    else      logax <- ""
}
else {
    ylog <- (log=='y' || log=='xy')
    xlog <- (log=='x' || log=='xy')
    logax  <- log
}

if (!missing(fun)) {
    if (is.character(fun)) {
        if (fun=='log'|| fun=='logpct') ylog <- TRUE
        if (fun=='cloglog') {
            xlog <- TRUE
            if (ylog) logax <- 'xy'
            else logax <- 'x'
        }
        if (fun=="cumhaz" && missing(cumhaz)) cumhaz <- TRUE
    }
}
@ 
 
<<plot-data>>=
# 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)) {
    legal <- c("log",  "log-log",  "plain", 
                                  "logit", "arcsin", "none")
    if (!(conf.type %in% legal)) stop("conf.type must be one of ", 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 confidence limits are to be plotted, and they were not part of the
data that is passed in, create them.  Confidence limits for the 
cumulative hazard must always be created, and they don't use transforms.
<<plot-confint>>=
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=TRUE)
        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
@ 

The functional form of the fun argument can be whatever the user wants.
For the character form we try to thin out the obvious mistakes.
If fun=='cumhaz', the code above has already replaced ssurv with the
cumulative hazard, so this part of the code should plug in an identity
function.

<<plot-transform>>=
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)
    }
}
@ 
 
The \code{mark} argument is a holdover from S, when pch could not have
numeric values; mark has since disappeared from the manual page for
\code{par}.  We honor it for backwards compatability.
To be consistent with matplot and others, we allow pch to be a character
string or a vector of characters.

<<plot-setup-marks>>=
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)
}
@ 


Create the frame for the plot. 
We draw an empty figure, letting R figure out the limits.

<<plot-makebox>>=
# check consistency
if (!missing(xlim)) {
    if (!missing(xmax)) warning("cannot have both xlim and xmax arguments, xmax ignored")
    if (!missing(firstx)) stop("cannot have both xlim and firstx arguments")
}
if (!missing(ylim)) {
    if (!missing(ymin)) stop("cannot have both ylim and ymin arguments")
}

# Do axis range computations
if (!missing(xlim) && !is.null(xlim)) {
    tempx <- xlim
    xmax <- xlim[2]
    if (xaxs == 'S') tempx[2] <- tempx[1] + diff(tempx)*1.04
}
else {
    temp <-  stime[is.finite(stime)]
    if (!missing(xmax) && missing(xlim)) temp <- pmin(temp, xmax)
    else xmax <- NULL
    
    if (xaxs=='S') {
        rtemp <- range(temp)
        delta <- diff(rtemp)
        #special x- axis style for survival curves
        if (xlog) tempx <- c(min(rtemp[rtemp>0]), min(rtemp)+ delta*1.04)
        else tempx <- c(min(rtemp), min(rtemp)+ delta*1.04)
    }
    else if (xlog) tempx <- range(temp[temp > 0])
    else tempx <- range(temp)
}  
if (!missing(xlim) || !missing(xmax)) 
    options(plot.survfit = list(xmax=tempx[2]))
else options(plot.survfit = NULL)

if (!missing(ylim) && !is.null(ylim)) tempy <- ylim
else {
    skeep <- is.finite(stime) & stime >= tempx[1] & stime <= tempx[2]

    if (ylog) {
        if (!is.null(supper))
            tempy <- range(c(slower[is.finite(slower) & slower>0 & skeep], 
                             supper[is.finite(supper) & skeep]))
        else tempy <-  range(ssurv[is.finite(ssurv)& ssurv>0 & skeep])
        if (tempy[2]==1) tempy[2] <- .99   # makes for a prettier axis
        if (any(c(ssurv, slower)[skeep] ==0)) {
            tempy[1] <- tempy[1]*.8
            ssurv[ssurv==0] <- tempy[1]
            if (!is.null(slower))  slower[slower==0] <- tempy[1]
        }
    }
    else {
        if (!is.null(supper)) 
            tempy <- range(c(supper[skeep], slower[skeep]), finite=TRUE, na.rm=TRUE)
        else tempy <- range(ssurv[skeep], finite=TRUE, na.rm= TRUE)
        if (yzero) tempy <- range(c(0, tempy))
    }
}

if (!missing(ymin)) tempy[1] <- ymin

#
# Draw the basic box
#
temp <- if (xaxs=='S') 'i' else xaxs
plot(range(tempx, finite=TRUE, na.rm=TRUE)/xscale, 
     range(tempy, finite=TRUE, na.rm=TRUE)*yscale, 
     type='n', log=logax, xlab=xlab, ylab=ylab, xaxs=temp,...)
if(yscale != 1) {
    if (ylog) par(usr =par("usr") -c(0, 0, log10(yscale), log10(yscale))) 
    else par(usr =par("usr")/c(1, 1, yscale, yscale))   
}
if (xscale !=1) {
    if (xlog) par(usr =par("usr") -c(log10(xscale), log10(xscale), 0,0)) 
    else par(usr =par("usr")*c(xscale, xscale, 1, 1))   
}  
@ 
The use of [[par(usr)]] just above is a bit sneaky.  I want the
lines and points routines to be able to add to the plot, \emph{without}
passing them a global parameter that determines the y-scale or forcing
the user to repeat it.

The next functions do the actual drawing.
<<plot-functions>>=
# 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, ...)
}
@ 

The code to draw the lines and confidence bands.
<<plot-draw>>=
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)
@