{"id":513,"date":"2015-08-01T11:31:15","date_gmt":"2015-08-01T10:31:15","guid":{"rendered":"http:\/\/pcool.dyndns.org:8080\/statsbook\/?page_id=513"},"modified":"2025-07-01T19:12:03","modified_gmt":"2025-07-01T18:12:03","slug":"survival-plot","status":"publish","type":"page","link":"https:\/\/pcool.dyndns.org\/index.php\/survival-plot\/","title":{"rendered":"Survival Plot"},"content":{"rendered":"\n

For a more detailed discussion about survival analysis, please refer to the relevant section<\/a>. <\/p>\n\n\n\n

When the failure time is know accurately,  Kaplan-Meier analysis is the preferred method of analysis and Life Table analysis is mostly obsolete. Here, several methods on how to create a survival plot in R are discussed.<\/p>\n\n\n\n

Download the plotsurvival.rda<\/a> dataset for this example. The data frame is called plotsurvival and the variables are: number (patient number), fu (continuous; follow up time), group (categorical data; Hip 1 or Hip 2) and censor (binary outcome variable). To show the data:<\/p>\n\n\n\n

plotsurvival<\/mark>\n   number   fu group censor\n1       1  5.2 Hip 1      0\n2       2  6.3 Hip 1      0\n3       3  7.0 Hip 1      0\n4       4  8.5 Hip 1      0\n5       5  9.4 Hip 1      0\n6       6 10.1 Hip 1      0\n7       7 11.2 Hip 1      0\n8       8 12.0 Hip 1      0\n9       9 12.5 Hip 1      1\n10     10 12.7 Hip 1      0\n11     11  4.5 Hip 2      1\n12     12  5.8 Hip 2      0\n13     13  6.9 Hip 2      0\n14     14  7.8 Hip 2      1\n15     15  8.2 Hip 2      1\n16     16  8.9 Hip 2      1\n17     17  9.5 Hip 2      1\n18     18 10.2 Hip 2      0\n19     19 11.5 Hip 2      0\n20     20 12.9 Hip 2      0<\/mark><\/em><\/code><\/pre>\n\n\n\n
To perform survival analysis, the survival package1<\/a><\/sup> should be installed<\/a>.<\/p>\n\n\n\n
library(survival)<\/span><\/em><\/code><\/pre>\n\n\n\n
Method 1<\/strong>, in R console:<\/p>\n\n\n\n
First load the survival library and create a survival object called ‘hipsurvival’:<\/em><\/p>\n\n\n\n
library(survival)<\/em><\/span>\nhipsurvival <- Surv(plotsurvival<\/em><\/span><\/em><\/span>$<\/mark><\/em>fu, plotsurvival<\/em><\/span><\/em><\/span>$<\/mark><\/em>censor)<\/em><\/span>\nhipsurvival<\/mark>\n [1]  5.2+  6.3+  7.0+  8.5+  9.4+ 10.1+ 11.2+ 12.0+ 12.5  12.7+  4.5   5.8+  6.9+  7.8   8.2   8.9   9.5  10.2+ 11.5+ 12.9+<\/mark><\/em><\/code><\/pre>\n\n\n\n
A summary can be requested:<\/p>\n\n\n\n
summary(hipsurvival)<\/mark>\n      time            status   \n Min.   : 4.500   Min.   :0.0  \n 1st Qu.: 6.975   1st Qu.:0.0  \n Median : 9.150   Median :0.0  \n Mean   : 9.055   Mean   :0.3  \n 3rd Qu.:11.275   3rd Qu.:1.0  \n Max.   :12.900   Max.   :1.0  <\/mark><\/em><\/code><\/pre>\n\n\n\n
Next, create a survival curve with 95% confidence intervals from the survival object and call this ‘hipsurvivalcurve’ :<\/p>\n\n\n\n
hipsurvivalcurve<-survfit(hipsurvival~1,conf.int=0.95)<\/em><\/span><\/code><\/pre>\n\n\n\n
The ~1 (tilde one) is required for the function; and means no grouping; otherwise, the grouping variable comes after the tilde (see later).<\/p>\n\n\n\n
To plot the curve:<\/p>\n\n\n\n
plot(hipsurvivalcurve)<\/em><\/span><\/code><\/pre>\n\n\n\n
and add a title and axis labels:<\/p>\n\n\n\n
title(main='Hip Survival Curve',xlab='Follow Up [years]',ylab='probability'<\/em>)<\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
To create a plot grouped by type of hip replacement is very straight forward. Instead of using the ~1 in the survival curve object, use:<\/p>\n\n\n\n
hipsurvivalcurve2 <- survfit(hipsurvival~plotsurvival$group<\/strong>,conf.int=0.95)<\/span><\/em><\/code><\/pre>\n\n\n\n
A new plot can be created by plotting the object, with the first curve in red and the second in green :<\/p>\n\n\n\n
plot(hipsurvivalcurve2,col=c('red','green'))<\/span><\/em>\ntitle(main='Hip Survival Curve',xlab='Follow Up [years]',ylab='probability')<\/span><\/em><\/code><\/pre>\n\n\n\n
(if there are more groups, more colours can be added)<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Method 2<\/strong>, curves with numbers at risk<\/p>\n\n\n\n
It is preferable to have the number of patients at risk indicated in the plot. The KMatRisk<\/a> function can be used for this. Copy the function to the clipboard, paste it into the console and hit return. The function should now be available in R (the library survival should also be loaded).<\/p>\n\n\n\n
Create the survival curve object directly (in one step) from the data and call this hipfit:<\/p>\n\n\n\n
hipfit <- survfit(Surv(fu,censor)~group,data=plotsurvival)<\/span><\/em><\/code><\/pre>\n\n\n\n
Next call the loaded function to create the survival curve with numbers at risk:<\/p>\n\n\n\n
ggkm(hipfit,timeby=1, ystratalabs=c(\"Hip 1\",\"Hip 2\"), main=\"KM plot, Hip by type\",xlabs=\"FU [Years]\")<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Method 3, <\/strong>curves using the product limit package2<\/a><\/sup>:<\/p>\n\n\n\n
The prodlim package creates nice survival curves. First install the package:<\/p>\n\n\n\n
library(prodlim)<\/span><\/em><\/code><\/pre>\n\n\n\n
Create a Kaplan-Meier object using the product limit method:<\/p>\n\n\n\n
kmhipfit <- prodlim(Hist(fu, censor) ~ 1, data = plotsurvival)<\/span><\/em><\/code><\/pre>\n\n\n\n
To create a plot that shows numbers at risk at yearly intervals with a shadowed 95% confidence interval:<\/p>\n\n\n\n
plot(kmhipfit, percent=FALSE, axes=TRUE, axis1.at=seq(0,kmhipfit$maxtime+1,1), axis1.lab=seq(0,kmhipfit$maxtime+1,1), marktime=TRUE, xlab=\"Years\", atrisk=TRUE, atrisk.labels=\"At Risk:\", confint=TRUE, confint.citype=\"shadow\", col=4)<\/span><\/em><\/code><\/pre>\n\n\n\n
And add a title:<\/p>\n\n\n\n
title(main= \"Kaplan Meier Survival Analysis - All Hips\")<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Similarly, create a plot grouped by type of hip replacement (instead of ~1, use ~group):<\/p>\n\n\n\n
kmhipfit2 <- prodlim(Hist(fu, censor) ~ group, data = plotsurvival)<\/em><\/span>\nplot(kmhipfit2, percent=FALSE, axes=TRUE, axis1.at=seq(0,kmhipfit2$maxtime+1,1), axis1.lab=seq(0,kmhipfit2$maxtime+1,1), marktime=TRUE, xlab=\"Years\", atrisk=TRUE, atrisk.labels=\"At Risk:\", confint=TRUE, confint.citype=\"shadow\", col=c(4,3), legend=TRUE, legend.x=0, legend.y=0.4, legend.cex=1)<\/em><\/span>\ntitle(main= \"Kaplan Meier Survival Analysis - By Hip Type\")<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
For more information re prodlim options:<\/p>\n\n\n\n
??prodlim<\/span><\/em><\/p>\n\n\n\n
Please note the syntax to place the legend position.<\/p>\n\n\n\n
Method 4<\/strong>, using the survminer package3<\/a><\/sup> to create survival curves with ggplot24<\/a><\/sup>:<\/p>\n\n\n\n
Still using the plotsurvival.rda<\/a> dataset in this example, make sure the appropriate packages are available:<\/p>\n\n\n\n
library(survival)<\/span><\/em>\nlibrary(ggplot2)<\/span><\/em>\nlibrary(survminer)<\/span><\/em><\/code><\/pre>\n\n\n\n
First create a survival object and examine it:<\/p>\n\n\n\n
fit <- survfit(Surv(fu, censor) ~ group, data = plotsurvival)<\/span><\/em>\nfit<\/span><\/em>\nCall: survfit(formula = Surv(fu, censor) ~ group, data = plotsurvival)\n\n             n events median 0.95LCL 0.95UCL\ngroup=Hip 1 10      1   12.5    12.5      NA\ngroup=Hip 2 10      5    9.5     8.2      NA<\/mark><\/em>\n\nsummary(fit)<\/em><\/span>\nCall: survfit(formula = Surv(fu, censor) ~ group, data = plotsurvival)\n\n                group=Hip 1 \n        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI \n      12.500        2.000        1.000        0.500        0.354        0.125        1.000 \n\n                group=Hip 2 \n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  4.5     10       1    0.900  0.0949        0.732         1.00\n  7.8      7       1    0.771  0.1442        0.535         1.00\n  8.2      6       1    0.643  0.1679        0.385         1.00\n  8.9      5       1    0.514  0.1769        0.262         1.00\n  9.5      4       1    0.386  0.1732        0.160         0.93<\/mark><\/em><\/code><\/pre>\n\n\n\n
To create a basic survival plot:<\/p>\n\n\n\n
ggsurvplot(fit)<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Or a more advanced plot:<\/p>\n\n\n\n
ggsurvplot(fit, legend = 'bottom', <\/span><\/em>\nlegend.title = 'Group', <\/span><\/em>\nlegend.labs = c('Hip 1', 'Hip 2'), <\/span><\/em>\npalette = c(\"#E7B800\", \"#2E9FDF\"), <\/span><\/em>\nconf.int = TRUE, <\/span><\/em>\nconf.int.style = 'ribbon', <\/span><\/em>\npval = TRUE, <\/span><\/em>\nrisk.table = TRUE, <\/span><\/em>\nmain = 'Survival Curve', <\/span><\/em>\nxlab = 'Follow Up [Months]', <\/span><\/em>\nylab = 'Cumulative Survival Probability [%]', <\/span><\/em>\nsurv.scale = 'percent')<\/span><\/em><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Method 5<\/strong>, using the ggfortify package5<\/a><\/sup> to create survival curves with ggplot26<\/a><\/sup>:<\/p>\n\n\n\n
library(survival)<\/em><\/span>\nlibrary(ggfortify)<\/em><\/span>\nfit <- survfit(Surv(fu, censor) ~ group, data = plotsurvival)<\/em><\/span>\nautoplot(fit)<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Or for a more advanced plot, ggplot27<\/a><\/sup> calls can be added:<\/p>\n\n\n\n
autoplot(fit) + <\/em><\/span>\nxlab(\"Tme [years]\") + <\/em><\/span>\nylab(\"Survival Probability [%]\") + <\/em><\/span>\nggtitle(\"Kaplan Meier Survival\") +<\/em><\/span>\ntheme_bw(base_size=14, base_family=\"Roboto\")<\/em><\/span><\/code><\/pre>\n\n\n\n
\"\"<\/figure>\n\n\n\n
Method 6<\/strong>, using the ggfortify package8<\/a><\/sup> to create bespoke survival curves<\/em> with ggplot2 9<\/a><\/sup> and the RcmdrPlugin.KMggplot2 package10<\/a><\/sup>:<\/p>\n\n\n\n
install.packages(\"RcmdrPlugin.KMggplot2\")<\/em><\/span><\/code><\/pre>\n\n\n\n
However, there is no need to load the package with the library command. The survival and ggfortify packages should be loaded:<\/p>\n\n\n\n
library(survival)<\/span><\/em>\nlibrary(ggfortify)<\/span><\/em>\n\nfit <- survfit(Surv(fu, censor) ~ group, data = plotsurvival)<\/em><\/span>\nfit<\/em><\/span>\nCall: survfit(formula = Surv(fu, censor) ~ group, data = plotsurvival)\n\n             n events median 0.95LCL 0.95UCL\ngroup=Hip 1 10      1   12.5    12.5      NA<\/mark><\/em>\n\nsummary(fit)<\/em><\/span>\nCall: survfit(formula = Surv(fu, censor) ~ group, data = plotsurvival)\n\n                group=Hip 1 \n        time       n.risk      n.event     survival      std.err lower 95% CI upper 95% CI \n      12.500        2.000        1.000        0.500        0.354        0.125        1.000 \n\n                group=Hip 2 \n time n.risk n.event survival std.err lower 95% CI upper 95% CI\n  4.5     10       1    0.900  0.0949        0.732         1.00\n  7.8      7       1    0.771  0.1442        0.535         1.00\n  8.2      6       1    0.643  0.1679        0.385         1.00\n  8.9      5       1    0.514  0.1769        0.262         1.00\n  9.5      4       1    0.386  0.1732        0.160         0.93<\/mark><\/em><\/code><\/pre>\n\n\n\n
Convert the survival object to a data frame with the fortify() function:<\/p>\n\n\n\n
plot_data <- fortify(fit, surv.connect = TRUE)<\/em><\/span>\nplot_data<\/em><\/span>\n   time n.risk n.event n.censor      surv   std.err     upper     lower strata\n1   0.0     10       0        0 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n2   0.0     10       0        0 1.0000000 0.0000000 1.0000000 1.0000000  Hip 2\n3   5.2     10       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n4   6.3      9       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n5   7.0      8       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n6   8.5      7       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n7   9.4      6       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n8  10.1      5       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n9  11.2      4       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n10 12.0      3       0        1 1.0000000 0.0000000 1.0000000 1.0000000  Hip 1\n11 12.5      2       1        0 0.5000000 0.7071068 1.0000000 0.1250488  Hip 1\n12 12.7      1       0        1 0.5000000 0.7071068 1.0000000 0.1250488  Hip 1\n13  4.5     10       1        0 0.9000000 0.1054093 1.0000000 0.7320116  Hip 2\n14  5.8      9       0        1 0.9000000 0.1054093 1.0000000 0.7320116  Hip 2\n15  6.9      8       0        1 0.9000000 0.1054093 1.0000000 0.7320116  Hip 2\n16  7.8      7       1        0 0.7714286 0.1868706 1.0000000 0.5348489  Hip 2\n17  8.2      6       1        0 0.6428571 0.2612546 1.0000000 0.3852425  Hip 2\n18  8.9      5       1        0 0.5142857 0.3438807 1.0000000 0.2621155  Hip 2\n19  9.5      4       1        0 0.3857143 0.4489847 0.9299128 0.1599887  Hip 2\n20 10.2      3       0        1 0.3857143 0.4489847 0.9299128 0.1599887  Hip 2\n21 11.5      2       0        1 0.3857143 0.4489847 0.9299128 0.1599887  Hip 2\n22 12.9      1       0        1 0.3857143 0.4489847 0.9299128 0.1599887  Hip 2<\/mark><\/code><\/pre>\n\n\n\n
Plot:<\/p>\n\n\n\n
ggplot(plot_data, aes(x = time, y = surv * 100)) + <\/em><\/span>\n# plot the survival curve(s)<\/em><\/span>\ngeom_step(aes(colour = strata), size = 2) + <\/em><\/span>\n# add the censor indicators<\/em><\/span>\ngeom_point(data = subset(plot_data, n.censor > 0), shape = 3, size = 3, stroke = 2, aes(colour = strata)) + <\/em><\/span>\n# add confidence intervals using the stepribbon geom in the RcmdrPlugin.KMggplot2 package<\/em><\/span>\nRcmdrPlugin.KMggplot2::geom_stepribbon(data = plot_data, aes(x = time, ymin = lower * 100, ymax = upper * 100, fill = strata), alpha = 0.15, colour = \"transparent\", show.legend = FALSE, kmplot = TRUE) + <\/em><\/span>\n# set the colours of the confidence intervals<\/em><\/span>\nscale_fill_manual(values = c(\"steelblue\", \"darkorange\")) +<\/em><\/span>\n# set the colours of the survival curves<\/em><\/span>\nscale_colour_manual(values = c(\"steelblue\", \"darkorange\"), name = \"Hip\") + <\/em><\/span>\n# apply a theme based on the Roboto font size 16<\/em><\/span>\ntheme_bw(base_size = 16, base_family = \"Roboto\") +<\/em><\/span>\n# add a title, x label, y label and footnote<\/em><\/span>\nlabs(title = \"Custom Made Survival Curve\", subtitle = \"with ggfortify and ggplot2\", x = \"Time [years]\", y = \"Survival Probability [%]\", caption = \"countcool.com\") +<\/em><\/span>\n# set the caption to a smaller italic font in darkolivergreen colour<\/em><\/span>\ntheme(plot.caption = element_text(colour = \"darkolivegreen\", face = \"italic\", size = 10)) +<\/em><\/span>\n# set the y axis limit in steps of 20 %<\/em><\/span>\nscale_y_continuous(breaks = seq(0, 100, 20)) +<\/em><\/span>\n# set the x axis in steps of 2 years<\/em><\/span>\nscale_x_continuous(breaks = seq(0, max(plot_data$time) + 2, 2)) +<\/em><\/span>\n# add numbers at risk when failed (label creates a border)<\/em><\/span>\nannotate(geom = \"label\", x = subset(plot_data$time, plot_data$n.censor == 0), y = -15, label = subset(plot_data$n.risk, plot_data$n.censor == 0)) +<\/em><\/span>\nannotate(geom = \"text\", x = 0.2, y = -8, label = \"At Risk\", fontface = \"italic\")<\/em><\/span><\/code><\/pre>\n\n\n\n
This should produce the plot below:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/p>\n\n\n\n
Method 7<\/strong>, using the ggfortify package11<\/a><\/sup> to create a function<\/em> for bespoke survival curves<\/em> with ggplot212<\/a><\/sup> and the RcmdrPlugin.KMggplot2 package13<\/a><\/sup>:<\/p>\n\n\n\n
install.packages(\"RcmdrPlugin.KMggplot2\")<\/em><\/span><\/code><\/pre>\n\n\n\n
However, there is no need to load the package with the library command.<\/p>\n\n\n\n
Load the following function:<\/p>\n\n\n\n
cool survival function<\/a><\/p>\n\n\n\n
by copying and pasting the contents in the R console.<\/p>\n\n\n\n
Run it with the following command:<\/p>\n\n\n\n
cool_survival(<\/span>plotsurvival, <\/span>time = <\/span>plotsurvival$<\/span>fu, <\/span>censor = <\/span>plotsurvival$<\/span>censor, <\/span>group = <\/span>plotsurvival$<\/span>group)<\/span><\/span><\/em><\/code><\/pre>\n\n\n\n
This should produce the same plot as before:<\/p>\n\n\n\n
\"\"<\/figure>\n\n\n\n
<\/p>\n\n\n\n
 <\/p>\n","protected":false},"excerpt":{"rendered":"
For a more detailed discussion about survival analysis, please refer to the relevant section. When the failure time is know accurately,  Kaplan-Meier analysis is the preferred method of analysis and Life Table analysis is mostly obsolete. Here, several methods on how to create a survival plot in R are discussed. Download the plotsurvival.rda dataset for […]<\/p>\n","protected":false},"author":2,"featured_media":0,"parent":0,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"inline_featured_image":false,"footnotes":""},"class_list":["post-513","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/513","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/comments?post=513"}],"version-history":[{"count":4,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/513\/revisions"}],"predecessor-version":[{"id":4724,"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/pages\/513\/revisions\/4724"}],"wp:attachment":[{"href":"https:\/\/pcool.dyndns.org\/index.php\/wp-json\/wp\/v2\/media?parent=513"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}