This chapter revisits survival estimation through tidy tools. The companion materials show how to extract estimates, produce tables with gtsummary, draw Kaplan–Meier curves with ggsurvfit, and analyze competing events with tidycmprsk.
Under construction
The narrative and worked explanations are under construction. The original slides and code are available for reference; the code collection has not yet received the same review as Chapters 1 and 2.
Further resources
ggsurvfit documents survival curves, confidence intervals, and risk tables. For competing events, consult tidycmprsk: its event variable is a factor whose first level denotes censoring. Event-specific cumulative incidence requires distinguishing competing events rather than combining them into one binary endpoint.
R code
Show the code
# -------------------------------------------# Survival Analysis: Module 3 Code# -------------------------------------------# Load necessary packageslibrary(tidyverse) # For data manipulationlibrary(survival) # For survival analysis functionslibrary(broom) # For tidying model outputlibrary(gtsummary) # For tabulating survival estimateslibrary(ggsurvfit) # For tidy survival plottinglibrary(tidycmprsk) # For competing risks analysislibrary(patchwork) # For combining ggplots# ---------------------------# 1. Load and Prepare Data# ---------------------------# Load mortality + relapse datagbc <-read.table("data/gbc.txt", header =TRUE)# Prepare dataset: collapse to time-to-first eventdf <- gbc |># Start with raw datasetgroup_by(id) |># Group by subject IDarrange(time) |># Sort by timeslice(1) |># Keep the first event per subjectungroup() # Remove grouping# ---------------------------# 2. Kaplan-Meier Estimation# ---------------------------# Fit Kaplan-Meier curves by hormone therapy groupkm_fit <-survfit(Surv(time, status >0) ~ hormone, data = df)# Summarize KM estimates at specific time pointssummary(km_fit, times =c(6, 12, 24, 36))# ---------------------------# 3. Tabulate Survival Estimates# ---------------------------# Create summary table of survival estimatestbl_survfit( km_fit, # survfit objectlabel ="Hormone", # Row labeltimes =c(6, 12, 24, 36), # Time pointslabel_header ="Month {time}"# Column label)# Tabulate estimates grouped by menopause and tumor gradedf |>tbl_survfit(y =Surv(time, status), # Survival objectinclude =c(meno, grade), # Grouping variableslabel =list(meno ="Menopause", grade ="Tumor grade"), # Row labelstimes =c(6, 12, 24, 36), # Time pointslabel_header ="Month {time}"# Column label )# Tabulate quantiles (median, quartiles)tbl_survfit( km_fit, # survfit objectlabel ="Hormone", # Row labelprobs =c(0.25, 0.5, 0.75), # Quantileslabel_header ="{100 * prob}% quantile"# Column label)# ---------------------------# 4. Visualize KM Curves (Base R)# ---------------------------# Plot KM curves with confidence intervalsplot(km_fit, ylim =c(0, 1),xlab ="Time (months)", # x-axis labelylab ="Relapse-free survival probability", # y-axis labelcol =c("red", "blue"), # Line colorsconf.int =TRUE) # Show confidence intervals# Add legend to the plotlegend("bottomleft", col =c("red", "blue"), lty =1,legend =c("No hormone", "Hormone"))# ---------------------------# 5. Enhanced KM Plot with ggsurvfit# ---------------------------# Fit KM curves with relabeled hormone variablekm_fit2 <-survfit2(Surv(time, status >0) ~ hormone,data = df |>mutate(hormone =if_else(hormone ==1, "No hormone", "Hormone")))# Customize and plot KM curveskm_fig <- km_fit2 |>ggsurvfit() +# Base plotadd_risktable() +# Add risk tableadd_confidence_interval() +# Add confidence bandsadd_pvalue(caption ="Log-rank {p.value}") +# Add log-rank p-valuescale_x_continuous("Time (months)", breaks =seq(0, 84, 12)) +# x-axisscale_y_continuous("Relapse-free survival probability", limits =c(0, 1)) +# y-axistheme_classic() +# Classic ggplot themetheme(legend.position ="top") # Move legend to top# Save KM plot to fileggsave("images/km_fig.png", km_fig, width =7.5, height =5)# ---------------------------# 6. Competing Risks Analysis# ---------------------------# Load sample trial data from tidycmprsk packagedata("trial", package ="tidycmprsk")# Fit CIF using Gray's estimator for competing riskscif_fit <-cuminc(Surv(ttdeath, death_cr) ~ trt, trial)# ---------------------------# 7. Tabulate CIF Estimates# ---------------------------# Tabulate CIF estimates with confidence intervals and p-valuestbl_cuminc( cif_fit, # cuminc objectoutcomes =c("death from cancer", "death other causes"), # Outcomestimes =c(10, 15, 20), # Time pointslabel_header ="Month {time}"# Column label) |>add_p() # Add Gray's test p-values# ---------------------------# 8. Plot CIF Curves# ---------------------------# Plot CIF for cancer-related death with customizationcif_fit |>ggcuminc(outcome ="death from cancer") +# Specify outcomeadd_confidence_interval() +# Confidence bandsadd_risktable() +# Risk tableadd_pvalue(caption ="Gray's test {p.value}") +# P-value annotationscale_x_continuous("Time (months)", breaks =seq(0, 24, 6)) +# x-axisscale_y_continuous("Cumulative incidence function", limits =c(0, 0.5)) +# y-axisggtitle("Death from cancer") +# Titletheme_classic() +# Themetheme(legend.position ="top") # Legend at top# Save CIF plot to fileggsave("images/cif_fig.png", width =7.5, height =5)# ---------------------------# 9. Combine CIF Plots by Outcome# ---------------------------# Define helper function to generate CIF plotscif_plot <-function(cif_fit, outcome){ cif_fit |>ggcuminc(outcome = outcome) +add_confidence_interval() +scale_x_continuous("Time (months)", breaks =seq(0, 24, 6)) +scale_y_continuous("Cumulative incidence function", limits =c(0, 0.5)) +ggtitle(str_to_sentence(outcome)) +theme_classic() +theme(legend.position ="top")}# Create plots for each event typecif_cancer_plot <-cif_plot(cif_fit, "death from cancer")cif_other_plot <-cif_plot(cif_fit, "death other causes")# Combine plots into single displaycif_trial <- cif_cancer_plot + cif_other_plot +plot_layout(guides ="collect") &theme(legend.position ="top")# Save combined plot to fileggsave("images/cif_combined_fig.png", cif_trial, width =8, height =4)