• API
• FAQ
• Tools
• Trends
• Archive
SHARE
TWEET

# Kaplan Meier - sub groups 2

nzcoops Jun 21st, 2012 5,644 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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.     if( !is.null(subs) ) pval <- FALSE
55.
56.     ##################################
57.     # data manipulation pre-plotting #
58.     ##################################
59.
60.     if(is.null(ystratalabs)) ystratalabs <- as.character(sub("group=*","",names(sfit\$strata))) #[subs1]
61.     if(is.null(ystrataname)) ystrataname <- "Strata"
62.     m <- max(nchar(ystratalabs))
63.     times <- seq(0, max(sfit\$time), by = timeby)
64.
65.     .df <- data.frame(                      # data to be used in the survival plot
66.         time = sfit\$time[subs2],
67.         n.risk = sfit\$n.risk[subs2],
68.         n.event = sfit\$n.event[subs2],
69.         surv = sfit\$surv[subs2],
70.         strata = factor(summary(sfit, censored = T)\$strata[subs2]),
71.         upper = sfit\$upper[subs2],
72.         lower = sfit\$lower[subs2]
73.     )
74.
75.     levels(.df\$strata) <- ystratalabs       # final changes to data for survival plot
76.     zeros <- data.frame(time = 0, surv = 1,
77.                         strata = factor(ystratalabs, levels=levels(.df\$strata)),
78.                         upper = 1, lower = 1)
79.     .df <- rbind.fill(zeros, .df)
80.     d <- length(levels(.df\$strata))
81.
82.     ###################################
83.     # specifying plot parameteres etc #
84.     ###################################
85.
86.     p <- ggplot( .df, aes(time, surv)) +
87.         geom_step(aes(linetype = strata), size = 0.7) +
88.         theme_bw() +
89.         opts(axis.title.x = theme_text(vjust = 0.5)) +
90.         scale_x_continuous(xlabs, breaks = times, limits = xlims) +
91.         scale_y_continuous(ylabs, limits = ylims) +
92.         opts(panel.grid.minor = theme_blank()) +
93.         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]
94.         opts(legend.key = theme_rect(colour = NA)) +
95.         labs(linetype = ystrataname) +
96.         opts(plot.margin = unit(c(0, 1, .5,ifelse(m < 10, 1.5, 2.5)),"lines")) +
97.         opts(title = main)
98.
99.     ## Create a blank plot for place-holding
100.     ## .df <- data.frame()
101.     blank.pic <- ggplot(.df, aes(time, surv)) +
102.         geom_blank() + theme_bw() +
103.         opts(axis.text.x = theme_blank(),axis.text.y = theme_blank(),
104.              axis.title.x = theme_blank(),axis.title.y = theme_blank(),
105.              axis.ticks = theme_blank(),
106.              panel.grid.major = theme_blank(),panel.border = theme_blank())
107.
108.     #####################
109.     # p-value placement #
110.     #####################a
111.
112.     if(pval) {
113.         sdiff <- survdiff(eval(sfit\$call\$formula), data = eval(sfit\$call\$data))
114.         pval <- pchisq(sdiff\$chisq,length(sdiff\$n) - 1,lower.tail = FALSE)
115.         pvaltxt <- ifelse(pval < 0.0001,"p < 0.0001",paste("p =", signif(pval, 3)))
116.         p <- p + annotate("text",x = 0.6 * max(sfit\$time),y = 0.1,label = pvaltxt)
117.     }
118.
119.     ###################################################
120.     # Create table graphic to include at-risk numbers #
121.     ###################################################
122.
123.     if(table) {
124.         risk.data <- data.frame(
125.             strata = factor(summary(sfit,times = times,extend = TRUE)\$strata[subs3]),
126.             time = summary(sfit,times = times,extend = TRUE)\$time[subs3],
127.             n.risk = summary(sfit,times = times,extend = TRUE)\$n.risk[subs3]
128.         )
129.         risk.data\$strata <- factor(risk.data\$strata, levels=rev(levels(risk.data\$strata)))
130.
131.         data.table <- ggplot(risk.data,aes(x = time, y = strata, label = format(n.risk, nsmall = 0))) +
132.             #, color = strata)) +
133.             geom_text(size = 3.5) + theme_bw() +
134.             scale_y_discrete(breaks = as.character(levels(risk.data\$strata)),
135.                              labels = rev(ystratalabs)) +
136.                                  # scale_y_discrete(#format1ter = abbreviate,
137.                                  # breaks = 1:3,
138.                                  # labels = ystratalabs) +
139.                                  scale_x_continuous("Numbers at risk", limits = xlims) +
140.                                  opts(axis.title.x = theme_text(size = 10, vjust = 1),
141.                                       panel.grid.major = theme_blank(), panel.grid.minor = theme_blank(),
142.                                       panel.border = theme_blank(),axis.text.x = theme_blank(),
143.                                       axis.ticks = theme_blank(),axis.text.y = theme_text(face = "bold",hjust = 1))
144.
145.         data.table <- data.table +
146.             opts(legend.position = "none") + xlab(NULL) + ylab(NULL)
147.
148.         data.table <- data.table +
149.             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
150.
151.         #######################
152.         # Plotting the graphs #
153.         #######################
154.
155.         ## p <- ggplotGrob(p)
156.         ## p <- addGrob(p, textGrob(x = unit(.8, "npc"), y = unit(.25, "npc"), label = pvaltxt,
157.         ## gp = gpar(fontsize = 12)))
158.         grid.arrange(p, blank.pic, data.table, clip = FALSE, nrow = 3,
159.                      ncol = 1, heights = unit(c(2, .1, .25),c("null", "null", "null")))
160.
161.         if(returns) {
162.             a <- arrangeGrob(p, blank.pic, data.table, clip = FALSE, nrow = 3,
163.                              ncol = 1, heights = unit(c(2, .1, .25), c("null", "null", "null")))
164.             return(a)
165.         }
166.     } else {
167.         ## p <- ggplotGrob(p)
168.         ## p <- addGrob(p, textGrob(x = unit(0.5, "npc"), y = unit(0.23, "npc"),
169.         ## label = pvaltxt, gp = gpar(fontsize = 12)))
170.
171.         if(returns) return(p)
172.     }
173. }
RAW Paste Data
Top