Breast Cancer: RWE Survival Analysis

1. Environment Setup & Data Ingestion

This analysis evaluates the real-world overall survival (OS) of patients with metastatic breast cancer. This initial section establishes the R environment, loads necessary statistical packages (including WeightIt, survival, and margins), and imports the primary dataset.

library(data.table)
library(readxl)
library(survival)
library(ggsurvfit)
library(survminer)
library(WeightIt)
library(cobalt)
library(ggplot2)
library(car)
library(ResourceSelection)
library(pROC)
library(margins)
library(EValue)
library(ebal)
library(DT)
library(ggpubr)

# Ensure the file path is correct for your local machine
tryCatch({
  dt <- as.data.table(read_excel("C:/Users/samuel.baird/OneDrive - Cardinal Health/Documents/Code_AACR/Data/AKT1-mutatedBreastCancer_Final dataset.xlsx"))
}, error = function(e) {
  stop("File not found. Please check your path.")
})

setnames(dt, old = c("Patient ID", "Primary Race", "Center", "Overall Tumor Grade",
                     "Has this Patient experienced a Loco-Regional Recurrence (LRR)?",
                     "Was a biopsy performed of metastatic site?", "Sequencing Method",
                     "Discordant receptor status with primary",
                     "Interval between sequencing and metastatic diagnosis (months)",
                     "Additional Sequencing Performed", "Fraction Genome Altered",
                     "HER2 Status", "AJCC Stage", "AKT Mutation Status",
                     "Tumor Mutation Count", "Hormone Receptor Status",
                     "ER/PR Receptor Change", "mTOR Inhibitor Overall", 
                     "CDK4/6 Inhibitor Overall", "AKT Inhibitor Overall", 
                     "Fulvestrant Overall", "Aromatase Inhibitor Overall",
                     "Age at Initial Diagnosis", "Age at Metastatic Diagnosis",
                     "Therapeutic clinical trial overall", "Met Site: Bone",
                     "Met Site: Lymph Node", "Met Site: Multiple", 
                     "Met Site: Other", "Met Site: Soft Tissue",
                     "Met Site: Visceral", "Met Site: Bone Only",
                     "Met Site: Lung", "Met Site: Brain", "Met Site: Liver",
                     "Cancer Type Detailed", "Overall Survival (Months)",
                     "Overall Survival Status"),
             new = c("patient_id", "race", "center", "tumor_grade", "lrr",
                     "met_biopsy", "seq_meth", "discordant", "seq_int",
                     "add_seq", "frac_alt", "HER2", "stage", "akt", "TMB",
                     "HR_status", "HR_change", "mTOR", "CDK4/6", 
                     "AKT_inhibitor", "fulvestrant", "aromatase_inhibitor",
                     "prim_age", "met_age", "clinical_trial", "bone_mets",
                     "lymph_mets", "multiple_mets", "other_mets", "soft_mets",
                     "visceral_mets", "boneonly_mets", "lung_mets", "brain_mets",
                     "liver_mets", "cancer_type", "survival_time", "survival_status"))

2. Cohort Definition

To emulate a clinical trial population, the cohort is restricted to patients who are HER2-negative, Hormone Receptor-positive (HR+), and Stage 4.

Patients are stratified into two cohorts: * Group A (Intensive Regimen): Received a sequence of Fulvestrant, an Aromatase Inhibitor, and a targeted therapy (mTOR, CDK4/6, or AKT inhibitor). * Group B (Control): Received alternative standard-of-care regimens.

eligible <- dt[HER2 == "Negative" & (HR_status %in% c("ER+ & PR-", "ER+ & PR+")) & (stage == "4")]

eligible[, fat_flag := (fulvestrant == "Yes" & aromatase_inhibitor == "Yes") & 
           (mTOR == "Yes" | get("CDK4/6") == "Yes" | AKT_inhibitor == "Yes")]

eligible[, therapy_flag := factor(ifelse(fat_flag == TRUE, "Group A", "Group B"), 
                                  levels = c("Group B", "Group A"))]

eligible[, survival_status := ifelse(survival_status == "1:DECEASED", 1, 0)]

3. Exploratory Covariate Testing

Real-world data often contains missingness and sparse sub-categories that can destabilize causal models. This section recategorizes key clinical variables (e.g., grouping age at metastasis, stabilizing tumor grade, and consolidating metastatic sites) to ensure robust statistical balancing later in the pipeline.

eligible[is.na(met_biopsy) | met_biopsy == "NA", met_biopsy := "Unknown"]
eligible[, met_biopsy := as.factor(met_biopsy)]
eligible[is.na(bone_mets) | bone_mets == "NA", bone_mets := "Unknown"]
eligible[is.na(lymph_mets) | lymph_mets == "NA", lymph_mets := "Unknown"]
eligible[, lymph_mets := as.factor(lymph_mets)]
eligible[, multiple_mets := as.factor(multiple_mets)]
eligible[is.na(other_mets) | other_mets == "NA", other_mets := "Unknown"]
eligible[is.na(soft_mets) | soft_mets == "NA", soft_mets := "Unknown"]
eligible[, soft_mets := as.factor(soft_mets)]
eligible[, visceral_mets := as.factor(visceral_mets)]
eligible[, seq_meth := as.factor(seq_meth)]
eligible[is.na(discordant) | discordant == "NA", discordant := "Unknown"]
eligible[, discordant := as.factor(discordant)]

eligible[, TMB := as.numeric(TMB)]
eligible[, TMB := cut(TMB, c(-Inf, 10, Inf), c("low", "high"), right = T)]
eligible[TMB == "NA" | is.na(TMB), TMB := "Unknown"]

eligible[, frac_alt := as.numeric(frac_alt)]
eligible[, frac_alt := cut(frac_alt, c(-Inf, 0.35, Inf), c("low fra", "high fra"), right = T)]
eligible[frac_alt == "NA" | is.na(frac_alt), frac_alt := "Unknown"]

eligible[race %in% c("Black", "Other", "NA"), race := "Non-White"] 
eligible[, cancer_type := ifelse(cancer_type %in% c("Breast Invasive Lobular Carcinoma", "Invasive Breast Carcinoma", "Breast Mixed Ductal and Lobular Carcinoma", "Breast Invasive Carcinoma, NOS"), "Non-Ductal", "Ductal")]
eligible[, liver_mets := ifelse(liver_mets == "NA" | is.na(liver_mets), "Unknown", liver_mets)]
eligible[, brain_mets := ifelse(brain_mets == "NA" | is.na(brain_mets), "Unknown", brain_mets)]
eligible[, lung_mets := ifelse(lung_mets == "NA" | is.na(lung_mets), "Unknown", lung_mets)]
eligible[, met_age := cut(met_age, c(-Inf,55,Inf), labels = c("<=55",">55"), right = TRUE)]

eligible[tumor_grade %in% c("1","2"), tumor_grade := "Grade 1-2"]
eligible[tumor_grade == "3", tumor_grade := "Grade 3"]
eligible[tumor_grade %in% c("NA", "Unk"), tumor_grade := "Unknown"]

eligible[, seq_int := cut(seq_int, c(-Inf, 3, Inf), c("seq int 1-3 months", "seq int 4+ months"), right = T)]
eligible[, add_seq := as.factor(add_seq)]
eligible[, boneonly_mets := as.factor(boneonly_mets)]
eligible[, akt := as.factor(akt)]
eligible[HR_change == "NA" | is.na(HR_change) | HR_change == "Unk", HR_change := "Unknown"]
eligible[, HR_change := as.factor(HR_change)]
eligible[, clinical_trial := as.factor(clinical_trial)]

# 3-Tier Global Hierarchy for Stability across AME and WeightIt
eligible[, met_pattern_stable := fcase(
  liver_mets == "Yes", "Liver Involved",
  boneonly_mets == "Yes", "Bone-Only",
  default = "Other Sites (Non-Liver)" 
)]

eligible[, met_pattern_stable := factor(met_pattern_stable, 
                                        levels = c("Bone-Only", "Other Sites (Non-Liver)", "Liver Involved"))]

eligible[, tumor_grade_stable := factor(ifelse(tumor_grade == "Grade 1-2", "Grade 1-2", "Grade 3/Unknown"),
                                        levels = c("Grade 1-2", "Grade 3/Unknown"))]

# --- 1L Chemotherapy Data Prep ---
# Create the factor and handle the long column name
eligible[, first_chemo := as.factor(`Chemotherapy in first-line treatment`)]

# Define the binary outcome for the GLM (1 = Yes, 0 = No)
# Note: Check if your data uses "Yes" or "1:Yes"
eligible[, first_chemo_binary := ifelse(first_chemo == "Yes", 1, 0)]

# Clean up the old column name to keep the DT tidy
eligible[, `Chemotherapy in first-line treatment` := NULL]

4. Unweighted Survival Analysis

A naive Kaplan-Meier analysis is generated to observe the crude survival outcomes before adjusting for clinical confounding factors.

# 1. Standard KM analysis
surv_object <- Surv(time = eligible[, survival_time], event = eligible[, survival_status])
km_fit <- survfit2(surv_object ~ therapy_flag, data = eligible)

ggsurvplot(km_fit, data = eligible, pval = TRUE, risk.table = TRUE, risk.table.col = "strata",
           surv.median.line = "hv", linetype = "solid", break.time.by = 12,
           title = "Unweighted Overall Survival", 
           legend.title = "Cohort",
           legend.labs = c("Group B (Control)", "Group A (Intensive Regimen)"),
           palette = c("#549ed1", "#f79c59"),
           ggtheme = theme_classic())

# 2. Reverse KM for Median Follow-Up
eligible[, follow_up_status := ifelse(survival_status == 1, 0, 1)]
follow_up_fit <- survfit(Surv(survival_time, follow_up_status) ~ 1, data = eligible)
print(follow_up_fit)

## Call: survfit(formula = Surv(survival_time, follow_up_status) ~ 1, 
##     data = eligible)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 83     52   58.9    45.6    74.5

5. Entropy Balancing (WeightIt)

Due to the lack of randomization in real-world data, treatment assignment is often confounded by patient prognosis (e.g., sicker patients receiving different regimens). We utilize Entropy Balancing to generate patient-level weights, ensuring exact covariate balance across age, tumor grade, AKT status, and metastatic pattern to simulate a randomized average treatment effect (ATE).

weightit_model <- weightit(therapy_flag ~ met_age + tumor_grade_stable + akt + met_pattern_stable,
                           data = eligible, method = "ebal", estimand = "ATE")

summary(weightit_model)

##                   Summary of weights
## 
## - Weight ranges:
## 
##           Min                                 Max
## treated 0.143 |---------------------------| 3.038
## control 0.41     |------------|             1.752
## 
## - Units with the 5 most extreme weights by group:
##                                       
##             22    25    18    16     3
##  treated  1.48  1.48 2.383 2.692 3.038
##             29    19     8    14    17
##  control 1.532 1.677 1.752 1.752 1.752
## 
## - Weight statistics:
## 
##         Coef of Var   MAD Entropy # Zeros
## treated       0.713 0.537   0.221       0
## control       0.349 0.270   0.058       0
## 
## - Effective Sample Sizes:
## 
##            Control Treated
## Unweighted   53.     30.  
## Weighted     47.33   20.12

eligible[, ws := weightit_model$weights]

final_labels_weight <- c(
  "met_age_>55" = "Age > 55 at Diagnosis",
  "akt_Wildtype" = "AKT Wildtype",
  "tumor_grade_stable_Grade 3/Unknown" = "Tumor Grade: 3 or Unknown",
  "met_pattern_stable_Bone-Only" = "Bone-Only",
  "met_pattern_stable_Liver Involved" = "Liver Involved",
  "met_pattern_stable_Other Sites (Non-Liver)" = "Other Sites (Non-Liver)"
)

love.plot(weightit_model, thresholds = 0.1, abs = TRUE, var.order = "unadjusted", 
          line = TRUE, shapes = c(16, 17), limits = c(0, 0.6), colors = c("#f79c59", "#549ed1"), 
          var.names = final_labels_weight,
          theme = theme_classic() + theme(plot.title = element_text(face = "bold")))

6. Adjusted (Weighted) Survival Analysis

Using the entropy weights, an adjusted Cox Proportional Hazards model and corresponding Kaplan-Meier curves are generated. This represents the primary outcome of the study, estimating the isolated effect of Group A therapy on overall survival. We also compute an E-value to quantify the minimum strength of an unmeasured confounder required to nullify these results.

km_fit_weighted <- survfit2(Surv(survival_time, survival_status) ~ therapy_flag,
                            data = eligible, weights = ws)

cox_model_weighted <- coxph(Surv(survival_time, survival_status) ~ therapy_flag, 
                            data = eligible, weights = ws, robust = TRUE)
summary(cox_model_weighted)

## Call:
## coxph(formula = Surv(survival_time, survival_status) ~ therapy_flag, 
##     data = eligible, weights = ws, robust = TRUE)
## 
##   n= 83, number of events= 31 
## 
##                        coef exp(coef) se(coef) robust se      z Pr(>|z|)  
## therapy_flagGroup A -0.9458    0.3884   0.4264    0.4588 -2.062   0.0392 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A    0.3884      2.575     0.158    0.9544
## 
## Concordance= 0.628  (se = 0.046 )
## Likelihood ratio test= 5.51  on 1 df,   p=0.02
## Wald test            = 4.25  on 1 df,   p=0.04
## Score (logrank) test = 5.27  on 1 df,   p=0.02,   Robust = 4.99  p=0.03
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

# Test Proportional Hazards Assumption
ph_test <- cox.zph(cox_model_weighted, transform = "km", terms = TRUE)
print(ph_test)

##              chisq df    p
## therapy_flag  1.01  1 0.32
## GLOBAL        1.01  1 0.32

ggsurvplot(km_fit_weighted, data = eligible, pval = TRUE, risk.table = TRUE, risk.table.col = "strata",
           surv.median.line = "hv", linetype = "solid", break.time.by = 12,
           title = "Adjusted Overall Survival (Entropy Balanced)", 
           legend.title = "Cohort",
           legend.labs = c("Group B (Control)", "Group A (Intensive Regimen)"),
           palette = c("#549ed1", "#f79c59"),
           ggtheme = theme_classic())

evalues.HR(est = 0.39, lo = 0.16, hi = 0.95, true = 1, rare = FALSE)

##              point     lower     upper
## RR       0.5240857 0.2941386 0.9650718
## E-values 3.2244070        NA 1.2298470

7. Univariable Cox Models

To explore the independent prognostic value of clinical covariates, we run univariable Cox models on both the raw and the entropy-balanced cohorts.

covariates_filtered <- c("met_age", "tumor_grade_stable", "akt", "met_pattern_stable")

full_label_map <- c(
  "met_age>55" = "Age > 55 at Diagnosis",
  "aktWildtype" = "AKT Wildtype",
  "tumor_grade_stableGrade 3/Unknown" = "Tumor Grade: 3 or Unknown",
  "met_pattern_stableLiver Involved" = "Liver Involved",
  "met_pattern_stableOther Sites (Non-Liver)" = "Other Sites (Non-Liver)",
  "met_pattern_stableBone-Only" = "Bone-Only"
)

# --- UNWEIGHTED UNIVARIATE ---
uni_cox_unweighted <- lapply(covariates_filtered, function(cov) {
  formula <- as.formula(paste("Surv(survival_time, survival_status) ~", cov))
  fit <- coxph(formula, data = eligible) 
  return(summary(fit))
})
names(uni_cox_unweighted) <- covariates_filtered

uni_table_unweighted <- rbindlist(lapply(names(uni_cox_unweighted), function(cov_name) {
  s <- uni_cox_unweighted[[cov_name]]
  res <- as.data.table(s$conf.int, keep.rownames = "Raw_Level")
  res[, p_val_numeric := s$coefficients[, 5]] # index 5 for unweighted standard
  return(res)
}))

uni_table_unweighted[, Variable := ifelse(Raw_Level %in% names(full_label_map), full_label_map[Raw_Level], Raw_Level)]

datatable(uni_table_unweighted[, .(Variable, HR = round(`exp(coef)`, 2), Lower = round(`lower .95`, 2), Upper = round(`upper .95`, 2), P = round(p_val_numeric, 3))], 
          rownames = FALSE, options = list(dom = 't'), caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 1: Unweighted Univariate Analysis'))

# --- WEIGHTED UNIVARIATE ---
uni_cox_weighted <- lapply(covariates_filtered, function(cov) {
  formula <- as.formula(paste("Surv(survival_time, survival_status) ~", cov))
  fit <- coxph(formula, data = eligible, weights = ws, robust = TRUE) 
  return(summary(fit))
})
names(uni_cox_weighted) <- covariates_filtered

uni_table_weighted <- rbindlist(lapply(names(uni_cox_weighted), function(cov_name) {
  s <- uni_cox_weighted[[cov_name]]
  res <- as.data.table(s$conf.int, keep.rownames = "Raw_Level")
  res[, p_val_numeric := s$coefficients[, 6]] # index 6 for robust models
  return(res)
}))

uni_table_weighted[, Variable := ifelse(Raw_Level %in% names(full_label_map), full_label_map[Raw_Level], Raw_Level)]

datatable(uni_table_weighted[, .(Variable, HR = round(`exp(coef)`, 2), Lower = round(`lower .95`, 2), Upper = round(`upper .95`, 2), P = round(p_val_numeric, 3))], 
          rownames = FALSE, options = list(dom = 't'), caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 2: Weighted Univariate Analysis'))

8. Landmark OS Analysis

To assess the durability of the treatment effect over time and account for potential immortal time bias in early stages, we calculate Hazard Ratios at fixed annual landmarks (1 to 6 years) post-metastatic diagnosis.

landmarks <- seq(12, 72, by = 12)

# --- UNWEIGHTED LANDMARK ---
landmark_unweighted <- rbindlist(lapply(landmarks, function(t) {
  d_land <- eligible[survival_time >= t][, st_adj := survival_time - t]
  if(nrow(d_land) < 5 || length(unique(d_land$therapy_flag)) < 2 || sum(d_land$survival_status) < 1) return(NULL)
  
  fit <- coxph(Surv(st_adj, survival_status) ~ therapy_flag, data = d_land)
  s <- summary(fit)
  
  data.table(Time = paste(t, "Months"), N = nrow(d_land), Events = sum(d_land$survival_status),
             HR = round(s$conf.int[1, 1], 2), Lower = round(s$conf.int[1, 3], 2), Upper = round(s$conf.int[1, 4], 2), P = round(s$coefficients[1, 5], 3))
}))

datatable(landmark_unweighted, rownames = FALSE, options = list(dom = 't'), 
          caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 3: Unweighted Landmark Analysis'))

# --- WEIGHTED LANDMARK ---
landmark_weighted <- rbindlist(lapply(landmarks, function(t) {
  d_land <- eligible[survival_time >= t][, st_adj := survival_time - t]
  if(nrow(d_land) < 5 || length(unique(d_land$therapy_flag)) < 2 || sum(d_land$survival_status) < 1) return(NULL)
  
  fit <- coxph(Surv(st_adj, survival_status) ~ therapy_flag, data = d_land, weights = d_land$ws, robust = TRUE)
  s <- summary(fit)
  
  data.table(Time = paste(t, "Months"), N = nrow(d_land), Events = sum(d_land$survival_status),
             HR = round(s$conf.int[1, 1], 2), Lower = round(s$conf.int[1, 3], 2), Upper = round(s$conf.int[1, 4], 2), P = round(s$coefficients[1, 6], 3))
}))

datatable(landmark_weighted, rownames = FALSE, options = list(dom = 't'), 
          caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 4: Weighted Landmark Analysis'))

9. Clinical Drivers of Selection (Diagnostics)

Understanding why a physician selected a specific therapy in the real world is critical. This section uses logistic regression models and calculates Average Marginal Effects (AME) to quantify how baseline clinical characteristics influenced the probability of being assigned to the Control arm versus the Intensive Regimen arm, as well as the probability of initiating first-line chemotherapy (therapeutic diversion). Goodness-of-fit is assessed via Hosmer-Lemeshow and AUC.

# ==============================================================================
# GROUP B (CONTROL) DIAGNOSTICS & AME
# ==============================================================================
eligible[, therapy_flag_B := ifelse(therapy_flag == "Group B", 1, 0)]

# Set factors to Low-Risk baseline for model
eligible[, met_age_B := factor(met_age, levels = c("<=55", ">55"))]
eligible[, tumor_grade_B := factor(tumor_grade_stable, levels = c("Grade 1-2", "Grade 3/Unknown"))]
eligible[, akt_B := factor(akt, levels = c("Mutant", "Wildtype"))]
eligible[, met_pattern_B := factor(met_pattern_stable, levels = c("Bone-Only", "Liver Involved", "Other Sites (Non-Liver)"))]

glm_B <- glm(therapy_flag_B ~ met_age_B + tumor_grade_B + akt_B + met_pattern_B, data = eligible, family = binomial())

# --- Model Validation ---
y_B <- glm_B$y 
yhat_B <- fitted(glm_B)

# Calibration
hl_B <- hoslem.test(y_B, yhat_B, g = 6)
# Discrimination
roc_B <- roc(y_B, yhat_B, quiet = TRUE)
auc_val_b <- round(auc(roc_B), 3)

cat("--- Group B Assignment Prediction ---\n")

## --- Group B Assignment Prediction ---

# --- AME Generation ---
ame_B <- as.data.table(summary(margins(glm_B)))
labels_B <- c("met_pattern_BLiver Involved" = "Liver Involved", "met_pattern_BOther Sites (Non-Liver)" = "Other Sites (Non-Liver)", 
              "tumor_grade_BGrade 3/Unknown" = "Tumor Grade: 3 or Unknown", "met_age_B>55" = "Age > 55 at Diagnosis", "akt_BWildtype" = "AKT Wildtype")
ame_B[, factor := ifelse(factor %in% names(labels_B), labels_B[factor], factor)]

ame_table_b_display <- ame_B[, .(Variable = factor, "Marginal Effect" = round(AME, 3), Lower = round(lower, 3), Upper = round(upper, 3), "P-Value" = round(p, 3))]

datatable(ame_table_b_display, rownames = FALSE, options = list(dom = 't'), 
          caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 5A: Average Marginal Effects for Group B Assignment'))

ggplot(ame_B, aes(x = AME, y = reorder(factor, AME))) +
  geom_point(size = 3, color = "black") + geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") + scale_x_continuous(labels = scales::percent) +
  labs(
    title = "Clinical Drivers of Group B (Control) Selection",
    x = "Change in Probability of Receiving Group B",
    y = NULL,
    caption = paste0(
      "Note: Values represent Average Marginal Effects (AME). A positive % indicates a higher absolute probability\n",
      "of assignment to the Control arm relative to the baseline, holding other factors constant.\n",
      "Baseline references reflect 'low-risk' criteria: Bone-Only, AKT Mutant, Age <= 55, and Grade 1-2."
    )
  ) +
  theme_pubr() + theme(plot.caption = element_text(hjust = 0, size = 9, face = "italic", color = "grey30"))

# ==============================================================================
# GROUP A (Intensive Regimen) DIAGNOSTICS & AME
# ==============================================================================
eligible[, therapy_flag_A := ifelse(therapy_flag == "Group A", 1, 0)]

# Set factors to High-Risk baseline for model
eligible[, met_age_A := factor(met_age, levels = c(">55", "<=55"))]
eligible[, tumor_grade_A := factor(tumor_grade_stable, levels = c("Grade 3/Unknown", "Grade 1-2"))]
eligible[, akt_A := factor(akt, levels = c("Wildtype", "Mutant"))]
eligible[, met_pattern_A := factor(met_pattern_stable, levels = c("Liver Involved", "Bone-Only", "Other Sites (Non-Liver)"))]

glm_A <- glm(therapy_flag_A ~ met_age_A + tumor_grade_A + akt_A + met_pattern_A, data = eligible, family = binomial())

# --- Model Validation ---
y_A <- glm_A$y 
yhat_A <- fitted(glm_A)

# Calibration
hl_A <- hoslem.test(y_A, yhat_A, g = 6)
# Discrimination
roc_A <- roc(y_A, yhat_A, quiet = TRUE)
auc_val_a <- round(auc(roc_A), 3)

cat("--- Group A Assignment Prediction ---\n")

## --- Group A Assignment Prediction ---

# --- AME Generation ---
ame_A <- as.data.table(summary(margins(glm_A)))
labels_A <- c("met_age_A<=55" = "Age <= 55 at Diagnosis", "tumor_grade_AGrade 1-2" = "Tumor Grade: 1-2", "akt_AMutant" = "AKT Mutant",
              "met_pattern_ABone-Only" = "Bone-Only", "met_pattern_AOther Sites (Non-Liver)" = "Other Sites (Non-Liver)")
ame_A[, factor := ifelse(factor %in% names(labels_A), labels_A[factor], factor)]

ame_table_a_display <- ame_A[, .(Variable = factor, "Marginal Effect" = round(AME, 3), Lower = round(lower, 3), Upper = round(upper, 3), "P-Value" = round(p, 3))]

datatable(ame_table_a_display, rownames = FALSE, options = list(dom = 't'), 
          caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 5B: Average Marginal Effects for Group A Assignment'))

ggplot(ame_A, aes(x = AME, y = reorder(factor, AME))) +
  geom_point(size = 3, color = "black") + geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") + scale_x_continuous(labels = scales::percent, limits = c(-0.6, 0.6)) +
  labs(
    title = "Clinical Drivers of Intensive Regimen Therapy (Group A) Selection",
    x = "Change in Probability of Receiving Group A",
    y = NULL,
    caption = paste0(
      "Note: Values represent Average Marginal Effects (AME). A positive % indicates a higher absolute probability\n",
      "of assignment to the Intensive Regimen arm relative to the baseline, holding other factors constant.\n",
      "Baseline references reflect 'high-risk' criteria: Liver Involved, AKT Wildtype, Age > 55, and Grade 3/Unknown."
    )
  ) +
  theme_pubr() + theme(plot.caption = element_text(hjust = 0, size = 9, face = "italic", color = "grey30"))

Therapeutic Diversion (1L Chemotherapy)

To further contextualize the clinical decision-making process, we evaluate the baseline clinical features driving the initiation of first-line systemic chemotherapy. This analysis aims to confirm whether specific high-risk presentations systematically diverted patients away from targeted standard-of-care options.

# ==============================================================================
# 1L CHEMOTHERAPY SELECTION (THERAPEUTIC DIVERSION)
# ==============================================================================

# 1. Define the Outcome within the chunk to keep Section 3 unchanged
eligible[, first_chemo_binary := ifelse(first_chemo == "Yes", 1, 0)]

# 2. Fit the 1L Chemo Selection Model 
# We use the same 'Stable' variables and references as Group B/A for consistency
glm_chemo_1L <- glm(
  first_chemo_binary ~ met_age + tumor_grade_stable + akt + met_pattern_stable,
  data = eligible,
  family = binomial()
)

# 3. Model Diagnostics (Discrimination & Calibration)
y_chemo <- glm_chemo_1L$y
yhat_chemo <- fitted(glm_chemo_1L)
auc_chemo <- round(auc(roc(y_chemo, yhat_chemo, quiet = TRUE)), 3)
hl_chemo <- hoslem.test(y_chemo, yhat_chemo, g = 6)

cat("--- 1L Chemotherapy Prediction Diagnostics ---\n")

## --- 1L Chemotherapy Prediction Diagnostics ---

cat("AUC: ", auc_chemo, "\n")

## AUC:  0.764

cat("Hosmer-Lemeshow p-value: ", round(hl_chemo$p.value, 3), "\n\n")

## Hosmer-Lemeshow p-value:  0.876

# 4. Generate AME Table & Labels
ame_chemo <- as.data.table(summary(margins(glm_chemo_1L)))

labels_chemo <- c(
  "met_pattern_stableLiver Involved" = "Liver Involved",
  "met_pattern_stableOther Sites (Non-Liver)" = "Other Sites (Non-Liver)", 
  "tumor_grade_stableGrade 3/Unknown" = "Tumor Grade: 3 or Unknown",
  "met_age>55" = "Age > 55 at Diagnosis",
  "aktWildtype" = "AKT Wildtype"
)

ame_chemo[, factor := ifelse(factor %in% names(labels_chemo), labels_chemo[factor], factor)]

# 5. Display the Data Table
ame_table_chemo_display <- ame_chemo[, .(Variable = factor, "Marginal Effect" = round(AME, 3), Lower = round(lower, 3), Upper = round(upper, 3), "P-Value" = round(p, 3))]

datatable(ame_table_chemo_display, rownames = FALSE, options = list(dom = 't'), 
          caption = htmltools::tags$caption(style = 'caption-side: top; text-align: left; font-weight: bold; color: black;', 'Table 5C: Average Marginal Effects for 1L Chemotherapy Selection'))

# 6. Final AME Plot for 1L Chemotherapy (No Reference Row)
ggplot(ame_chemo, aes(x = AME, y = reorder(factor, AME))) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray60") +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.2, color = "black") +
  geom_point(size = 3.5, color = "black") +
  scale_x_continuous(labels = scales::percent_format(accuracy = 1), limits = c(-0.6, 0.8)) +
  labs(
    title = "Clinical Drivers of 1L Chemotherapy Selection",
    x = "Change in Probability of Receiving 1L Chemotherapy",
    y = NULL,
    caption = paste0(
      "Note: Values represent Average Marginal Effects (AME). A positive % indicates a higher absolute probability\n",
      "of receiving 1L chemotherapy relative to the baseline (Bone-Only, AKT Mutant, Age <= 55, Grade 1-2)."
    )
  ) +
  theme_pubr() + 
  theme(plot.caption = element_text(hjust = 0, size = 9, face = "italic", color = "grey30"),
        plot.title = element_text(face = "bold"))

10. Sensitivity Analyses

To validate the robustness of the primary entropy balancing model, we estimate the ATE using Covariate Balancing Propensity Scores (CBPS) and also conduct a sensitivity check to assess if the sequencing interval (‘seq_int’) acts as an unmeasured confounder influencing the treatment effect.

# ==============================================================================
# Sensitivity 1: CBPS Weights
# ==============================================================================
cbps_model <- weightit(therapy_flag ~ met_age + tumor_grade_stable + akt + met_pattern_stable,
                       data = eligible, 
                       method = "cbps", 
                       estimand = "ATE")

eligible[, ws_cbps := cbps_model$weights]

cox_cbps <- coxph(Surv(survival_time, survival_status) ~ therapy_flag, 
                  data = eligible, 
                  weights = ws_cbps, 
                  robust = TRUE)

summary(cox_cbps)

## Call:
## coxph(formula = Surv(survival_time, survival_status) ~ therapy_flag, 
##     data = eligible, weights = ws_cbps, robust = TRUE)
## 
##   n= 83, number of events= 31 
## 
##                        coef exp(coef) se(coef) robust se      z Pr(>|z|)  
## therapy_flagGroup A -1.0226    0.3597   0.2906    0.4415 -2.316   0.0206 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                     exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A    0.3597       2.78    0.1514    0.8545
## 
## Concordance= 0.654  (se = 0.051 )
## Likelihood ratio test= 12.83  on 1 df,   p=3e-04
## Wald test            = 5.36  on 1 df,   p=0.02
## Score (logrank) test = 13.42  on 1 df,   p=2e-04,   Robust = 6.12  p=0.01
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

km_fit_cbps <- survfit2(Surv(survival_time, survival_status) ~ therapy_flag,
                        data = eligible, weights = ws_cbps)

ggsurvplot(km_fit_cbps, data = eligible, pval = TRUE, risk.table = TRUE, risk.table.col = "strata",
           surv.median.line = "hv", linetype = "solid", break.time.by = 12,
           title = "Sensitivity OS Analysis (CBPS Weighted)", 
           legend.title = "Cohort",
           legend.labs = c("Group B (Control)", "Group A (Intensive Regimen)"),
           palette = c("#549ed1", "#f79c59"),
           ggtheme = theme_classic())

# ==============================================================================
# Sensitivity 2: Adjusting for Sequencing Interval (seq_int)
# ==============================================================================
cox_sensitivity <- coxph(
  Surv(survival_time, survival_status) ~ therapy_flag + seq_int, 
  data = eligible, 
  weights = ws, 
  robust = TRUE
)

summary(cox_sensitivity)

## Call:
## coxph(formula = Surv(survival_time, survival_status) ~ therapy_flag + 
##     seq_int, data = eligible, weights = ws, robust = TRUE)
## 
##   n= 83, number of events= 31 
## 
##                             coef exp(coef) se(coef) robust se      z Pr(>|z|)
## therapy_flagGroup A      -0.9775    0.3763   0.4279    0.4680 -2.089 0.036734
## seq_intseq int 4+ months -1.4559    0.2332   0.4486    0.3795 -3.836 0.000125
##                             
## therapy_flagGroup A      *  
## seq_intseq int 4+ months ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                          exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A         0.3763      2.658    0.1504    0.9415
## seq_intseq int 4+ months    0.2332      4.288    0.1108    0.4906
## 
## Concordance= 0.713  (se = 0.041 )
## Likelihood ratio test= 15.5  on 2 df,   p=4e-04
## Wald test            = 20.75  on 2 df,   p=3e-05
## Score (logrank) test = 16.68  on 2 df,   p=2e-04,   Robust = 10.8  p=0.005
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

11. Oncology Summit Provider Analysis (ARS Results)

In this section, we analyze clinician response data (ARS) from the Oncology Summits to investigate the subjective rationale behind treatment selection. This provides a “real-world consensus” layer to complement the survival analysis.

#### Prepare the Data ####
file_path_day1 <- "~./Code_AACR/Data/2026_Feb27_Onc_Summit_Day 1_ARSResults.xlsx"
dt2 <- as.data.table(read_xlsx(file_path_day1, sheet = 3, skip = 13, na = "-"))
dt2 <- dt2[-c(1,.N)] #Remove first and last row
dt2 <- dt2[!apply(dt2, 1, function(x) any(grepl("I do not manage patients with breast cancer", x, fixed = TRUE)))] #Remove non-providers

file_path_day2 <- "~./Code_AACR/Data/2026_Mar14_Onc_Summit_Day 2_ARSResults.xlsm"
dt3 <- as.data.table(read_xlsx(file_path_day2, sheet = 4, skip = 12, na = "-"))
dt3 <- dt3[-c(1,.N)] #Remove first and last row
dt3 <- dt3[!apply(dt3, 1, function(x) any(grepl("I do not manage patients with breast cancer", x, fixed = TRUE)))] #Remove non-providers

#### Untie Formula ####
untie <- function(question, data) {
  # 1. Hard-calculate the denominator
  valid_rows <- data[!is.na(get(question)) & get(question) != "", ]
  denom_value <- uniqueN(valid_rows, by = NULL) 
  
  # 2. Parse the strings
  a <- valid_rows[, .(question_clean = trimws(unlist(strsplit(as.character(get(question)), "\r?\n")))), 
                  by = .I] 
  
  # 3. Aggregate counts
  res <- a[, .(N = .N), by = .(Category = question_clean)]
  
  # 4. Final Math
  res[, prop := N / denom_value]
  return(res[order(-N)])
}

#### 4. Dorinda Choice (Q4 & Q22) ####
q_dor_choice <- rbind(dt2[!is.na(Q4) & Q4 != "I do not manage patients with breast cancer", .(Response = Q4)], 
                      dt3[!is.na(Q22) & Q22 != "I do not manage patients with breast cancer", .(Response = Q22)])
q_dor_choice_dis <- q_dor_choice[, .(N = .N), by = Response][, prop := N / sum(N)]
q_dor_choice_dis[grep("Single-agent endocrine", Response), Response := "Single-agent endocrine therapy"]

ggplot(q_dor_choice_dis, aes(x = reorder(Response, -N), y = N)) +
  geom_col(fill = "steelblue", color = "black") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  geom_text(aes(label = paste0(N, " (", round(prop * 100, 1), "%)")), vjust = -0.5) +
  labs(title = "Dorinda Patient Case: Preferred 1L Therapy (Combined)",
       x = "Choice", y = "Count",
       caption = paste0("Total N = ", sum(q_dor_choice_dis$N)))

#### 5. Dorinda Influencing Factors (Q5 & Q23) ####
# 1. Run untie
q5_res <- untie("Q5", dt2)
q23_res <- untie("Q23", dt3)

# 2. Denominator math
total_n_dor <- round(q5_res$N[1] / q5_res$prop[1]) + round(q23_res$N[1] / q23_res$prop[1])

# 3. Standardize and MERGE the two Liver categories
# This treats "Visceral crisis" and "Liver metastasis" as the same factor
q5_res[grepl("Visceral|metastasis", Category, ignore.case = T), Category := "Visceral crisis/Liver metastasis"]
q23_res[grepl("Visceral|metastasis", Category, ignore.case = T), Category := "Visceral crisis/Liver metastasis"]

# Standardize other discrepancies
q5_res[grepl("novo", Category), Category := "Stage/prior systemic therapy"]
q23_res[grepl("novo", Category), Category := "Stage/prior systemic therapy"]
q5_res[grepl("Age", Category),  Category := "Age (55-year-old)"]
q23_res[grepl("Age", Category),  Category := "Age (55-year-old)"]
q5_res[grepl("Menopausal", Category), Category := "Postmenopausal status"]
q23_res[grepl("Menopausal", Category), Category := "Postmenopausal status"]

# 4. Merge and SUM (This collapses the repeats into single bars)
q_dor_inf_dis <- rbind(q5_res[, .(Category, N)], q23_res[, .(Category, N)])
q_dor_inf_dis <- q_dor_inf_dis[, .(N = sum(N)), by = Category]
q_dor_inf_dis[, prop := N / total_n_dor]

# 5. LOCK THE ORDER (Forces Visceral/Liver to the top)
q_dor_inf_dis <- q_dor_inf_dis[order(-N)]
q_dor_inf_dis$Category <- factor(q_dor_inf_dis$Category, levels = q_dor_inf_dis$Category)

# 6. Plot
ggplot(q_dor_inf_dis, aes(x = Category, y = N)) +
  geom_col(fill = "steelblue", color = "black") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  geom_text(aes(label = paste0(N, " (", round(prop * 100, 1), "%)")), vjust = -0.5) +
  labs(title = "Factors Influencing Treatment for Dorinda (Combined)",
       x = "Reason", y = "Count",
       caption = paste0("Total N = ", total_n_dor, " respondents; multi-select."))

# 1. Run untie
q58_res <- untie("Q58", dt2)
q67_res <- untie("Q67", dt3)

# 2. Denominator math
total_n_reasons <- round(q58_res$N[1]/q58_res$prop[1]) + round(q67_res$N[1]/q67_res$prop[1])

# 3. Standardize Strings (Merging Visceral Crisis and Cost/Access)
# This ensures "Cost/access" and "Cost to patient/access" are treated as the same
q58_res[grepl("Visceral crisis", Category), Category := "Presence of visceral crisis"]
q67_res[grepl("Visceral crisis", Category), Category := "Presence of visceral crisis"]

q58_res[grepl("high visceral tumor burden", Category), Category := "High visceral tumor burden"]
q67_res[grepl("high visceral tumor burden", Category), Category := "High visceral tumor burden"]

# MERGE COST/ACCESS BARRIERS
q58_res[grepl("Cost.*access", Category, ignore.case = T), Category := "Cost/Access barriers"]
q67_res[grepl("Cost.*access", Category, ignore.case = T), Category := "Cost/Access barriers"]

# 4. Merge and SUM (Collapsing the rows)
q_reasons_dis <- rbind(q58_res[, .(Category, N)], q67_res[, .(Category, N)])
q_reasons_dis <- q_reasons_dis[, .(N = sum(N)), by = Category]
q_reasons_dis[, prop := N / total_n_reasons]

# 5. LOCK THE ORDER (Highest N on the left)
q_reasons_dis <- q_reasons_dis[order(-N)]
q_reasons_dis$Category <- factor(q_reasons_dis$Category, levels = q_reasons_dis$Category)

# 6. Plot
# 6. Plot (Ensure you run this whole block)
ggplot(q_reasons_dis, aes(x = Category, y = N)) +
  geom_col(fill = "steelblue", color = "black") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1)) +
  geom_text(aes(label = paste0(N, " (", round(prop * 100, 1), "%)")), vjust = -0.5) +
  labs(title = "Reasons for 1L Chemotherapy vs Targeted Therapy (Combined)",
       x = "Reason", y = "Count",
       caption = paste0("Total N = ", total_n_reasons, " respondents; multi-select."))

12. Manuscript Drafts

Pitch

This study identifies a significant “effectiveness gap” in metastatic breast cancer by demonstrating that the 34-month survival advantage associated with an intensive endocrine/targeted regimen (characterized by the use of an aromatase inhibitor, fulvestrant, and at least one targeted therapy [mTOR, CDK4/6, or AKT1 inhibitor], and potentially tamoxifen or ovarian suppression) (89.9 vs. 55.9 months compared to conventional endocrine/targeted therapy) is not merely an artifact of patient selection. Utilizing a mature real-world cohort (N=83) from various institutions, with nearly five years of follow-up, we identified a distinct clinical trend where high-risk visceral patients were differentially allocated to alternative care—including first-line chemotherapy. Specifically, the presence of liver metastases drove a nearly 44% absolute increase in the probability of initiating chemotherapy compared to patients with bone-only disease. However, through rigorous Entropy Balancing of baseline characteristics, we confirmed that even after neutralizing these initial clinical differences, the intensive endocrine/targeted regimen provided a robust 61% reduction in the risk of death (\(aHR: 0.39; p=0.039\)). While the relatively small sample size (\(N=83\), with an effective sample size of 20.1 for the intensive regimen group after weighting) necessitates cautious interpretation, these findings strongly suggest that the current clinical practice of prioritizing chemotherapy over this intensive regimen for high-risk visceral patients is associated with a substantial loss of survival years. This highlights the urgent need to re-evaluate current treatment paradigms, suggesting that this intensive endocrine/targeted regimen should be considered for high-risk visceral patients, irrespective of initial metastatic burden.

Statistical Methods

To compare clinical outcomes across treatment groups, we first defined overall survival (OS) as the primary endpoint. Survival distributions were estimated using the Kaplan-Meier methodology, and the association between treatment and survival was evaluated via Cox proportional hazards regression. To address the lack of randomization and potential selection bias inherent to this real-world cohort, Entropy Balancing¹ was employed to calculate patient-level weights. This method has been increasingly applied in medical research to address confounding in observational studies by achieving exact covariate balance². The covariates included in the balancing model—age, tumor grade, AKT mutation status, and metastatic presentation—were selected a priori based on their established prognostic significance and their likelihood to influence clinical treatment selection. By ensuring balance across these factors, this approach aims to estimate the Average Treatment Effect (ATE), simulating a scenario where the treatment assignment is independent of observed baseline characteristics.

To mitigate immortal time bias, all treatment assignments were required to be initiated subsequent to the baseline assessment (defined as the date of metastatic diagnosis). The intensive regimen group (Group A) was defined by the cumulative exposure to an aromatase inhibitor, fulvestrant, and at least one targeted agent (mTOR, CDK4/6, or AKT inhibitor), representing a high-intensity strategy rather than a strictly concurrent triple-therapy protocol, while the control group (Group B) received conventional endocrine/targeted therapy that did not comprise all three agent types.

Baseline characteristics were compared between groups using Fisher’s exact test, Chi-square test, or Wilcoxon rank-sum test. Statistical analyses were performed using R (v. 4.5.2) utilizing the WeightIt, survival, and margins packages. All statistical tests were two-sided, and a \(p\)-value < 0.05 was considered statistically significant. The proportional hazards assumption was verified for all models using the Global Schoenfeld residual test. Following weighting, an adjusted Cox model with Huber-White robust (sandwich) standard errors³ was utilized to estimate the primary hazard ratio (HR). To quantify the potential impact of unmeasured confounding, an E-value⁴ was computed for the primary ATE. Furthermore, sensitivity analyses were conducted utilizing Covariate Balancing Propensity Scores (CBPS)⁵ and adjusting for sequencing intervals. Finally, to evaluate clinical drivers of treatment selection and identify potential therapeutic diversion, logistic regression models with Average Marginal Effects (AME) were generated. The discrimination and calibration of these models were evaluated using the Area Under the Receiver Characteristic Curve (AUC) and the Hosmer-Lemeshow goodness-of-fit test⁶, respectively.

To validate the clinical decision-making patterns observed in the retrospective cohort, a supplemental cross-sectional analysis of provider preferences was conducted using anonymized audience response system (ARS) data from two recent oncology summits. Descriptive statistics were utilized to evaluate treatment preferences and influencing clinical factors in response to a standardized high-risk patient vignette, as well as the general rationales for prioritizing first-line chemotherapy over targeted standard-of-care options. This analysis aimed to subjectively contextualize the objective predictive modeling of therapeutic diversion.

Results

From an initial dataset of 185 patients, 83 met final eligibility criteria (Control Group, N=53; Intensive Regimen Group, N=30). The median follow-up for the study was 58.9 months (95% CI: 45.6–74.5), with a maximum observation period of 111.7 months for the control group and 166.6 months for the intensive regimen group; 31 deaths were observed during the study period. Following entropy balancing, robust covariate balance was achieved (ASMD < 0.1 across all factors), and the pseudo-population maintained a stable effective sample size (ESS) of 20.1 in the intensive group and 47.3 in the control group.

In the unweighted analysis, the intensive regimen group demonstrated improved median overall survival compared to the control group (121.6 [95% CI: 77.7, NE] vs. 55.9 [95% CI: 46.4, NE] months; HR: 0.35 [95% CI 0.15-81]; \(p=0.013\)). The proportional hazards assumption for the unweighted model was satisfied (Global \(p=0.39\)). In the primary entropy-balanced model, the survival benefit remained robust, with a median OS of 89.9 months for the intensive cohort versus 55.9 months for the control cohort (Adjusted HR: 0.39; 95% CI: 0.16–0.95; \(p=0.039\)). The proportional hazards assumption was satisfied for the weighted model (Global \(p=0.32\)).

Robustness was validated in sensitivity analyses (CBPS Adjusted HR: 0.36; 95% CI, 0.15–0.85; \(p=0.021\); E-value = 3.22), indicating that an unmeasured-confounder-outcome association would need a hazard ratio of 3.22 to nullify the observed effect. Similarly, multivariable adjustment for the interval between sequencing and metastatic diagnosis confirmed that the survival benefit remained independent of sequencing timing (Adjusted HR: 0.38; 95% CI, 0.15–0.94; \(p = 0.037\)). Univariate analyses demonstrated that liver involvement was the sole significant independent predictor of mortality in the balanced cohort (Adjusted HR: 2.59; 95% CI, 1.04–6.48; \(p=0.041\)), and landmark OS analyses confirmed a durable treatment effect at the 12-month post-diagnosis landmark (Adjusted HR: 0.39; 95% CI: 0.16–0.95; \(p=0.039\)).

Treatment assignment was significantly influenced by metastatic site, revealing systematic patterns of selection bias between the two cohorts. Patients with liver involvement had a significantly higher absolute probability of receiving the control therapy (+25.2% AME; 95% CI: 3.8%–46.5%; \(p = 0.021\)). Conversely, the intensive regimen was disproportionately selected for patients with bone-only disease (+25.2% AME; 95% CI: 3.8%–46.5%; \(p = 0.021\)) or other non-liver metastatic sites (+34.9% AME; 95% CI: 11.9%–58.0%; \(p = 0.003\)) when compared to those with liver involvement. These models demonstrated adequate discrimination (AUC = 0.753) and calibration (Hosmer-Lemeshow \(p \ge 0.534\) for both).

To further contextualize these selection patterns and identify why certain high-risk patients may have been diverted from endocrine-based strategies, we evaluated the drivers of first-line (1L) chemotherapy initiation. The selection of 1L therapy followed a distinct visceral risk gradient based on the site of metastatic disease. This was captured by a prediction model with adequate discrimination (\(AUC = 0.764\)) and calibration (Hosmer-Lemeshow \(p = 0.876\)), which identified liver involvement as the dominant driver of 1L chemotherapy initiation. Patients with liver metastases experienced a 43.6% higher absolute probability of receiving 1L chemotherapy compared to the those with bone-only disease (95% CI: 16.9%–70.3%; \(p = 0.001\)), while even non-liver visceral involvement (Other Sites) triggered a significant \(30.2\%\) absolute increase relative to bone-only disease (95% CI: 10.5%–49.9%; \(p = 0.003\)).

To validate these real-world selection patterns, a cross-sectional analysis of clinician decision-making was conducted using ARS data from recent oncology summits. When presented with a clinical vignette of a high-risk patient (\(N=97\)), providers demonstrated highly polarized treatment preferences, split nearly evenly between an intensive targeted regimen (46.9%) and immediate chemotherapy (44.8%). A majority of participating clinicians (61.9%) explicitly identified “visceral crisis/liver metastasis” as the primary factor driving this treatment decision. Furthermore, when querying the general rationale for prioritizing 1L chemotherapy over targeted therapies (\(N=102\)), respondents predominantly cited “rapidly progressing disease” (64.7%) and the “presence of symptomatic visceral involvement/crisis” (49.0%).

References

¹Hainmueller J. Entropy Balancing for Causal Effects: A Multivariate Reweighting Method to Produce Balanced Samples in Observational Studies. Polit Anal. 2012;20(1):25-46.

²Chen R, Chen G, Yu M. Entropy Balancing for Causal Generalization with Target Sample Summary Information. Biometrics. 2023;79(1):16-28.

³White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817-838.

⁴VanderWeele TJ, Ding P. Sensitivity analysis in observational research: introducing the E-value. Ann Intern Med. 2017;167(4):268-274.

⁵Imai K, Ratkovic M. Covariate balancing propensity score. J R Stat Soc Series B Stat Methodol. 2014;76(1):243-263.

⁶Hosmer DW Jr, Lemeshow S. A goodness-of-fit test for the multiple logistic regression model. Commun Stat Theory Methods. 1980;9(10):1035-1049.