Royston-Parmar Model

A developer demonstrated the utility of Royston-Parmar (RP) models for causal inference in survival analysis, particularly in scenarios with non-proportional hazards where Cox proportional hazards models fail. The code, created with assistance from Google Gemini and reviewed by RB, simulates data with a hazard switch and calculates restricted mean survival time (RMST) to show the flexibility of RP models.

| 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" |