| # Title: Royston-Parmar Models for Causal Inference | | | # Description: RP models can be very useful for causal inference. The goal of this code is to | | | # demonstrate why they can be so useful. | | | # This example uses a non-proportional hazard, to demonstrate how they can be useful | | | # and how they are flexible. In this case, where a Cox PH model would fail. | | | # Additionally, the restricted mean survival time (RMST) is calculated to demonstrate this | | | # Note: This code was created with the help of Google Gemini, and reviewed by me (RB) | |
| # Setup ---- | |
| #... Packages ---- | |
| library(tidyverse) # ol faithful | |
| library(flexsurv) # flexsurvspline() | |
| library(survival) # Surv(), survfit() | |
| library(patchwork) # combining plots | | | # Simulate Data ---- | | | # Here we want to simulate data to demonstrate non-proportional hazards. | | | # Note: here we are not using any confounders. | | | set.seed(2028) # setting seed for reproducibility | |
| n <- 800 # sample size (arbitrary) | |
| group <- rbinom(n, 1, 0.5) # which group (only 2 for this example) | |
| tau_switch <- 0.5 # Hazard switch at 6 months | | | h0_early <- 0.20 | | | h0_late <- 0.40 | |
| b_early <- 1.2 # log(HR = 3.3) -> High early harm | |
| b_late <- -1.8 # log(HR = 0.16) -> High late benefit | |
| lp_early <- b_early * group | | | lp_late <- b_late * group | | | # Simulating the parts needed to estimate the event time. Once we | | | # have the event time, we can know if it was observed, or if it was censored. | |
| u <- runif(n) | |
| neglogu <- -log(u) | |
| H_at_switch <- h0_early * exp(lp_early) * tau_switch | |
| event_time <- ifelse( | |
| neglogu <= H_at_switch, | |
| neglogu / (h0_early * exp(lp_early)), | |
| tau_switch + (neglogu - H_at_switch) / (h0_late * exp(lp_late)) | |
| ) | | | # Based on the event time, we determine if a patient was censored or not. | | | # Note: there is another way to do this using the simsurv package, which | | | # makes simulating survival data much easier. | |
| time_obs <- pmin(event_time, 5, rexp(n, 0.05)) | |
| status <- as.numeric(event_time == time_obs) | |
| dat <- data.frame(time = time_obs, status = status, group = factor(group, labels = c("Standard", "New Rx"))) | |
| # True RMST Difference ---- | |
| # Note: there is a section at the end of this code titled "Bonus: Calculate RMST" | | | # that can be used to get this value. For now, we will use this. | | | # "True" RMST (based on our simulated data) | | | true_pop_diff <- 0.48717 | | | # Fitting Models ---- | | | # We will use 4 years as the cutoff for RMST | | | tau_eval <- 4.0 # cutoff for RMST | |
| t_grid <- seq(0.001, tau_eval, length.out = 300) | |
| #... Kaplan-Meier ---- | |
| km_fit <- survfit(Surv(time, status) ~ group, data = dat) | |
| #... Royston-Parmar Model (Flexible Spline) ---- | |
| # Here we are using the flexsurvspline() function to fit an RP model. | | | # The k is the number of knots used. Here we used 3 but we can use different numbers | | | # depending on our data and model fit. Additionally, we do this on the hazard scale but can be | | | # on other scales as well (i.e. odds, or normal) | |
| # Royston-Parmar Model (Flexible Spline) | |
| rp_nonph <- flexsurvspline( | |
| Surv(time, status) ~ group, | |
| data = dat, | | | k = 3, | | | scale = "hazard", | | | anc = list(gamma1 = ~ group) | | | ) | |
| # Getting Results from the Models ---- | |
| #.. RP model results ---- | |
| rp_rmst_res <- summary(rp_nonph, | |
| newdata = data.frame(group = levels(dat$group)), | |
| type = "rmst", | | | t = tau_eval, | | | ci = TRUE) | | | # Getting the point estimate and standard error | | | # For the control group | |
| est_std_rp <- rp_rmst_res[[1]]$est[1] | |
| se_std_rp <- (rp_rmst_res[[1]]$ucl[1] - rp_rmst_res[[1]]$lcl[1]) / (2 * 1.96) | |
| # For the "Rx" group | |
| est_rx_rp <- rp_rmst_res[[2]]$est[1] | |
| se_rx_rp <- (rp_rmst_res[[2]]$ucl[1] - rp_rmst_res[[2]]$lcl[1]) / (2 * 1.96) | |
| # Calculating the difference and 95% confidence interval | |
| diff_rp <- as.numeric(est_rx_rp - est_std_rp) | |
| se_rp_diff <- as.numeric(sqrt(se_std_rp^2 + se_rx_rp^2)) | |
| lcl_rp <- diff_rp - 1.96 * se_rp_diff | | | ucl_rp <- diff_rp + 1.96 * se_rp_diff | | | #... Format Annotation for Later ---- | | | # Formatted annotation string for plot | |
| rmst_annotation <- sprintf( | |
| "Mean Survival Benefit (Δ RMST at τ = %.1fy):\n • RP Model: %+.2f yrs [95%% CI: %+.2f, %+.2f]\n • True RMST: %+.2f yrs", | |
| tau_eval, diff_rp, lcl_rp, ucl_rp, true_pop_diff | | | ) | |
| # Plots ---- | |
| #... Setting up aesthetics ---- | |
| # Creating a theme that can be used across both plots. | |
| theme_pub <- function() { | |
| theme_minimal(base_size = 18) + | |
| theme( | |
| plot.title = element_text(face = "bold", size = 18, hjust = 0), | |
| plot.subtitle = element_text(color = "grey30", size = 16, margin = margin(b = 6)), | |
| legend.position = "top", | |
| text = element_text(size = 16), | |
| legend.title = element_blank(), | |
| panel.grid.minor = element_blank(), | |
| axis.title = element_text(face = "bold", size = 14) | |
| ) | | | } | | | # Setting the colors | |
| cols <- c("Standard" = "#2b5c8f", "New Rx" = "#d95f02") | |
| #... Plot A: Hazards over Time ---- | |
| # Getting the hazard over time, to plot. This shows how it is non-proportional. | |
| df_true_haz <- data.frame( | |
| time = rep(t_grid, 2), | |
| group = rep(c("Standard", "New Rx"), each = length(t_grid)) | |
| ) %>% | |
| mutate( | |
| haz = case_when( | |
| group == "Standard" & time <= tau_switch ~ h0_early, | |
| group == "Standard" & time > tau_switch ~ h0_late, | |
| group == "New Rx" & time <= tau_switch ~ h0_early * exp(b_early), | |
| group == "New Rx" & time > tau_switch ~ h0_late * exp(b_late) | |
| ) | | | ) | | | # Creating the plot. | |
| p_hazards <- ggplot(df_true_haz, aes(x = time, y = haz, color = group)) + | |
| geom_step(size = 1.1) + | |
| geom_vline(xintercept = tau_switch, linetype = "dashed", color = "grey40") + | |
| annotate("text", x = tau_switch + 0.1, y = max(df_true_haz$haz) * 0.88, | |
| label = "True Hazard Switch\n(Early Harm, Late Benefit)", | |
| hjust = 0, size = 5, fontface = "bold.italic", color = "firebrick") + | |
| scale_color_manual(values = cols) + | |
| scale_y_continuous(limits = c(0, max(df_true_haz$haz) * 1.05)) + | |
| labs( | | | title = "A. Hazards over Time", | |
| subtitle = "Hazards over time (non-proportional hazards)", | |
| x = "Time (Years)", | |
| y = "Hazard Rate h(t)" | |
| ) + | |
| theme_pub() | |
| #... Plot B: Survival Curves ---- | |
| # This has the Kaplan-Meier curves, and the RP models. Additionally it also has the | | | # estimated RMST, and the "true" RMST. | |
| # Kaplan-Meier ---- | |
| km_df <- data.frame( | |
| time = summary(km_fit)$time, | |
| surv = summary(km_fit)$surv, | |
| group = gsub("group=", "", summary(km_fit)$strata) | |
| ) | |
| #... RP model ---- | |
| rp_surv_obj <- summary(rp_nonph, newdata = data.frame(group = levels(dat$group)), t = t_grid, type = "survival", ci = TRUE) | |
| df_rp_surv <- bind_rows( | |
| rp_surv_obj[[1]] %>% mutate(group = "Standard"), | |
| rp_surv_obj[[2]] %>% mutate(group = "New Rx") | |
| ) | |
| #... Creating Plot ---- | |
| p_survival <- ggplot() + | |
| # Raw empirical KM steps | | | geom_step(data = km_df, aes(x = time, y = surv, color = group), alpha = 0.35, size = 0.8) + | | | # RP Model fit with 95% CI ribbons | |
| geom_ribbon(data = df_rp_surv, aes(x = time, ymin = lcl, ymax = ucl, fill = group), alpha = 0.18) + | |
| geom_line(data = df_rp_surv, aes(x = time, y = est, color = group), size = 1.1) + | |
| # Text Box Annotation Card in Lower Left | | | annotate( | | | "label", | | | x = 0.15, y = 0.08, | | | label = rmst_annotation, | | | hjust = 0, vjust = 0, | | | fill = "#f8f9fa", color = "#1a1a1a", | |
| fontface = "bold", size = 5, | |
| label.padding = unit(0.5, "lines"), | |
| label.size = 0.3 | | | ) + | |
| scale_color_manual(values = cols) + | |
| scale_fill_manual(values = cols) + | |
| scale_y_continuous(limits = c(0, 1), labels = scales::percent) + | |
| labs( | | | title = "B. Survival Curves and Royston-Parmar Model", | | | subtitle = "Solid lines (RP spline model) track empirical KM steps", | |
| x = "Time (Years)", | |
| y = "Survival S(t)" | |
| ) + | |
| theme_pub() | |
| # Combine Plots! ---- | |
| final_dashboard <- p_hazards | p_survival | |
| print(final_dashboard) | |
| # Bonus: Calculating the RMST! ---- | |
| # For this example, the "true" RMST is difficult to see. To do this, we need to | | | # integrate to get the area under the curves, and the difference. The below code | | | # (with the help of Google Gemini, does that) | | | # Closed-form formula for integral of piecewise survival S(t) from 0 to tau | |
| calc_true_rmst <- function(h_early, h_late, t_switch = 0.5, tau = 4.0) { | |
| # Piece 1: Area under S(t) from t = 0 to t = t_switch | |
| area_1 <- (1 - exp(-h_early * t_switch)) / h_early | |
| # Survival probability at the switch point S(t_switch) | |
| s_switch <- exp(-h_early * t_switch) | |
| # Piece 2: Area under S(t) from t = t_switch to t = tau | |
| area_2 <- s_switch * (1 - exp(-h_late * (tau - t_switch))) / h_late | |
| return(area_1 + area_2) | |
| } | | | # 1. Calculate true rates for both groups | | | h_std_early <- h0_early | | | h_std_late <- h0_late | |
| h_rx_early <- h0_early * exp(b_early) | |
| h_rx_late <- h0_late * exp(b_late) | |
| # 2. Compute exact RMSTs at tau = 4.0 | | | true_rmst_std <- calc_true_rmst(h_std_early, h_std_late, tau_switch, tau_eval) | | | true_rmst_rx <- calc_true_rmst(h_rx_early, h_rx_late, tau_switch, tau_eval) | | | exact_diff_pop <- true_rmst_rx - true_rmst_std | | | # 3. Print Results to Console | |
| cat("\n=================================================================\n") | |
| cat(sprintf(" EXACT POPULATION TRUTH AT TAU = %.1f YEARS \n", tau_eval)) | |
| cat("=================================================================\n") | |
| cat(sprintf("Standard Group RMST: %.5f years (%.2f months)\n", true_rmst_std, true_rmst_std * 12)) | |
| cat(sprintf("New Rx Group RMST: %.5f years (%.2f months)\n", true_rmst_rx, true_rmst_rx * 12)) | |
| cat("-----------------------------------------------------------------\n") | |
| cat(sprintf("True RMST Gain (Δ): %+.5f years (%+.2f months)\n", exact_diff_pop, exact_diff_pop * 12)) | |
| cat("=================================================================\n\n") |