# # TITLE: Conditional Probability Curves of Event Time Distributions # AUTHOR: Takahiro Hasegawa # ORIGINAL DATE: June 8, 2016 # MODIFIED DATE: # REFERENCE: Uno H, Hasegawa T, Cronin AM, Hassett MJ. # Clinically useful presentation of event time distributions # using conditional probability curves. #----------------------------------------------------------------------------- # PARAMETERS: object = Object returned by the survfit class of functions # containing one survival curve # type = Type indicator, # 1 = Type 1 CP curves, and 2 = Type 2 CP curves # t0 = Time from the initiation of treatment at which # patients are event-free for Type 1 # (scalar or vector) # delta = Conditional probability is calculated for # each specified delta value for Type 2 # (scalar or vector) # ci = Indicator variable to control the presence of # pointwise confidence intervals in the plot # (If TRUE and “delta” has only one element, # confidence bands will be presented) # cl = Confidence coefficient of the pointwise confidence intervals # (by default, 0.95 for 95% intervals) # npert = The number of resampling. The default is 1000. # #=========================================== CPcurve <- function(object, type = 2, t0, delta, ci = TRUE, cl = 0.95, npert = 1000, col = c(1 : 8), lty = 1, lwd = 3, xscale = 1, yscale = 1, bty = "l", xaxs = "i", yaxs = "i", xmin = 0, xmax, ymin = 0, ymax = 1, xlim, ylim, xlab = "", ylab = "", ...){ library(flexsurv) condp <- function(x, dt, bs) { (1 - pgengamma(x + dt, mu = bs[1], sigma = exp(bs[2]), Q = bs[3], lower.tail = FALSE) / pgengamma(x, mu = bs[1], sigma = exp(bs[2]), Q = bs[3], lower.tail = FALSE)) * yscale } time <- with(object, c(rep(time, times = n.event), rep(time, times = n.censor))) status <- with(object, c(numeric(sum(n.event)) + 1, numeric(sum(n.censor)))) fit <- flexsurvreg(formula = Surv(time, status) ~ 1, dist = "gengamma") if (missing(xmax)) { xmax <- max(time) } if (missing(xlim)) { xlim <- c(xmin, xmax) / xscale } if (missing(ylim)) { ylim <- c(ymin, ymax) * yscale } else { if (missing(ymax)){ ymax <- ylim[2] / yscale } } if (type == 1 & !missing(t0)){ t0 <- t0[t0 >= 0] if (length(col) < length(t0)) { col <- rep(col, length.out = length(t0)) } else { col <- col[1 : length(t0)] } if (length(lty) < length(t0)) { lty <- rep(lty, length.out = length(t0)) } else { lty <- lty[1 : length(t0)] } if (length(lwd) < length(t0)) { lwd <- rep(lwd, length.out = length(t0)) } else { lwd <- lwd[1 : length(t0)] } if (missing(xlab)){ xlab = "Time since the initiation of treatment" } if (missing(ylab)){ ylab = "Conditional probability of event occurrence" } curve(condp(t0[1], x * xscale - t0[1], coef(fit)), t0[1] / xscale, max(time) / xscale, col = col[1], lty = lty[1], lwd = lwd[1], xlim = xlim, ylim = ylim, bty = bty, xaxs = xaxs, yaxs = yaxs, xlab = xlab, ylab = ylab, ...) if (length(t0) > 1) { for (i in 2 : length(t0)){ curve(condp(t0[i], x * xscale - t0[i], coef(fit)), t0[i] / xscale, max(time) / xscale, col = col[i], lty = lty[i], lwd = lwd[i], add = TRUE) } } if (ci == TRUE){ for (i in 1 : length(t0)){ bootp <- normboot.flexsurvreg(fit, B = npert, transform = TRUE) times <- seq(t0[i], xmax, length = 1000) bootcp <- apply(bootp, 1, function(x) condp(t0[i], times - t0[i], x)) bootcp <- apply(bootcp, 1, function(x) quantile(x, c((1 - cl) / 2, 1 - (1 - cl) / 2))) COL <- col2rgb(col[i]) / 255 for (j in 1 : (length(times) - 1)){ polygon(c(times[j : (j + 1)], rev(times[j : (j + 1)])) / xscale, c(bootcp[1, j : (j + 1)], rev(bootcp[2, j : (j + 1)])), col = rgb(COL[1], COL[2], COL[3], alpha = 0.1), border = NA) } } } } if (type == 2 & !missing(delta)){ delta <- delta[delta > 0] if (length(col) < length(delta)) { col <- rep(col, length.out = length(delta)) } else { col <- col[1 : length(delta)] } if (length(lty) < length(delta)) { lty <- rep(lty, length.out = length(delta)) } else { lty <- lty[1 : length(delta)] } if (length(lwd) < length(delta)) { lwd <- rep(lwd, length.out = length(delta)) } else { lwd <- lwd[1 : length(delta)] } if (missing(xlab)){ xlab = "How long have you been event-free after the initiation of treatment?" } if (missing(ylab)){ ylab = "Probability of event occurrence in future" } curve(condp(x * xscale, delta[1], coef(fit)), xmin / xscale, (max(time) - delta[1]) / xscale, col = col[1], lty = lty[1], lwd = lwd[1], xlim = xlim, ylim = ylim, bty = bty, xaxs = xaxs, yaxs = yaxs, xlab = xlab, ylab = ylab, ...) if (length(delta) > 1) { for (i in 2 : length(delta)){ curve(condp(x * xscale, delta[i], coef(fit)), xmin / xscale, (max(time) - delta[i]) / xscale, col = col[i], lty = lty[i], lwd = lwd[i], add = TRUE) } } legend(xmax / xscale, ymax * yscale, rev(delta), col = rev(col), lty = rev(lty), lwd = rev(lwd), bty = "n", xjust = 1) if (ci == TRUE){ for (i in 1 : length(delta)){ bootp <- normboot.flexsurvreg(fit, B = npert, transform = TRUE) times <- seq(0, max(time) - delta[i], length = 1000) bootcp <- apply(bootp, 1, function(x) condp(times, delta[i], x)) bootcp <- apply(bootcp, 1, function(x) quantile(x, c((1 - cl) / 2, 1 - (1 - cl) / 2))) COL <- col2rgb(col[i]) / 255 for (j in 1 : (length(times) - 1)){ polygon(c(times[j : (j + 1)], rev(times[j : (j + 1)])) / xscale, c(bootcp[1, j : (j + 1)], rev(bootcp[2, j : (j + 1)])), col = rgb(COL[1], COL[2], COL[3], alpha = 0.1), border = NA) } } } } } #=========================================== # # #--- EXAMPLE ---# # library(survival) # # DFS <- read.csv("Example_CP.csv") # DFS1 <- survfit(Surv(time, status) ~ 1, data = DFS[DFS$group == 1,]) # DFS0 <- survfit(Surv(time, status) ~ 1, data = DFS[DFS$group == 0,]) # # # Type 1 CP curves for XELOX # CPcurve(DFS1, type = 1, t0 = c(0 : 6)) # CPcurve(DFS1, type = 1, t0 = 3) # CPcurve(DFS1, type = 1, t0 = 3, ci = FALSE) # # # Type 1 CP curves for FU/FA # CPcurve(DFS0, type = 1, t0 = c(0 : 6), # ylab = "Conditional Prob. (%)", yscale = 100, yaxt = "n") # axis(2, las = 2) # CPcurve(DFS0, type = 1, t0 = 3, ylim = c(0, 0.4)) # CPcurve(DFS0, type = 1, t0 = 3, ci = FALSE, xlim = c(2, 7)) # # # Type 2 CP curves for XELOX # CPcurve(DFS1, delta = c(1 : 6), xscale = 1 / 12, col = rainbow(6)) # CPcurve(DFS1, delta = 3, ymax = 0.5) # CPcurve(DFS1, delta = 3, ci = FALSE) # # # Type 2 CP curves for FU/FA # CPcurve(DFS0, delta = c(1 : 6), yscale = 100) # CPcurve(DFS0, delta = 3, ymax = 0.5) # CPcurve(DFS0, delta = 3, ci = FALSE, xmax = 5) # # #