library(ggplot2)
library(dplyr)
library(tidyr)
# 1. Initialize cycle vector: 0 to 1056 months (88 years)
cycles <- 0:1056
# 2. Define Gompertz-based background mortality and SMR multiplier functions
calc_bg_mortality_hazard <- function(t, smr = 4.0) {
# Starting age is 12
age <- 12 + t / 12
# Gompertz approximation for England and Wales (2018-2020)
a <- 0.00003
b <- 0.09
qx <- a * exp(b * age)
qx <- pmin(qx, 0.99) # limit to avoid qx > 1 at extremely high age
# Convert annual probability to monthly probability, then to monthly hazard
q_monthly <- 1 - (1 - qx)^(1/12)
h_monthly <- -log(1 - q_monthly)
# Apply SMR multiplier
return(smr * h_monthly)
}
# Pre-calculate cumulative background SMR-adjusted survival
h_bg_smr <- sapply(cycles, calc_bg_mortality_hazard)
cum_h_bg_smr <- c(0, cumsum(h_bg_smr[-length(h_bg_smr)]))
S_bg_smr <- exp(-cum_h_bg_smr)
# 3. Survival prediction function for parametric models
calc_surv_mcm <- function(t, fit, dist, is_mcm = TRUE) {
if (t == 0) return(1.0)
# Extract parameters from fit
if (is_mcm && inherits(fit, "flexsurvcure")) {
theta <- fit$res["theta", "est"]
if (dist == "llogis") {
shape <- fit$res["shape", "est"]
scale <- fit$res["scale", "est"]
S_u <- flexsurv::pllogis(t, shape = shape, scale = scale, lower.tail = FALSE)
} else if (dist == "lnorm") {
meanlog <- fit$res["meanlog", "est"]
sdlog <- fit$res["sdlog", "est"]
S_u <- plnorm(t, meanlog = meanlog, sdlog = sdlog, lower.tail = FALSE)
}
return(theta + (1 - theta) * S_u)
} else {
# Standard parametric model or fallback flexsurvreg
if (dist == "llogis") {
shape <- fit$res["shape", "est"]
scale <- fit$res["scale", "est"]
return(flexsurv::pllogis(t, shape = shape, scale = scale, lower.tail = FALSE))
} else if (dist == "lnorm") {
meanlog <- fit$res["meanlog", "est"]
sdlog <- fit$res["sdlog", "est"]
return(plnorm(t, meanlog = meanlog, sdlog = sdlog, lower.tail = FALSE))
}
}
return(1.0)
}
# Calculate raw parametric survival for all arms
S_tisa_os_raw <- sapply(cycles, function(t) calc_surv_mcm(t, mcm_tisa_os, "llogis", is_mcm = TRUE))
S_tisa_efs_raw <- sapply(cycles, function(t) calc_surv_mcm(t, mcm_tisa_efs, "llogis", is_mcm = TRUE))
S_blin_os_raw <- sapply(cycles, function(t) calc_surv_mcm(t, mcm_blin_os, "lnorm", is_mcm = TRUE))
S_chemo_os_raw <- sapply(cycles, function(t) calc_surv_mcm(t, mcm_chemo_os, "lnorm", is_mcm = TRUE))
# 4. Generate trace matrices with constraints
build_trace <- function(S_os_raw, S_efs_raw, S_bg_smr) {
# Constrain OS by general population SMR-adjusted survival
S_os_constrained <- pmin(S_os_raw, S_bg_smr)
# Constrain EFS to be <= OS and background survival
S_efs_constrained <- pmin(S_efs_raw, S_os_constrained)
# Compute state occupancies
efs <- S_efs_constrained
pd <- S_os_constrained - S_efs_constrained
dead <- 1 - S_os_constrained
data.frame(EFS = efs, PD = pd, Dead = dead)
}
# Comparator EFS derivation: S_EFS = (S_OS)^(1/0.83) for t <= 60, then flat relative to OS
derive_comparator_efs <- function(S_os_raw) {
S_efs_derived <- numeric(length(cycles))
for (i in seq_along(cycles)) {
t <- cycles[i]
if (t <= 60) {
S_efs_derived[i] <- S_os_raw[i]^(1 / 0.83)
} else {
# Keep flat relative to OS after 5 years (S_EFS(t) = min(S_EFS(60), S_OS(t)))
S_efs_derived[i] <- min(S_os_raw[61]^(1 / 0.83), S_os_raw[i])
}
}
return(S_efs_derived)
}
S_blin_efs_raw <- derive_comparator_efs(S_blin_os_raw)
S_chemo_efs_raw <- derive_comparator_efs(S_chemo_os_raw)
# Build PSM traces (infused/active cohort trace)
trace_tisa_infused <- build_trace(S_tisa_os_raw, S_tisa_efs_raw, S_bg_smr)
trace_blin <- build_trace(S_blin_os_raw, S_blin_efs_raw, S_bg_smr)
trace_chemo <- build_trace(S_chemo_os_raw, S_chemo_efs_raw, S_bg_smr)
# 5. Apply ITT correction for CAR-T group
P1 <- 0.814 # Success proceed to infusion
P2 <- 0.113 # Apheresis only, switch to 50% blina + 50% chemo
P3 <- 1 - P1 - P2 # pre-infusion death (0.073)
trace_tisa_itt <- data.frame(
EFS = P1 * trace_tisa_infused$EFS + P2 * (0.5 * trace_blin$EFS + 0.5 * trace_chemo$EFS) + P3 * 0,
PD = P1 * trace_tisa_infused$PD + P2 * (0.5 * trace_blin$PD + 0.5 * trace_chemo$PD) + P3 * 0,
Dead = P1 * trace_tisa_infused$Dead + P2 * (0.5 * trace_blin$Dead + 0.5 * trace_chemo$Dead) + P3 * 1.0
)
# 6. Apply Half-cycle correction (Trapezoidal Method)
apply_hcc <- function(trace) {
n <- nrow(trace)
trace_hcc <- trace
for (col in c("EFS", "PD")) {
val <- trace[[col]]
val_hcc <- numeric(n)
val_hcc[1] <- 0.5 * val[1]
for (i in 2:n) {
val_hcc[i] <- 0.5 * val[i] + 0.5 * val[i - 1]
}
trace_hcc[[col]] <- val_hcc
}
# Dead is computed as 1 - EFS - PD to ensure consistency
trace_hcc$Dead <- 1 - trace_hcc$EFS - trace_hcc$PD
return(trace_hcc)
}
trace_tisa_hcc <- apply_hcc(trace_tisa_itt)
trace_blin_hcc <- apply_hcc(trace_blin)
trace_chemo_hcc <- apply_hcc(trace_chemo)
# 7. Visualization: Area plots for all three arms
prep_plot_data <- function(trace, arm_name) {
trace %>%
mutate(Cycle = cycles) %>%
pivot_longer(cols = c(EFS, PD, Dead), names_to = "State", values_to = "Proportion") %>%
mutate(
State = factor(State, levels = c("Dead", "PD", "EFS")),
Arm = arm_name
)
}
plot_data <- rbind(
prep_plot_data(trace_tisa_itt, "Tisagenlecleucel (ITT)"),
prep_plot_data(trace_blin, "Blinatumomab"),
prep_plot_data(trace_chemo, "Salvage Chemotherapy")
)
# Render area plot
ggplot(plot_data, aes(x = Cycle / 12, y = Proportion, fill = State)) +
geom_area(alpha = 0.85, color = "white", size = 0.1) +
facet_wrap(~Arm) +
scale_fill_manual(values = c("Dead" = "#E06666", "PD" = "#F6B26B", "EFS" = "#6AA84F")) +
labs(
title = "State Occupancy Over Lifetime Horizon",
subtitle = "Partitioned Survival Model (PSM) with ITT correction & background mortality constraints",
x = "Years in Model",
y = "Proportion of Cohort",
fill = "Health State"
) +
theme_minimal(base_family = "sans") +
theme(
strip.text = element_text(size = 11, face = "bold"),
legend.position = "bottom",
plot.title = element_text(face = "bold", size = 13),
panel.spacing = unit(1.5, "lines")
)
