Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Kaplan Meier - sub groups

By: nzcoops on May 7th, 2012  |  syntax: R  |  size: 7.18 KB  |  views: 276  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. ggkm <- function(sfit,
  2.                  table = TRUE,
  3.                  returns = FALSE,
  4.                  xlabs = "Time",
  5.                  ylabs = "Survival Probability",
  6.                  xlims = c(0,max(sfit$time)),
  7.                  ylims = c(0,1),
  8.                  ystratalabs = NULL,
  9.                  ystrataname = NULL,
  10.                  timeby = 100,
  11.                  main = "Kaplan-Meier Plot",
  12.                  pval = TRUE,
  13.                  subs = NULL,
  14.                  ...) {
  15.  
  16.     #############
  17.     # libraries #
  18.     #############
  19.  
  20.     require(ggplot2)
  21.     require(survival)
  22.     require(gridExtra)
  23.  
  24.     #################################
  25.     # sorting the use of subsetting #
  26.     #################################
  27.  
  28.     times <- seq(0, max(sfit$time), by = timeby)
  29.  
  30.     if(is.null(subs)){
  31.         subs1 <- 1:length(levels(summary(sfit)$strata))
  32.         subs2 <- 1:length(summary(sfit,censored=T)$strata)
  33.         subs3 <- 1:length(summary(sfit,times = times,extend = TRUE)$strata)
  34.     } else{
  35.         for(i in 1:length(subs)){
  36.             if(i==1){
  37.                 ssvar <- paste("(?=.*\\b=",subs[i],sep="")
  38.             }
  39.             if(i==length(subs)){
  40.                 ssvar <- paste(ssvar,"\\b)(?=.*\\b=",subs[i],"\\b)",sep="")
  41.             }
  42.             if(!i %in% c(1, length(subs))){
  43.                 ssvar <- paste(ssvar,"\\b)(?=.*\\b=",subs[i],sep="")
  44.             }
  45.             if(i==1 & i==length(subs)){
  46.                 ssvar <- paste("(?=.*\\b=",subs[i],"\\b)",sep="")
  47.             }
  48.         }
  49.         subs1 <- which(regexpr(ssvar,levels(summary(sfit)$strata), perl=T)!=-1)
  50.         subs2 <- which(regexpr(ssvar,summary(sfit,censored=T)$strata, perl=T)!=-1)
  51.         subs3 <- which(regexpr(ssvar,summary(sfit,times = times,extend = TRUE)$strata, perl=T)!=-1)
  52.     }
  53.  
  54.     ##################################
  55.     # data manipulation pre-plotting #
  56.     ##################################
  57.  
  58.     if(is.null(ystratalabs)) ystratalabs <- as.character(levels(summary(sfit)$strata)[subs1])
  59.     if(is.null(ystrataname)) ystrataname <- "strata"
  60.     m <- max(nchar(ystratalabs))
  61.     times <- seq(0, max(sfit$time), by = timeby)
  62.  
  63.     .df <- data.frame(                      # data to be used in the survival plot
  64.         time = sfit$time[subs2],
  65.         n.risk = sfit$n.risk[subs2],
  66.         n.event = sfit$n.event[subs2],
  67.         surv = sfit$surv[subs2],
  68.         strata = factor(summary(sfit, censored = T)$strata[subs2]),
  69.         upper = sfit$upper[subs2],
  70.         lower = sfit$lower[subs2]
  71.         )
  72.  
  73.     levels(.df$strata) <- ystratalabs       # final changes to data for survival plot
  74.     zeros <- data.frame(time = 0, surv = 1,
  75.                         strata = factor(ystratalabs, levels=levels(.df$strata)),
  76.                         upper = 1, lower = 1)
  77.     .df <- rbind.fill(zeros, .df)
  78.     d <- length(levels(.df$strata))
  79.  
  80.     ###################################
  81.     # specifying plot parameteres etc #
  82.     ###################################
  83.  
  84.     p <- ggplot( .df, aes(time, surv)) +
  85.         geom_step(aes(linetype = strata), size = 0.7) +
  86.         theme_bw() +
  87.         opts(axis.title.x = theme_text(vjust = 0.5)) +
  88.         scale_x_continuous(xlabs, breaks = times, limits = xlims) +
  89.         scale_y_continuous(ylabs, limits = ylims) +
  90.         opts(panel.grid.minor = theme_blank()) +
  91.         opts(legend.position = c(ifelse(m < 10, .28, .35),ifelse(d < 4, .25, .35))) +    # MOVE LEGEND HERE [first is x dim, second is y dim]
  92.         opts(legend.key = theme_rect(colour = NA)) +
  93.         labs(linetype = ystrataname) +
  94.         opts(plot.margin = unit(c(0, 1, .5,ifelse(m < 10, 1.5, 2.5)),"lines")) +
  95.         opts(title = main)
  96.  
  97.     ## Create a blank plot for place-holding
  98.     ## .df <- data.frame()
  99.     blank.pic <- ggplot(.df, aes(time, surv)) +
  100.         geom_blank() + theme_bw() +
  101.         opts(axis.text.x = theme_blank(),axis.text.y = theme_blank(),
  102.              axis.title.x = theme_blank(),axis.title.y = theme_blank(),
  103.              axis.ticks = theme_blank(),
  104.              panel.grid.major = theme_blank(),panel.border = theme_blank())
  105.  
  106.     #####################
  107.     # p-value placement #
  108.     #####################a
  109.  
  110.     if(pval) {
  111.         sdiff <- survdiff(eval(sfit$call$formula), data = eval(sfit$call$data))
  112.         pval <- pchisq(sdiff$chisq,length(sdiff$n) - 1,lower.tail = FALSE)
  113.         pvaltxt <- ifelse(pval < 0.0001,"p < 0.0001",paste("p =", signif(pval, 3)))
  114.         p <- p + annotate("text",x = 0.6 * max(sfit$time),y = 0.1,label = pvaltxt)
  115.     }
  116.  
  117.     ###################################################
  118.     # Create table graphic to include at-risk numbers #
  119.     ###################################################
  120.  
  121.     if(table) {
  122.         risk.data <- data.frame(
  123.             strata = factor(summary(sfit,times = times,extend = TRUE)$strata[subs3]),
  124.             time = summary(sfit,times = times,extend = TRUE)$time[subs3],
  125.             n.risk = summary(sfit,times = times,extend = TRUE)$n.risk[subs3]
  126.             )
  127.         risk.data$strata <- factor(risk.data$strata, levels=rev(levels(risk.data$strata)))
  128.  
  129.         data.table <- ggplot(risk.data,aes(x = time, y = strata, label = format(n.risk, nsmall = 0))) +
  130.             #, color = strata)) +
  131.             geom_text(size = 3.5) + theme_bw() +
  132.             scale_y_discrete(breaks = as.character(levels(risk.data$strata)),
  133.                              labels = rev(ystratalabs)) +
  134.                                  # scale_y_discrete(#format1ter = abbreviate,
  135.                                  # breaks = 1:3,
  136.                                  # labels = ystratalabs) +
  137.                                  scale_x_continuous("Numbers at risk", limits = xlims) +
  138.                                  opts(axis.title.x = theme_text(size = 10, vjust = 1),
  139.                                       panel.grid.major = theme_blank(), panel.grid.minor = theme_blank(),
  140.                                       panel.border = theme_blank(),axis.text.x = theme_blank(),
  141.                                       axis.ticks = theme_blank(),axis.text.y = theme_text(face = "bold",hjust = 1))
  142.  
  143.         data.table <- data.table +
  144.             opts(legend.position = "none") + xlab(NULL) + ylab(NULL)
  145.  
  146.         data.table <- data.table +
  147.             opts(plot.margin = unit(c(-1.5, 1, 0.1, ifelse(m < 10, 2.5, 3.5) - 0.28 * m), "lines")) # ADJUST POSITION OF TABLE FOR AT RISK
  148.  
  149.         #######################
  150.         # Plotting the graphs #
  151.         #######################
  152.  
  153.         ## p <- ggplotGrob(p)
  154.         ## p <- addGrob(p, textGrob(x = unit(.8, "npc"), y = unit(.25, "npc"), label = pvaltxt,
  155.         ## gp = gpar(fontsize = 12)))
  156.         grid.arrange(p, blank.pic, data.table, clip = FALSE, nrow = 3,
  157.                      ncol = 1, heights = unit(c(2, .1, .25),c("null", "null", "null")))
  158.  
  159.         if(returns) {
  160.             a <- arrangeGrob(p, blank.pic, data.table, clip = FALSE, nrow = 3,
  161.                              ncol = 1, heights = unit(c(2, .1, .25), c("null", "null", "null")))
  162.             return(a)
  163.         }
  164.     } else {
  165.         ## p <- ggplotGrob(p)
  166.         ## p <- addGrob(p, textGrob(x = unit(0.5, "npc"), y = unit(0.23, "npc"),
  167.         ## label = pvaltxt, gp = gpar(fontsize = 12)))
  168.  
  169.         if(returns) return(p)
  170.     }
  171. }