# Shameless copy of the original plot.survfit funciton from Terry Therneau with
# one modification to reproduce the behaviour of SAS in plotting censored event(s)
# event when they co-occur with real event(s)
# Change that: (x$n.event[who] ==0 & x$n.censor[who] >0) # censoring ties
# to that: (x$n.event[who] ==0 | x$n.censor[who] >0) #NOT censoring ties
.sibylKMplot <- function(x, conf.int, mark.time=FALSE,
mark=3, col=1,lty=1, lwd=1,
cex=1, log=FALSE,
xscale=1, yscale=1,
firstx=0, firsty=1,
xmax, ymin=0,
fun, xlab="", ylab="", xaxs='S',
conf.times, conf.cap=.005, conf.offset=.012, ...) {
dotnames <- names(list(...))
if (any(dotnames=='type'))
stop("The graphical argument 'type' is not allowed")
if (missing(mark.time) & !missing(mark)) mark.time <- TRUE
if (inherits(x, "survfitms")) {
x$surv <- 1- x$pstate
if (is.matrix(x$surv)) {
dimnames(x$surv) <- list(NULL, x$states)
if (ncol(x$surv) > 1 && any(x$states == '')) {
x$surv <- x$surv[, x$states != '']
if (is.matrix(x$p0)) x$p0 <- x$p0[, x$states != '']
else x$p0 <- x$p0[x$states != '']
}
}
if (!is.null(x$lower)) {
x$lower <- 1- x$lower
x$upper <- 1- x$upper
}
if (missing(fun)) fun <- "event"
}
if (missing(firsty) && !is.null(x$p0)) firsty <- 1-x$p0
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'
}
}
}
# The special x axis style only applies when firstx is not given
if (missing(xaxs) && (firstx!=0 || !missing(fun) ||
(missing(fun) && inherits(x, "survfitms"))))
xaxs <- par("xaxs") #use the default
ssurv <- as.matrix(x$surv)
stime <- x$time
if( !is.null(x$upper)) {
supper <- as.matrix(x$upper)
slower <- as.matrix(x$lower)
}
else {
conf.int <- FALSE
supper <- NULL #marker for later code
}
# 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)
firsty <- matrix(firsty, nrow=nstrat, ncol=ncol(ssurv))
if (!missing(xmax) && any(x$time>xmax)) {
# prune back the survival curves
# I need to replace x's over the limit with xmax, and y's over the
# limit with either the prior y value or firsty
keepx <- keepy <- NULL # lines to keep
tempn <- table(stemp)
offset <- cumsum(c(0, tempn))
for (i in 1:nstrat) {
ttime <-stime[stemp==i]
if (all(ttime <= xmax)) {
keepx <- c(keepx, 1:tempn[i] + offset[i])
keepy <- c(keepy, 1:tempn[i] + offset[i])
}
else {
bad <- min((1:tempn[i])[ttime>xmax])
if (bad==1) { #lost them all
if (!is.na(firstx)) { # and we are plotting lines
keepy <- c(keepy, 1+offset[i])
ssurv[1+offset[i],] <- firsty[i,]
}
}
else keepy<- c(keepy, c(1:(bad-1), bad-1) + offset[i])
keepx <- c(keepx, (1:bad)+offset[i])
stime[bad+offset[i]] <- xmax
x$n.event[bad+offset[i]] <- 1 #don't plot a tick mark
}
}
# ok, now actually prune it
stime <- stime[keepx]
stemp <- stemp[keepx]
x$n.event <- x$n.event[keepx]
if (!is.null(x$n.censor)) x$n.censor <- x$n.censor[keepx]
ssurv <- ssurv[keepy,,drop=FALSE]
if (!is.null(supper)) {
supper <- supper[keepy,,drop=FALSE]
slower <- slower[keepy,,drop=FALSE]
}
}
#stime <- stime/xscale #scaling is deferred until xmax processing is done
if (!missing(fun)) {
if (is.character(fun)) {
tfun <- switch(fun,
'log' = function(x) x,
'event'=function(x) 1-x,
'cumhaz'=function(x) -log(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,
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)
}
firsty <- tfun(firsty)
}
if (missing(firstx)) {
if (!is.null(x$start.time))
firstx <- x$start.time
else {
if (xlog) firstx <- min(stime[stime>0])
else firstx <- min(0, stime)
}
}
# The default for plot and lines is to add confidence limits
# if there is only one curve
if (missing(conf.int) && missing(conf.times)) conf.int <- (ncurve==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 (is.logical(conf.int)) plot.surv <- TRUE
else {
temp <- match.arg(conf.int, c("both", "only", "none"))
if (is.na(temp)) stop("invalid value for conf.int")
if (temp=="none") conf.int <- FALSE else conf.int <- TRUE
if (temp=="only") plot.surv <- FALSE else plot.surv <- TRUE
}
# Marks are not placed on confidence bands
mark <- rep(mark, 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)
}
#axis setting parmaters that depend on the fun argument
if (!missing(fun)) {
ymin <- tfun(ymin) #lines routine doesn't have it
}
# Do axis range computations
if (xaxs=='S') {
#special x- axis style for survival curves
xaxs <- 'i' #what S thinks
tempx <- max(stime) * 1.04
}
else tempx <- max(stime)
tempx <- c(firstx, tempx, firstx)
if (ylog) {
tempy <- range(ssurv[is.finite(ssurv)& ssurv>0])
if (tempy[2]==1) tempy[2] <- .99
if (any(ssurv==0)) {
tempy[1] <- tempy[1]*.8
ssurv[ssurv==0] <- tempy[1]
if (!is.null(supper)) {
supper[supper==0] <- tempy[1]
slower[slower==0] <- tempy[1]
}
}
tempy <- c(tempy, firsty)
}
else tempy <- range(ssurv, firsty, finite=TRUE, na.rm=TRUE)
if (missing(fun)) {
tempx <- c(tempx, firstx)
if (!ylog) tempy <- c(tempy, ymin)
}
#
# Draw the basic box
#
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=xaxs,...)
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))
}
# 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]
yy <- y[censor]
}
else { #interpolate
xx <- mark.time
yy <- approx(x, y, xx, method="constant", f=0)$y
}
points(xx, yy, cex=cex, ...)
}
plot.surv <- TRUE
type <- 's'
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)
censor <- if (is.null(x$n.censor))
(x$n.event[who] ==0) else
(x$n.event[who] ==0 & x$n.censor[who] >0) #censoring ties
xx <- c(firstx, stime[who])
censor <- c(FALSE, censor) #no mark at firstx
yy <- c(firsty[i,j], ssurv[who,j])
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=mark[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, c(firsty[i,j], slower[who,j]), x2,
method='constant', f=1)$y
temphigh<- approx(xx, c(firsty[i,j], supper[who,j]), 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, c(firsty[i,j], slower[who,j])), lty=lty[c2],
col=col[c2],lwd=lwd[c2])
c2 <- c2 +1
lines(dostep(xx, c(firsty[i,j], supper[who,j])), lty=lty[c2],
col=col[c2], lwd= lwd[c2])
c2 <- c2 + 1
}
else {
lines(xx, c(firsty[i,j], slower[who,j]), lty=lty[c2],
col=col[c2],lwd=lwd[c2], type=type)
c2 <- c2 +1
lines(xx, c(firsty[i,j], supper[who,j]), lty=lty[c2],
col=col[c2], lwd= lwd[c2], type= type)
c2 <- c2 + 1
}
}
}
}
invisible(list(x=xend, y=yend))
}
lines.survfit <- function(x, type='s',
mark=3, col=1, lty=1, lwd=1,
cex=1,
mark.time=FALSE, xscale=1,
firstx=0, firsty=1, xmax,
fun, conf.int=FALSE,
conf.times, conf.cap=.005, conf.offset=.012, ...) {
xlog <- par("xlog")
if (missing(mark.time) & !missing(mark)) mark.time <- TRUE
if (inherits(x, "survfitms")) {
x$surv <- 1- x$pstate
if (is.matrix(x$surv)) {
dimnames(x$surv) <- list(NULL, x$states)
if (ncol(x$surv) > 1 && any(x$states == '')) {
x$surv <- x$surv[, x$states != '']
if (is.matrix(x$p0)) x$p0 <- x$p0[, x$states != '']
else x$p0 <- x$p0[x$states != '']
}
}
if (!is.null(x$lower)) {
x$lower <- 1- x$lower
x$upper <- 1- x$upper
}
if (missing(fun)) fun <- "event"
}
if (missing(firsty) && !is.null(x$p0)) firsty <- 1-x$p0
ssurv <- as.matrix(x$surv)
stime <- x$time
if( !is.null(x$upper)) {
supper <- as.matrix(x$upper)
slower <- as.matrix(x$lower)
}
else {
conf.int <- FALSE
supper <- NULL #marker for later code
}
# 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)
firsty <- matrix(firsty, nrow=nstrat, ncol=ncol(ssurv))
if (!missing(xmax) && any(x$time>xmax)) {
# prune back the survival curves
# I need to replace x's over the limit with xmax, and y's over the
# limit with either the prior y value or firsty
keepx <- keepy <- NULL # lines to keep
tempn <- table(stemp)
offset <- cumsum(c(0, tempn))
for (i in 1:nstrat) {
ttime <-stime[stemp==i]
if (all(ttime <= xmax)) {
keepx <- c(keepx, 1:tempn[i] + offset[i])
keepy <- c(keepy, 1:tempn[i] + offset[i])
}
else {
bad <- min((1:tempn[i])[ttime>xmax])
if (bad==1) { #lost them all
if (!is.na(firstx)) { # and we are plotting lines
keepy <- c(keepy, 1+offset[i])
ssurv[1+offset[i],] <- firsty[i,]
}
}
else keepy<- c(keepy, c(1:(bad-1), bad-1) + offset[i])
keepx <- c(keepx, (1:bad)+offset[i])
stime[bad+offset[i]] <- xmax
x$n.event[bad+offset[i]] <- 1 #don't plot a tick mark
}
}
# ok, now actually prune it
stime <- stime[keepx]
stemp <- stemp[keepx]
x$n.event <- x$n.event[keepx]
if (!is.null(x$n.censor)) x$n.censor <- x$n.censor[keepx]
ssurv <- ssurv[keepy,,drop=FALSE]
if (!is.null(supper)) {
supper <- supper[keepy,,drop=FALSE]
slower <- slower[keepy,,drop=FALSE]
}
}
#stime <- stime/xscale #scaling is deferred until xmax processing is done
if (!missing(fun)) {
if (is.character(fun)) {
tfun <- switch(fun,
'log' = function(x) x,
'event'=function(x) 1-x,
'cumhaz'=function(x) -log(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,
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)
}
firsty <- tfun(firsty)
}
if (missing(firstx)) {
if (!is.null(x$start.time))
firstx <- x$start.time
else {
if (xlog) firstx <- min(stime[stime>0])
else firstx <- min(0, stime)
}
}
# The default for plot and lines is to add confidence limits
# if there is only one curve
if (missing(conf.int) && missing(conf.times)) conf.int <- (ncurve==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 (is.logical(conf.int)) plot.surv <- TRUE
else {
temp <- match.arg(conf.int, c("both", "only", "none"))
if (is.na(temp)) stop("invalid value for conf.int")
if (temp=="none") conf.int <- FALSE else conf.int <- TRUE
if (temp=="only") plot.surv <- FALSE else plot.surv <- TRUE
}
# Marks are not placed on confidence bands
mark <- rep(mark, 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 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]
yy <- y[censor]
}
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)
censor <- if (is.null(x$n.censor))
(x$n.event[who] ==0) else
(x$n.event[who] ==0 | x$n.censor[who] >0) #NOT censoring ties
xx <- c(firstx, stime[who])
censor <- c(FALSE, censor) #no mark at firstx
yy <- c(firsty[i,j], ssurv[who,j])
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=mark[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, c(firsty[i,j], slower[who,j]), x2,
method='constant', f=1)$y
temphigh<- approx(xx, c(firsty[i,j], supper[who,j]), 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, c(firsty[i,j], slower[who,j])), lty=lty[c2],
col=col[c2],lwd=lwd[c2])
c2 <- c2 +1
lines(dostep(xx, c(firsty[i,j], supper[who,j])), lty=lty[c2],
col=col[c2], lwd= lwd[c2])
c2 <- c2 + 1
}
else {
lines(xx, c(firsty[i,j], slower[who,j]), lty=lty[c2],
col=col[c2],lwd=lwd[c2], type=type)
c2 <- c2 +1
lines(xx, c(firsty[i,j], supper[who,j]), lty=lty[c2],
col=col[c2], lwd= lwd[c2], type= type)
c2 <- c2 + 1
}
}
}
}
invisible(list(x=xend, y=yend))
}
points.survfit <- function(x, xscale, xmax, fun, censor=FALSE,
col=1, pch, ...) {
# this function is used rarely
conf.int <- FALSE # never draw these with 'points'
if (inherits(x, "survfitms")) {
x$surv <- 1- x$pstate
if (is.matrix(x$surv)) {
dimnames(x$surv) <- list(NULL, x$states)
if (ncol(x$surv) > 1 && any(x$states == '')) {
x$surv <- x$surv[, x$states != '']
if (is.matrix(x$p0)) x$p0 <- x$p0[, x$states != '']
else x$p0 <- x$p0[x$states != '']
}
}
if (!is.null(x$lower)) {
x$lower <- 1- x$lower
x$upper <- 1- x$upper
}
if (missing(fun)) fun <- "event"
}
firstx <- firsty <- NA # part of the common args, but irrelevant for points
ssurv <- as.matrix(x$surv)
stime <- x$time
if( !is.null(x$upper)) {
supper <- as.matrix(x$upper)
slower <- as.matrix(x$lower)
}
else {
conf.int <- FALSE
supper <- NULL #marker for later code
}
# 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)
firsty <- matrix(firsty, nrow=nstrat, ncol=ncol(ssurv))
if (!missing(xmax) && any(x$time>xmax)) {
# prune back the survival curves
# I need to replace x's over the limit with xmax, and y's over the
# limit with either the prior y value or firsty
keepx <- keepy <- NULL # lines to keep
tempn <- table(stemp)
offset <- cumsum(c(0, tempn))
for (i in 1:nstrat) {
ttime <-stime[stemp==i]
if (all(ttime <= xmax)) {
keepx <- c(keepx, 1:tempn[i] + offset[i])
keepy <- c(keepy, 1:tempn[i] + offset[i])
}
else {
bad <- min((1:tempn[i])[ttime>xmax])
if (bad==1) { #lost them all
if (!is.na(firstx)) { # and we are plotting lines
keepy <- c(keepy, 1+offset[i])
ssurv[1+offset[i],] <- firsty[i,]
}
}
else keepy<- c(keepy, c(1:(bad-1), bad-1) + offset[i])
keepx <- c(keepx, (1:bad)+offset[i])
stime[bad+offset[i]] <- xmax
x$n.event[bad+offset[i]] <- 1 #don't plot a tick mark
}
}
# ok, now actually prune it
stime <- stime[keepx]
stemp <- stemp[keepx]
x$n.event <- x$n.event[keepx]
if (!is.null(x$n.censor)) x$n.censor <- x$n.censor[keepx]
ssurv <- ssurv[keepy,,drop=FALSE]
if (!is.null(supper)) {
supper <- supper[keepy,,drop=FALSE]
slower <- slower[keepy,,drop=FALSE]
}
}
#stime <- stime/xscale #scaling is deferred until xmax processing is done
if (!missing(fun)) {
if (is.character(fun)) {
tfun <- switch(fun,
'log' = function(x) x,
'event'=function(x) 1-x,
'cumhaz'=function(x) -log(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,
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)
}
firsty <- tfun(firsty)
}
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
}
}
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.