# Wrapper for the 'survival' package plots and P-Values (event=TRUE for deceased) using time / event arguments # Author : Sylvain Mareschal # License : GPL 3 http://www.gnu.org/licenses/gpl.html # Update : 2013/05/13 11:58 srv <- function( time, event, group, subset, limit = NA, graph = TRUE, test = TRUE, ... ) { require("survival", quietly=TRUE) # Checks if(!is.numeric(time) || length(time) == 0) stop("'time' must be a numeric vector") if(!is.logical(event) || length(event) == 0) stop("'event' must be a logical vector") if(length(time) != length(event) || length(event) != length(group)) stop("'time' (", length(time), "), 'event' (", length(event), ") and 'group' (", length(group), ") must have the same length") # List for return output <- NULL # Apply subset if(!missing(subset)) { time <- time[ subset ] event <- event[ subset ] group <- group[ subset ] } # Remove NA Values sel <- (!is.na(time) & !is.na(event) & !is.na(group)) time <- time[ sel ] event <- event[ sel ] group <- group[ sel ] # Remove missing level in factors if(is.factor(group)) group <- factor(group) # Time limitation if(!is.na(limit)) { event[ time > limit ] <- FALSE time[ time > limit ] <- limit } # Survival computation Srv <- Surv(time=time, event=event) # Statistical difference if(isTRUE(test)) { Sdf <- survdiff(Srv~group, rho=0) Pval <- 1 - pchisq( Sdf$chisq , (sum(1 * (Sdf$exp > 0))) - 1 ) output$p.value <- Pval output$method <- "Log-rank" } # Curve object survfit <- survfit(Srv~group) output$survfit <- survfit # Graphics if(isTRUE(graph)) { # Argument extraction args <- list(...) # Default palette if(! "col" %in% names(args)) args$col <- rainbow(length(unique(group)), v=0.5) # Default legend if(! "x" %in% names(args)) args$x <- "topright" if(! "inset" %in% names(args)) args$inset <- 0.02 if(! "bg" %in% names(args)) args$bg <- "#EEEEEE" if(! "legend" %in% names(args)) args$legend <- sprintf("%s (%i)", names(table(group)), table(group)) if(! "pch" %in% names(args)) args$pch <- NA if(! "lty" %in% names(args)) args$lty <- "solid" if(isTRUE(test)) { args$legend = c(args$legend, paste("P", signif(Pval, 3), sep=" = ")) args$col = c(rep(args$col, length.out=length(unique(group))), "#000000") args$pch = c(rep(args$pch, length.out=length(unique(group))), NA) args$lty = c(rep(args$lty, length.out=length(unique(group))), "blank") } # Argument dispatch plotFormals <- c(names(formals(plot.default)), names(formals(getAnywhere(plot.survfit)$objs[[1]])), names(par())) plotFormals <- unique(plotFormals) argPlot <- args[ intersect(names(args), plotFormals) ] argLegend <- args[ intersect(names(args), names(formals(legend))) ] # Plot argPlot$x <- survfit do.call(what="plot", args=argPlot) if(!is.na(limit)) abline(v=limit, lty="dotted") box() # Legend of group do.call(what="legend", args=argLegend) } # Return return(output) }