This chapter brings regression models into the tidy workflow. The companion materials cover Cox regression tables, forest plots, predicted survival, diagnostic graphics, and Fine–Gray regression for competing risks.
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
The official gtsummary regression tutorial covers model tables and their customization. The tidycmprsk regression reference describes Fine–Gray models; their subdistribution hazard ratios have a different interpretation from Cox cause-specific hazard ratios.
R code
Show the code
# -------------------------------------------# Survival Analysis: Module 4 Code# -------------------------------------------# Load necessary packageslibrary(tidyverse) # Load tidyverse packageslibrary(survival) # Load survival packagelibrary(broom) # Load broom packagelibrary(gtsummary) # Load gtsummary packagelibrary(ggsurvfit) # Load ggsurvfit packagelibrary(survminer) # Load survminer packagelibrary(tidycmprsk) # Load tidycmprsk package# ---------------------------# 1. Load and Prepare Data# ---------------------------# Load GBC datasetgbc <-read.table("data/gbc.txt", header =TRUE) # Load GBC dataset# Reformat the datadf <- gbc |># calculate time to first event (relapse or death)group_by(id) |># group by idarrange(time) |># sort rows by timeslice(1) |># get the first row within each idungroup() |># remove groupingmutate(age40 =ifelse(age >=40, 1, 0), # create binary variable for age >= 40grade =factor(grade), # convert grade to factorprog = prog /100, # rescale progesterone receptorestrg = estrg /100# rescale estrogen receptor ) # ---------------------------# 2. Fit Cox PH Model# ---------------------------# Model fitting: survival::coxph()cox_fit <-coxph(Surv(time, status) ~ hormone + meno + age40 + grade + size + prog + estrg, data = df)summary(cox_fit) # Print model summary# ---------------------------# 3. Tidy Output with broom# ---------------------------tidy_cox <-tidy(cox_fit) # Tidy the coxph outputtidy_cox # Display the tidy output# ---------------------------# 4. Tabulate Results with gtsummary# ---------------------------cox_tbl <- cox_fit |>tbl_regression( # Create a regression tableexponentiate =TRUE, # Exponentiate coefficients to get hazard ratioslabel =list(hormone ~"Hormone Therapy", # Custom labels meno ~"Menopausal", age40 ~"Older than 40", grade ~"Tumor Grade", size ~"Tumor Size (mm)", prog ~"Progesterone Receptor (100 fmol/ml)", estrg ~"Estrogen Receptor (100 fmol/ml)") ) |>add_global_p() # Add global p-value for categorical variablescox_tbl # Display the regression table# ---------------------------# 5. Customize the Regression Table# ---------------------------cox_tbl |># Start with the regression tablemodify_caption("Cox regression analysis of the German breast cancer study") |># Add captionbold_p() |># Bold significant p-valuesitalicize_levels() # Italicize variable levels# ---------------------------# 6. Fit Accelerated Failure Time (AFT) Model# ---------------------------# Fit a Weibull AFT modelaft_fit <-survreg(Surv(time, status) ~ hormone + meno + age + grade + size + prog + estrg, data = df, dist ="weibull") # specify the Weibull model# ---------------------------# 7. Forest Plot of Hazard Ratios# ---------------------------# Tidy with exponentiated coeffs (HR) and CItidy_cox <-tidy(cox_fit, exponentiate =TRUE, conf.int =TRUE) tidy_cox$term <-recode(tidy_cox$term, # Relabel the variableshormone ="Hormone Therapy", meno ="Menopausal", age40 ="Older than 40", grade2 ="Tumor Grade II vs I",grade3 ="Tumor Grade III vs I",size ="Tumor Size (mm)", prog ="Progesterone (100 fmol/ml)", estrg ="Estrogen (100 fmol/ml)")tidy_cox |># plot of hazard ratios and 95% CIsggplot(aes(y=term, x=estimate, xmin=conf.low, xmax=conf.high)) +geom_pointrange(shape =15) +# plots center point (x) in square and range (xmin, xmax)geom_vline(xintercept=1, linetype =2) +# vertical line at HR=1scale_x_log10("Hazard ratio (95% CI)", # log scale for x-axisbreaks =c(0.5, 1, 2, 4)) +# log scale for x-axisggtitle("Cox Regression Results for GBC Data") +# Add titletheme_classic() +# classic theme for clean looktheme(axis.line.y =element_blank(), # remove y-axis lineaxis.ticks.y =element_blank(), # remove y-axis ticksaxis.text.y =element_text(size =11), # set variable label sizeaxis.title.y =element_blank() # remove y-axis title )# ---------------------------# 8. Prediction and Visualization# ---------------------------# Create new data for predictionnew_data <-data.frame(hormone =2, meno =2, age40 =1, grade =factor(2), size =20, prog =1, estrg =1)# Predict survival probabilities for `newdata`pred_surv <-survfit(cox_fit, newdata = new_data[1, ])tidy_pred_surv <-tidy(pred_surv) # Tidy the survival prediction outputhead(tidy_pred_surv) # Display the first few rows of the tidy output# Visualize predicted survivalpred_fig <- pred_surv |># Pass the survfit objectggsurvfit() +# Main functionadd_confidence_interval() +# Add confidence intervalscale_x_continuous("Time (months)", breaks =seq(0, 84, by =12)) +# x-axis formatscale_y_continuous("Relapse-free survival probability", limits =c(0, 1)) +# y-axis formatggtitle("Predicted Relapse-Free Survival for a GBC Patient") +# Add titletheme_classic() # Classic theme for clean look# Add horizontal grid linespred_fig +theme(panel.grid.major.y =element_line()) # Add horizontal grid lines# ---------------------------# 9. Cox Model Diagnostics# ---------------------------# PH assumption: Schoenfeld residualsph_test <-cox.zph(cox_fit) # Test proportional hazards assumptionggcoxzph(ph_test) # Visualize Schoenfeld residuals# Functional form: Martingale residuals vs linear predictorggcoxdiagnostics(cox_fit, type ="martingale", # martingale on y-axiox.scale ="linear.predictions") # linear predictor on x-axis# Influential points: Deviance residualsggcoxdiagnostics(cox_fit, type ="deviance", # deviance on y-axisox.scale ="observation.id", # observation ID on x-axissline =FALSE) # no smoothed line# ---------------------------# 10. Fine-Gray Model (Competing Risks)# ---------------------------# Load trial datasetdata("trial", package ="tidycmprsk") # Load trial data from tidycmprsk packagehead(trial) # Display the first few rows of the data# Fit Fine-Gray modelfg_fit <-crr(Surv(ttdeath, death_cr) ~ trt + age + marker + stage, # fit FG modelfailcode ="death from cancer", trial) # for death from cancerfg_fit # print the Fine-Gray model fit summary# Extract coefficients and variancecoef(fg_fit) # Extract coefficientsvcov(fg_fit) |>head() # Extract variance-covariance matrix# Tidy FG model outputtidy_fg <-tidy(fg_fit, exponentiate =TRUE, conf.int =TRUE) # Tidy model outputtidy_fg # Display the tidy output# Forest plot for sub-distribution hazard ratiostidy_fg |># plot of sub-distribution hazard ratios and 95% CIsggplot(aes(y=term, x=estimate, xmin=conf.low, xmax=conf.high)) +geom_pointrange() +# plots center point (x) and range (xmin, xmax)geom_vline(xintercept=1, linetype =2) +# vertical line at HR=1scale_x_log10("Sub-distribution hazard ratio (95% CI)") +# log scale for x-axistheme_classic() +# classic theme for clean looktheme(axis.line.y =element_blank(), # remove y-axis lineaxis.ticks.y =element_blank(), # remove y-axis ticksaxis.text.y =element_text(size =11), # set variable label sizeaxis.title.y =element_blank() # remove y-axis title )# FG regression tablefg_tbl <- fg_fit |>tbl_regression(exponentiate =TRUE) |># Create a regression table add_global_p() # Add global p-value for categorical variablesfg_tbl # display the regression table# Model-based prediction for CIFfg_pred <-predict(fg_fit, newdata= trial[1:10, ], times =c(6, 12, 18)) # Predict CIFfg_pred # Display the predicted CIF