Nothing
survplotp <- function(fit, ...) UseMethod("survplotp")
survplotp.npsurv <-
function(fit, xlim,
ylim, xlab, ylab, time.inc, state=NULL,
conf=c("bands", "none"), mylim=NULL, abbrev.label=FALSE,
col=colorspace::rainbow_hcl,
levels.only=TRUE,
loglog=FALSE, fun=function(y) y, aehaz=FALSE, times=NULL,
logt=FALSE, pr=FALSE, ...) {
conf <- match.arg(conf)
conf.int <- fit$conf.int
if(!length(conf.int) | conf == "none") conf.int <- 0
if(loglog)
fun <- function(y) logb(-logb(ifelse(y == 0 | y == 1, NA, y)))
mstate <- inherits(fit, 'survfitms')
if(mstate) fun <- function(y) 1 - y
## Multi-state model for competing risks
z <- seq(.1, .9, by = .1)
funtype <- if(loglog) 'loglog'
else if(all(fun(z) == z)) 'identity'
else if(all(abs(fun(z) - (1 - z)) < 1e-6)) 'inverse'
else if(loglog) 'loglog'
else 'other'
cylim <- function(ylim)
if(length(mylim)) c(min(ylim[1], mylim[1]), max(ylim[2], mylim[2]))
else ylim
mu <- markupSpecs$html
nbsp <- htmlSpecial('nbsp')
if(funtype %in% c('identity', 'inverse'))
survdiffplotp <-
function(fit, fun, xlim,
conf.int, convert=function(f) f, pobj)
{
if(length(fit$strata) != 2)
stop('must have exactly two strata')
h <- function(level, f) {
i <- f$strata == levels(f$strata)[level]
tim <- f$time[i]
surv <- fun(f$surv[i])
se <- f$std.err[i]
list(time=tim, surv=surv, se=se)
}
times <- sort(c(0, unique(fit$time)))
times <- times[times >= xlim[1] & times <= xlim[2]]
f <- convert(summary(fit, times=times, print.it=FALSE, extend=TRUE))
a <- h(1, f)
b <- h(2, f)
if(! identical(a$time, b$time)) stop('program logic error')
time <- a$time
surv <- (a$surv + b$surv) / 2
se <- sqrt(a$se^2 + b$se^2)
z <- qnorm((1 + conf.int) / 2)
lo <- pmax(0, surv - 0.5 * z * se)
hi <- pmin(1, surv + 0.5 * z * se)
k <- ! is.na(time + lo + hi)
list(times=time[k], lower=lo[k], upper=hi[k])
}
fit.orig <- fit
units <- fit$units
if(!length(units)) units <- "Day"
maxtime <- fit$maxtime
if(! length(maxtime)) maxtime <- max(fit$time)
mintime <- min(fit$time, 0)
pret <- pretty(c(mintime, maxtime))
maxtime <- max(pret)
mintime <- min(pret)
if(missing(time.inc)) {
time.inc <- switch(units, Day=30, Month=1, Year=1,
(maxtime - mintime) / 10)
if(time.inc > maxtime) time.inc <- (maxtime - mintime) / 10
}
if(mstate) {
## Multi-state model for competing risks
if(missing(state))
stop('state must be given when response is a multi-state/competing risk object from Surv()')
if(length(state) != 1) stop('at present state can only be a single state')
states <- fit$states
if(state %nin% states) stop(paste('state is not in',
paste(states, collapse=', ')))
}
if(missing(ylab))
ylab <- switch(funtype,
identity = 'Survival Probability',
inverse = if(mstate) paste('Cumulative Incidence of',
upFirst(state))
else 'Cumulative Incidence',
loglog = 'log(-log Survival Probability)',
other = '')
un <- fit$units
if(un != '') un <- paste0(un, 's')
if(missing(xlab))
xlab <- if(logt) paste0("log Follow-up Time in ", un)
else
mu$varlabel('Follow-up Time', un)
if(missing(xlim))
xlim <- if(logt) logb(c(maxtime / 100, maxtime)) else c(mintime, maxtime)
convert <- if(mstate) {
istate <- match(state, states)
conv <- function(f, istate) {
f$surv <- 1 - f$pstate [, istate]
f$lower <- 1 - f$lower [, istate]
f$upper <- 1 - f$upper [, istate]
f$std.err <- f$std.err[, istate]
icens <- which(states == '(s0)')
if(! length(icens))
stop('Program logic error: did not find (s0) column with competing risks')
f$n.risk <- f$n.risk[, icens]
if(all(f$n.risk == 0))
stop('program logic error: all n.risk are zero')
f
}
formals(conv) <- list(f=NULL, istate=istate)
conv
}
else function(f) f
fit <- convert(fit)
origsurv <- fit$surv
if(funtype != 'identity') {
fit$surv <- fun(fit$surv)
fit$surv[is.infinite(fit$surv)] <- NA
## handle e.g. logit(1) - Inf would mess up ylim in plot()
if(conf.int > 0) {
fit$lower <- fun(fit$lower)
fit$upper <- fun(fit$upper)
fit$lower[is.infinite(fit$lower)] <- NA
fit$upper[is.infinite(fit$upper)] <- NA
if(missing(ylim))
ylim <- cylim(range(c(fit$lower, fit$upper), na.rm=TRUE))
}
else if(missing(ylim)) ylim <- cylim(range(fit$surv, na.rm=TRUE))
}
else if(missing(ylim)) ylim <- c(0, 1)
olev <- slev <- names(fit$strata)
if(levels.only) slev <- gsub('.*=', '', slev)
sleva <- if(abbrev.label) abbreviate(slev) else slev
ns <- length(slev)
slevp <- ns > 0
ns <- max(ns, 1)
if(is.function(col)) col <- col(ns)
y <- 1 : ns
strat <- if(ns == 1) rep(1, length(fit$time))
else rep(1 : ns, fit$strata)
stime <- sort(unique(c(0, fit.orig$time)))
stime <- stime[stime >= mintime & stime <= maxtime]
# v <- convert(summary(fit.orig, times=stime, print.it=FALSE))
# vs <- if(ns > 1) as.character(v$strata)
## survival:::summary.survfit was not preserving order of strata levels
nevents <- totaltime <- numeric(ns)
cuminc <- character(ns)
p <- plotly::plot_ly(...)
pl <- function(p, x, y, n.risk=NULL, col, slev, type='est') {
sname <- if(ns == 1) '' else slev
snames <- if(sname == '') '' else paste0(sname, ' ')
d <- paste0('Difference<br>', mu$half(), ' ',
conf.int, ' CL')
nam <- switch(type,
est = sname,
lower = paste0(snames, conf.int, ' CL'),
upper = paste0(snames, conf.int, ' CL', type),
'diff lower' = d,
'diff upper' = paste0(d, type))
lg <- switch(type,
est = 'Estimates',
lower = paste0(snames, 'CL'),
upper = paste0(snames, 'CL'),
'diff lower' = 'Difference',
'diff upper' = 'Difference')
rx <- format(round(x, 3))
ry <- format(round(y, 3))
txt <- switch(type,
est = paste0('t=', rx,
'<br>Probability=', ry,
if(length(n.risk)) '<br>At risk:', n.risk),
lower = paste0('t=', rx, '<br>Lower:', ry),
upper = paste0('t=', rx, '<br>Upper:', ry),
'diff lower' = NULL,
'diff upper' = NULL)
## Mark in text the point estimates that correspond to every time.inc
if(type == 'est' && ! logt) {
nicet <- seq(mintime, maxtime, by=time.inc)
nicet <- nicet[nicet > 0]
for(ti in nicet) {
if(any(abs(ti - x) < 1e-6)) next # nice time already covered
k <- which(x < ti); k <- max(k)
txt[k] <- paste0(txt[k], '<br>(Also for t=', ti, ')')
}
}
fcol <- plotly::toRGB(col, 0.2)
vis <- if(ns == 2 && type %in% c('lower', 'upper'))
'legendonly' else TRUE
ln <- if(type == 'est') list(shape='hv', color=col)
else list(shape='hv', color=col, width=0)
dat <- if(length(txt)) data.frame(x, y, txt) else data.frame(x, y)
up <- type %in% c('upper', 'diff upper')
p <- if(length(txt))
plotly::add_lines(p, x=~ x, y=~ y, text=~ txt, data=dat,
hoverinfo='text', line=ln,
fillcolor=fcol,
fill=if(type %in% c('upper', 'diff upper'))
'tonexty' else 'none',
visible=vis, legendgroup=lg,
name=nam, showlegend=! up)
else
plotly::add_lines(p, x=~ x, y=~ y, data=dat,
hoverinfo='none', line=ln,
fillcolor=fcol,
fill=if(type %in% c('upper', 'diff upper'))
'tonexty' else 'none',
visible=vis, legendgroup=lg,
name=nam, showlegend=! up)
}
for(i in 1 : ns) {
st <- strat == i
time <- fit$time[st]
surv <- fit$surv[st]
lower <- fit$lower[st]
upper <- fit$upper[st]
osurv <- origsurv[st]
n.risk <- fit$n.risk[st]
if(! logt && xlim[1] ==0 && all(time > xlim[1])) {
time <- c(xlim[1], time)
surv <- c(fun(1), surv)
lower <- c(fun(1), lower)
upper <- c(fun(1), upper)
osurv <- c(1, osurv)
n.risk <- c(fit$n[i], n.risk)
}
## nevents[i] <- sum(fit$n.event[st])
## nrsk <- fit$n.risk[st]
## neachtime <- c(- diff(nrsk), min(nrsk))
## totaltime[i] <- sum(neachtime * time)
nevents[i] <-
if(mstate) {
if(ns == 1) fit$numevents[, state]
else fit$numevents[olev[i], state]
} else {
if(ns == 1) fit$numevents else fit$numevents[olev[i]]
}
totaltime[i] <- if(ns == 1) fit$exposure else fit$exposure[olev[i]]
if(length(times)) {
cumi <- 1. - approx(time, osurv, xout=times, method='constant')$y
noun <- units %in% c('', ' ')
cuminc[i] <- paste(round(cumi, 3), collapse=', ')
}
if(logt) time <- logb(time)
##don't let step function go beyond x-axis -
##this cuts it off but allows step to proceed axis end
if(max(time) > xlim[2]) {
srvl <- surv[time <= xlim[2] + 1e-6]
s.last <- srvl[length(srvl)]
k <- time < xlim[2]
time <- c(time[k], xlim[2])
surv <- c(surv[k], s.last)
n.risk <- c(n.risk[k], n.risk[length(srvl)])
if(conf.int > 0) {
low.last <- lower[time <= xlim[2] + 1e-6]
low.last <- low.last[length(low.last)]
up.last <- upper[time <= xlim[2] + 1e-6]
up.last <- up.last[length(up.last)]
lower <- c(lower[k], low.last)
upper <- c(upper[k], up.last)
}
}
if(logt) p <- pl(p, time, surv, n.risk, col=col[i], slev=sleva[i])
else {
xxx <- time
yyy <- surv
nr <- n.risk
if(mintime < min(time)) {
xxx <- c(mintime, time)
yyy <- c(fun(1), surv)
nr <- c(fit$n[i], n.risk)
}
p <- pl(p, xxx, yyy, nr, col=col[i], slev=sleva[i])
}
if(pr) {
zest <- rbind(time, surv)
dimnames(zest) <- list(c("Time", "Survival"),
rep("", length(time)))
if(slevp)cat("\nEstimates for ", slev[i], "\n\n")
print(zest, digits=3)
}
if(conf.int > 0) {
if(logt) {
p <- pl(p, time, lower, type='lower', col=col[i], slev=sleva[i])
p <- pl(p, time, upper, type='upper', col=col[i], slev=sleva[i])
}
else {
p <- pl(p, c(min(time), time), c(fun(1), lower),
col=col[i], slev=slev[i], type='lower') # see survplot ?max(tim)?
p <- pl(p, c(min(time), time), c(fun(1), upper),
col=col[i], slev=slev[i], type='upper')
}
}
}
if(funtype %in% c('identity', 'inverse') && ns == 2 && conf.int > 0) {
z <- survdiffplotp(fit.orig, fun=fun, conf.int=conf.int,
convert=convert, xlim=xlim, pobj=p)
g <- plotly::toRGB('gray')
p <- pl(p, z$time, z$lower, type='diff lower', col=g, slev='')
p <- pl(p, z$time, z$upper, type='diff upper', col=g, slev='')
}
slevat <- ifelse(sleva == '', '', paste0(sleva, ': '))
if(aehaz || length(times)) {
un <- if(units == ' ' | units == '') ''
else
paste('/', tolower(units), sep='')
haz <- round(nevents / totaltime, 4)
txt <- paste(nevents, 'events')
if(aehaz) txt <- paste0(txt, '<br>', htmlGreek('lambda'), '=',
haz, un, sep='')
z <- paste(paste0(slevat, txt), collapse='<br>')
if(length(times)) {
z2 <- paste0('<br>Cumulative<br>Incidence at<br>',
't=', paste(times, collapse=', '), ' ', units,
if(un !='') 's', '<br>',
paste0(slevat, cuminc, collapse='<br>'))
z <- paste0(z, z2)
}
}
else z <- paste(paste0(slevat, nevents, ' events'), collapse='<br>')
## Add empty trace just to add to bottom of legend. Used to have x=~NA y=~NA
## but plotly update made that point ignored in every way
dam <- data.frame(x=xlim[1], y=ylim[1])
p <- plotly::add_markers(p, x=~ x, y=~ y, # mode='markers',
marker=list(symbol='asterisk'), # suppresses pt
name=z, data=dam)
xaxis <- list(range=xlim, title=xlab)
if(! logt) xaxis <-
c(xaxis,
list(tickvals = seq(xlim[1], max(pretty(xlim)), time.inc)))
plotly::layout(p,
xaxis=xaxis,
yaxis=list(range=ylim, title=ylab))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.