Breast Cancer: RWE Survival Analysis
Samuel Baird
2026-05-28
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(ebal)
library(cobalt)
library(ggplot2)
library(car)
library(ResourceSelection)
library(pROC)
library(margins)
library(EValue)
library(DT)
library(ggpubr)
library(scales)
library(patchwork)
# Ensure the file path is correct for your local machine
# NOTE: Please update these file paths to match your local environment.
tryCatch({
dt <- as.data.table(read_excel("~/Code_AACR/Data/AKT1-mutatedBreastCancer_Final dataset.xlsx"))
}, error = function(e) {
stop("File not found. Please check your path for 'AKT1-mutatedBreastCancer_Final dataset.xlsx'.")
})
# Concise rename map (script-style)
rename_map <- c(
"Patient ID" = "patient_id", "Primary Race" = "race", "Center" = "center",
"Overall Tumor Grade" = "tumor_grade",
"Has this Patient experienced a Loco-Regional Recurrence (LRR)?" = "lrr",
"Was a biopsy performed of metastatic site?" = "met_biopsy",
"Sequencing Method" = "seq_meth",
"Discordant receptor status with primary" = "discordant",
"Interval between sequencing and metastatic diagnosis (months)" = "seq_int",
"Additional Sequencing Performed" = "add_seq",
"Fraction Genome Altered" = "frac_alt", "HER2 Status" = "HER2",
"AJCC Stage" = "stage", "AKT Mutation Status" = "akt",
"Tumor Mutation Count" = "TMB", "Hormone Receptor Status" = "HR_status",
"ER/PR Receptor Change" = "HR_change", "mTOR Inhibitor Overall" = "mTOR",
"CDK4/6 Inhibitor Overall" = "CDK4/6", "AKT Inhibitor Overall" = "AKT_inhibitor",
"Fulvestrant Overall" = "fulvestrant",
"Aromatase Inhibitor Overall" = "aromatase_inhibitor",
"Age at Initial Diagnosis" = "prim_age",
"Age at Metastatic Diagnosis" = "met_age",
"Therapeutic clinical trial overall" = "clinical_trial",
"Met Site: Bone" = "bone_mets", "Met Site: Lymph Node" = "lymph_mets",
"Met Site: Multiple" = "multiple_mets", "Met Site: Other" = "other_mets",
"Met Site: Soft Tissue" = "soft_mets", "Met Site: Visceral" = "visceral_mets",
"Met Site: Bone Only" = "boneonly_mets", "Met Site: Lung" = "lung_mets",
"Met Site: Brain" = "brain_mets", "Met Site: Liver" = "liver_mets",
"Cancer Type Detailed" = "cancer_type",
"Overall Survival (Months)" = "survival_time",
"Overall Survival Status" = "survival_status"
)
setnames(dt, old = names(rename_map), new = unname(rename_map))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.
dt[, eligibility_flag := HER2 == "Negative" &
(HR_status %in% c("ER+ & PR-", "ER+ & PR+")) &
(stage == "4")]
eligible <- dt[eligibility_flag == TRUE]
# fat = fulvestrant + aromatase inhibitor + targeted therapy
eligible[, fat_flag := (fulvestrant == "Yes" & aromatase_inhibitor == "Yes") &
(mTOR == "Yes" | `CDK4/6` == "Yes" | AKT_inhibitor == "Yes")]
eligible[, therapy_flag := factor(
fifelse(fat_flag == TRUE, "Group A", "Group B"),
levels = c("Group B", "Group A")
)]
eligible[, survival_status := ifelse(survival_status == "1:DECEASED", 1, 0)]
eligible[, .N, therapy_flag]## therapy_flag N
## <fctr> <int>
## 1: Group B 53
## 2: Group A 303. 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.
The streamlined helper below replaces ~25 repetitive blocks. Each
call: (a) recodes “NA” strings to “Unknown” when requested, (b) prints
group distribution, (c) runs Fisher’s exact or chi-square association vs
therapy_flag.
recode_na <- function(var, na_to_unknown = TRUE) {
if (na_to_unknown) eligible[get(var) == "NA", (var) := "Unknown"]
}
assoc_test <- function(var, test = c("fisher", "chisq"),
na_to_unknown = TRUE, simulate = FALSE, seed = NULL,
as_factor = FALSE) {
test <- match.arg(test)
recode_na(var, na_to_unknown)
if (as_factor) eligible[, (var) := as.factor(get(var))]
dis <- eligible[, .N, by = c("therapy_flag", var)][, prop := N / sum(N),
by = therapy_flag]
cat("\n--- ", var, " ---\n", sep = "")
print(dis)
if (!is.null(seed)) set.seed(seed)
res <- switch(
test,
fisher = if (simulate)
fisher.test(eligible[[var]], eligible$therapy_flag,
simulate.p.value = TRUE)
else fisher.test(eligible[[var]], eligible$therapy_flag),
chisq = chisq.test(eligible[[var]], eligible$therapy_flag)
)
print(res)
invisible(res)
}# 1 met_biopsy ~ Fisher, NA->Unknown, factor
assoc_test("met_biopsy", "fisher", as_factor = TRUE)
# 2 HR_status ~ Fisher (no NA)
assoc_test("HR_status", "fisher", na_to_unknown = FALSE)
# 3 bone_mets ~ Fisher, NA->Unknown
assoc_test("bone_mets", "fisher")
# 4 lymph_mets ~ Chisq, NA->Unknown, factor
assoc_test("lymph_mets", "chisq", as_factor = TRUE)
# 5 multiple_mets ~ Chisq, factor (no NA)
assoc_test("multiple_mets", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
# 6 other_mets ~ Chisq, NA->Unknown
assoc_test("other_mets", "chisq")
# 7 soft_mets ~ Fisher, NA->Unknown, factor
assoc_test("soft_mets", "fisher", as_factor = TRUE)
# 8 visceral_mets ~ Chisq, factor
assoc_test("visceral_mets", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
# 9 seq_meth ~ Fisher simulated, factor
assoc_test("seq_meth", "fisher", na_to_unknown = FALSE, simulate = TRUE,
seed = 100, as_factor = TRUE)
# 10 discordant ~ Fisher simulated, factor, NA->Unknown
assoc_test("discordant", "fisher", simulate = TRUE, seed = 101,
as_factor = TRUE)
# 11 TMB: continuous -> categorical, then test
eligible[, TMB := as.numeric(TMB)]
eligible[, TMB := cut(TMB, c(-Inf, 10, Inf), c("low", "high"), right = TRUE)]
eligible[TMB == "NA" | is.na(TMB), TMB := "Unknown"]
assoc_test("TMB", "fisher", na_to_unknown = FALSE, simulate = TRUE, seed = 102)
# 12 frac_alt: continuous -> categorical
eligible[, frac_alt := as.numeric(frac_alt)]
eligible[, frac_alt := cut(frac_alt, c(-Inf, 0.35, Inf),
c("low fra", "high fra"), right = TRUE)]
eligible[frac_alt == "NA" | is.na(frac_alt), frac_alt := "Unknown"]
assoc_test("frac_alt", "fisher", na_to_unknown = FALSE, simulate = TRUE,
seed = 103)
# 13 race
eligible[race %in% c("Black", "Other", "NA"), race := "Non-White"]
assoc_test("race", "chisq", na_to_unknown = FALSE)
# 14 cancer_type
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")]
assoc_test("cancer_type", "chisq", na_to_unknown = FALSE)
# 15-17 liver_mets, brain_mets, lung_mets (NA->Unknown via ifelse)
for (v in c("liver_mets", "brain_mets", "lung_mets")) {
eligible[, (v) := ifelse(get(v) == "NA", "Unknown", get(v))]
}
assoc_test("liver_mets", "fisher", na_to_unknown = FALSE, simulate = TRUE,
seed = 104)
assoc_test("brain_mets", "fisher", na_to_unknown = FALSE)
assoc_test("lung_mets", "chisq", na_to_unknown = FALSE)
# 18 met_age (<=55 vs >55)
eligible[, met_age := cut(met_age, c(-Inf, 55, Inf),
labels = c("<=55", ">55"), right = TRUE)]
assoc_test("met_age", "chisq", na_to_unknown = FALSE)
# 19 tumor_grade collapse
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"]
assoc_test("tumor_grade", "fisher", na_to_unknown = FALSE)
# 20 seq_int (categorical for balance / interpretability)
eligible[, seq_int := cut(seq_int, c(-Inf, 3, Inf),
c("seq int 1-3 months", "seq int 4+ months"),
right = TRUE)]
assoc_test("seq_int", "chisq", na_to_unknown = FALSE)
# 21 add_seq (NOT baseline — exclude from PS model)
assoc_test("add_seq", "fisher", na_to_unknown = FALSE, as_factor = TRUE)
# 22 boneonly_mets
assoc_test("boneonly_mets", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
# 23 akt
assoc_test("akt", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
# 24 HR_change
eligible[HR_change == "NA" | is.na(HR_change) | HR_change == "Unk",
HR_change := "Unknown"]
assoc_test("HR_change", "fisher", na_to_unknown = FALSE, as_factor = TRUE)
# 25 clinical_trial
assoc_test("clinical_trial", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
# ----------------------------------------------------------------------------
# Composite & stabilized variables for the PS model
# ----------------------------------------------------------------------------
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"))]
fisher.test(eligible$met_pattern_stable, eligible$therapy_flag)
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"))]
chisq.test(eligible$tumor_grade_stable, eligible$therapy_flag)
eligible[, therapy_flag_binary := ifelse(therapy_flag == "Group A", 1, 0)]
# 1L Chemotherapy: factor + binary outcome (used in Section 9 diversion analysis)
eligible[, first_chemo := as.factor(`Chemotherapy in first-line treatment`)]
eligible[, first_chemo_binary := ifelse(first_chemo == "Yes", 1, 0)]
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. The unweighted Cox model and proportional hazards assumption test are also presented.
# Standard KM analysis
surv_object <- Surv(time = eligible[, survival_time], event = eligible[, survival_status])
km_fit <- survfit2(surv_object ~ therapy_flag, data = eligible)
# Unweighted Cox model and PH test
cox_model_unweighted <- coxph(surv_object ~ therapy_flag, data = eligible)
summary(cox_model_unweighted) # HR: 0.35 [0.15-0.81], p=0.013## Call:
## coxph(formula = surv_object ~ therapy_flag, data = eligible)
##
## n= 83, number of events= 31
##
## coef exp(coef) se(coef) z Pr(>|z|)
## therapy_flagGroup A -1.041 0.353 0.421 -2.473 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A 0.353 2.833 0.1546 0.8056
##
## Concordance= 0.637 (se = 0.039 )
## Likelihood ratio test= 6.87 on 1 df, p=0.009
## Wald test = 6.12 on 1 df, p=0.01
## Score (logrank) test = 6.64 on 1 df, p=0.01## chisq df p
## therapy_flag 0.733 1 0.39
## GLOBAL 0.733 1 0.39# Kaplan-Meier plot
p1 <- ggsurvplot(
km_fit, data = eligible, pval = TRUE, pval.method = TRUE,
risk.table = TRUE, risk.table.col = "strata", surv.median.line = "hv",
linetype = "solid", linewidth = 0.8, 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_pubr()
)
p1$table <- p1$table + theme(legend.position = "none")
p1
# Reverse KM for Median Follow-Up (with 95% CI extraction)
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.5median_follow_up <- quantile(follow_up_fit, probs = 0.5)
cat("Median follow-up:", round(median_follow_up$quantile, 1),
"months (95% CI:", round(median_follow_up$lower, 1), "-",
round(median_follow_up$upper, 1), ")\n")## Median follow-up: 58.9 months (95% CI: 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). Balance is verified using a Love plot and a balance table, and weight distributions are inspected for extreme values.
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.12eligible[, ws := weightit_model$weights]
# Verify weight sums equal N per group
eligible[, .(Weight_Sum = sum(ws)), by = therapy_flag]## therapy_flag Weight_Sum
## <fctr> <num>
## 1: Group B 53
## 2: Group A 30final_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, binary = "std", stars = "std",
abs = TRUE, var.order = "unadjusted", limits = c(0, 0.6),
line = TRUE, shapes = c(16, 17),
colors = c("#f79c59", "#549ed1"),
sample.names = c("Unweighted", "Weighted"),
var.names = final_labels_weight, drop.distance = FALSE)
# Balance table complementing the Love plot
print(bal.tab(weightit_model, thresholds = 0.1, un = TRUE, abs = TRUE))## Balance Measures
## Type Diff.Un Diff.Adj
## met_age_>55 Binary 0.1239 0
## tumor_grade_stable_Grade 3/Unknown Binary 0.2270 0
## akt_Wildtype Binary 0.1302 0
## met_pattern_stable_Bone-Only Binary 0.0195 0
## met_pattern_stable_Other Sites (Non-Liver) Binary 0.1981 0
## met_pattern_stable_Liver Involved Binary 0.1786 0
## M.Threshold
## met_age_>55 Balanced, <0.1
## tumor_grade_stable_Grade 3/Unknown Balanced, <0.1
## akt_Wildtype Balanced, <0.1
## met_pattern_stable_Bone-Only Balanced, <0.1
## met_pattern_stable_Other Sites (Non-Liver) Balanced, <0.1
## met_pattern_stable_Liver Involved Balanced, <0.1
##
## Balance tally for mean differences
## count
## Balanced, <0.1 6
## Not Balanced, >0.1 0
##
## Variable with the greatest mean difference
## Variable Diff.Adj M.Threshold
## met_pattern_stable_Liver Involved 0 Balanced, <0.1
##
## Effective sample sizes
## Group B Group A
## Unadjusted 53. 30.
## Adjusted 47.33 20.12# Weight distribution histogram (extreme-weight inspection)
ggplot(eligible, aes(x = ws, fill = therapy_flag)) +
geom_histogram(binwidth = 0.05, colour = "black",
alpha = 0.7, position = "identity") +
geom_vline(xintercept = 1, linetype = "dashed", color = "grey30") +
scale_fill_manual(values = c("Group B" = "#549ed1",
"Group A" = "#f79c59"),
name = "Therapy Group") +
labs(title = "Distribution of Entropy Balancing Weights",
subtitle = "Checking for extreme weight outliers and cohort overlap",
x = "Calculated Weight Value (ws)", y = "Number of Patients",
caption = "Dashed line at 1.0 represents the unweighted baseline.") +
theme_classic() +
theme(legend.position = "top",
plot.title = element_text(face = "bold", size = 14),
axis.title = element_text(size = 11))
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.
# Overall weighted KM (for cohort-level median OS and landmark probabilities)
km_fit_weighted_overall <- survfit2(Surv(survival_time, survival_status) ~ 1,
data = eligible, weights = ws)
print(km_fit_weighted_overall)## Call: survfit(formula = Surv(survival_time, survival_status) ~ 1, data = eligible,
## weights = ws)
##
## records n events median 0.95LCL 0.95UCL
## [1,] 83 83 30.8 68.6 52.3 NA## Call: survfit(formula = Surv(survival_time, survival_status) ~ 1, data = eligible,
## weights = ws)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 12 71.7 0.00 1.000 0.0000 1.000 1.000
## 24 60.1 4.27 0.938 0.0278 0.885 0.994
## 36 46.1 7.07 0.819 0.0532 0.721 0.930
## 48 34.9 4.77 0.726 0.0652 0.608 0.865
## 60 21.1 7.23 0.558 0.0792 0.422 0.737
## 72 14.5 2.17 0.496 0.0807 0.361 0.682## strata median lower upper
## 1 All 68.58553 52.33553 NA# Stratified weighted KM
km_fit_weighted <- survfit2(Surv(survival_time, survival_status) ~ therapy_flag,
data = eligible, weights = ws)
print(km_fit_weighted)## Call: survfit(formula = Surv(survival_time, survival_status) ~ therapy_flag,
## data = eligible, weights = ws)
##
## records n events median 0.95LCL 0.95UCL
## therapy_flag=Group B 53 53 21.04 55.9 46.4 NA
## therapy_flag=Group A 30 30 9.78 89.9 52.3 NA## Call: survfit(formula = Surv(survival_time, survival_status) ~ therapy_flag,
## data = eligible, weights = ws)
##
## therapy_flag=Group B
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 53.00 0.00 1.000 0.0000 1.000 1.000
## 12 44.78 0.00 1.000 0.0000 1.000 1.000
## 24 34.55 4.27 0.900 0.0440 0.817 0.990
## 36 23.39 5.45 0.742 0.0778 0.605 0.912
## 48 16.30 4.77 0.579 0.0940 0.421 0.796
## 60 10.95 3.26 0.463 0.0966 0.308 0.697
## 72 6.16 1.60 0.383 0.1004 0.230 0.640
##
## therapy_flag=Group A
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 30.00 0.000 1.000 0.0000 1.000 1.000
## 12 26.96 0.000 1.000 0.0000 1.000 1.000
## 24 25.52 0.000 1.000 0.0000 1.000 1.000
## 36 22.68 1.627 0.933 0.0591 0.824 1.000
## 48 18.57 0.000 0.933 0.0591 0.824 1.000
## 60 10.15 3.977 0.687 0.1395 0.462 1.000
## 72 8.32 0.565 0.649 0.1392 0.426 0.988## strata median lower upper
## 1 therapy_flag=Group B 55.92105 46.38158 NA
## 2 therapy_flag=Group A 89.86842 52.33553 NA# Primary weighted Cox model
cox_model_weighted <- coxph(
Surv(survival_time, survival_status) ~ therapy_flag,
method = "breslow", data = eligible, weights = ws,
robust = TRUE, x = TRUE
)
summary(cox_model_weighted)## Call:
## coxph(formula = Surv(survival_time, survival_status) ~ therapy_flag,
## data = eligible, weights = ws, robust = TRUE, x = TRUE, method = "breslow")
##
## 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).# Proportional hazards assumption: numerical + graphical
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.32ggcoxzph(ph_test, var = "therapy_flag", resid = TRUE, se = TRUE, df = 2,
ggtheme = theme_survminer(),
title = "Proportional Hazards Assumption Test for Group A",
ylab = "Time-varying Log-Hazard Ratio for Group A")
# Build the weighted KM plot for the two-arm comparison
p2 <- ggsurvplot(
km_fit_weighted, data = eligible, pval = TRUE, pval.method = TRUE,
risk.table = TRUE, risk.table.col = "strata",
surv.median.line = "hv", linetype = "solid", size = 0.8,
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_pubr()
)
# Replace weighted risk-table counts with unweighted counts. The survival
# curve still reflects the entropy-balanced pseudo-population; the risk table
# below it shows actual patient counts.
km_fit_unweighted_for_table <- survfit2(
Surv(survival_time, survival_status) ~ therapy_flag, data = eligible)
p_unweighted_full <- ggsurvplot(
km_fit_unweighted_for_table, data = eligible,
risk.table = TRUE, risk.table.col = "strata", break.time.by = 12,
legend.labs = c("Group B (Control)", "Group A (Intensive Regimen)"),
palette = c("#549ed1", "#f79c59"), ggtheme = theme_pubr())
p2$table <- p_unweighted_full$table + theme(legend.position = "none")
p2
## point lower upper
## RR 0.5240857 0.2941386 0.9650718
## E-values 3.2244070 NA 1.22984707. Univariable & Multivariable Cox Regression Models
To explore the clinical features of the cohort, we evaluate the
prognostic impact of our covariates using univariable and multivariable
Cox proportional hazards models across the raw (unweighted) and
entropy-balanced (weighted) samples. A single helper function
(run_uni_cox / run_mv_cox) replaces repetitive
per-covariate fitting blocks.
# --- Define Global Variables and Mappings ---
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"
)7.1 Univariable Cox Models
run_uni_cox <- function(weighted = FALSE) {
fits <- lapply(covariates_filtered, function(cov) {
f <- as.formula(paste("Surv(survival_time, survival_status) ~", cov))
if (weighted) {
summary(coxph(f, data = eligible, weights = ws, robust = TRUE))
} else {
summary(coxph(f, data = eligible))
}
})
names(fits) <- covariates_filtered
tab <- rbindlist(lapply(names(fits), function(cov_name) {
s <- fits[[cov_name]]
if (is.null(s$conf.int) || nrow(s$conf.int) == 0) return(NULL)
res <- as.data.table(s$conf.int, keep.rownames = "Raw_Level")
cm <- as.matrix(s$coefficients)
res[, p_val_numeric := cm[, ncol(cm)]]
res
}))
tab[, `:=`(
Variable = ifelse(Raw_Level %in% names(full_label_map),
full_label_map[Raw_Level], Raw_Level),
HR = round(`exp(coef)`, 2),
Lower = round(`lower .95`, 2),
Upper = round(`upper .95`, 2),
P = ifelse(p_val_numeric < 0.001, "<0.001",
sprintf("%.3f", p_val_numeric))
)]
out <- tab[HR > 0 & Upper < 400, .(Variable, HR, Lower, Upper, P)]
setorder(out, -HR)
out
}
uni_unweighted_final <- run_uni_cox(weighted = FALSE)
datatable(uni_unweighted_final, rownames = FALSE,
options = list(dom = 't', pageLength = 20),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold; color: black;',
'Table 1: Unweighted Univariable OS Analysis'))uni_weighted_final <- run_uni_cox(weighted = TRUE)
datatable(uni_weighted_final, rownames = FALSE,
options = list(dom = 't', pageLength = 20),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold; color: black;',
'Table 2: Weighted Univariable OS Analysis'))7.2 Multivariable Cox Models
mv_formula <- as.formula(
paste("Surv(survival_time, survival_status) ~",
paste(covariates_filtered, collapse = " + "))
)
run_mv_cox <- function(weighted = FALSE) {
if (weighted) {
fit <- coxph(mv_formula, data = eligible, weights = ws, robust = TRUE)
} else {
fit <- coxph(mv_formula, data = eligible)
}
s <- summary(fit)
if (is.null(s$conf.int) || nrow(s$conf.int) == 0) return(NULL)
tab <- as.data.table(s$conf.int, keep.rownames = "Raw_Level")
cm <- as.matrix(s$coefficients)
tab[, p_val_numeric := cm[, ncol(cm)]]
tab[, `:=`(
Variable = ifelse(Raw_Level %in% names(full_label_map),
full_label_map[Raw_Level], Raw_Level),
HR = round(`exp(coef)`, 2),
Lower = round(`lower .95`, 2),
Upper = round(`upper .95`, 2),
P = ifelse(p_val_numeric < 0.001, "<0.001",
sprintf("%.3f", p_val_numeric))
)]
out <- tab[HR > 0 & Upper < 400, .(Variable, HR, Lower, Upper, P)]
setorder(out, -HR)
out
}
mv_unweighted_final <- run_mv_cox(weighted = FALSE)
datatable(mv_unweighted_final, rownames = FALSE,
options = list(dom = 't', pageLength = 20),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold; color: black;',
'Table 3: Unweighted Multivariable OS Analysis'))mv_weighted_final <- run_mv_cox(weighted = TRUE)
datatable(mv_weighted_final, rownames = FALSE,
options = list(dom = 't', pageLength = 20),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold; color: black;',
'Table 4: Weighted Multivariable OS 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 landmarks (3 to 48 months) post-metastatic
diagnosis. Entropy balancing is dynamically re-estimated within
each landmark sub-cohort rather than applying the original
full-cohort weights, because survivors at later landmarks represent a
distinct covariate distribution; static weights would no longer
guarantee balance in the depleted risk set. The landmark balancing model
uses the same four-covariate specification as the primary analysis
(met_age + tumor_grade_stable + akt + met_pattern_stable)
to preserve methodological consistency.
library(htmltools)
# ============================================================================
# DYNAMIC LANDMARK EXPLORATORY SWEEP (3–48 mo, 3-month intervals)
# ============================================================================
landmarks <- seq(3, 48, by = 3)
covars <- c("met_age", "tumor_grade_stable", "met_pattern_stable", "akt")
landmark_formula <- as.formula(paste("therapy_flag ~",
paste(covars, collapse = " + ")))
format_cox_results <- function(model) {
if (is.null(model)) return(list(ci = "NC", p = NA_character_))
s <- summary(model)
hr <- round(s$conf.int[1, 1], 2)
l <- round(s$conf.int[1, 3], 2)
u <- if (is.na(s$conf.int[1, 4]) || is.infinite(s$conf.int[1, 4])) "NE"
else round(s$conf.int[1, 4], 2)
# Prevent p-value rounding to "0" errors
p_raw <- s$coefficients[1, "Pr(>|z|)"]
p_fmt <- if (is.na(p_raw)) NA_character_
else if (p_raw < 0.001) "<0.001"
else as.character(round(p_raw, 3))
list(ci = paste0(hr, " (", l, "-", u, ")"), p = p_fmt)
}
landmark_dynamic_results <- rbindlist(lapply(landmarks, function(t) {
excluded <- eligible[survival_time < t]
n_excl <- nrow(excluded)
n_deaths_pre <- sum(excluded$survival_status == 1)
n_cens_pre <- sum(excluded$survival_status == 0)
d_land <- copy(eligible[survival_time >= t])[, st_adj := survival_time - t]
# keyby ensures consistent ordering and prevents cross-arm vector mismatch
n_per_arm <- d_land[, .N, keyby = therapy_flag]$N
events_per_arm <- d_land[, sum(survival_status), keyby = therapy_flag]$V1
# Relaxed safety thresholds to avoid blanking out late-timeline slices
fittable <- (
nrow(d_land) >= 20 &&
uniqueN(d_land$therapy_flag) == 2 &&
sum(d_land$survival_status) >= 5 &&
length(n_per_arm) == 2 && all(n_per_arm >= 5) &&
length(events_per_arm) == 2 && all(events_per_arm >= 1)
)
if (!fittable) {
return(data.table(
Time = paste0(t, " Mo"),
N = nrow(d_land),
Events = sum(d_land$survival_status),
`N Excluded` = n_excl,
`Deaths Pre-Landmark` = n_deaths_pre,
`Censored Pre-Landmark` = n_cens_pre,
`Unweighted HR (95% CI)` = "Underpowered", `Unweighted P-val` = NA_character_,
`Reweighted HR (95% CI)` = "Underpowered", `Reweighted P-val` = NA_character_
))
}
fit_unw <- tryCatch(
coxph(Surv(st_adj, survival_status) ~ therapy_flag, data = d_land),
error = function(e) NULL)
unw <- format_cox_results(fit_unw)
# Warnings retained (not suppressed) to catch structural balancing failures
eb_fit_lm <- tryCatch(
weightit(landmark_formula, data = d_land, method = "ebal", estimand = "ATE"),
error = function(e) NULL)
fit_w <- if (!is.null(eb_fit_lm)) {
d_land[, ws_dynamic := eb_fit_lm$weights]
tryCatch(coxph(Surv(st_adj, survival_status) ~ therapy_flag,
data = d_land, weights = ws_dynamic, robust = TRUE),
error = function(e) NULL)
} else NULL
wtd <- format_cox_results(fit_w)
data.table(
Time = paste0(t, " Mo"),
N = nrow(d_land),
Events = sum(d_land$survival_status),
`N Excluded` = n_excl,
`Deaths Pre-Landmark` = n_deaths_pre,
`Censored Pre-Landmark` = n_cens_pre,
`Unweighted HR (95% CI)` = unw$ci, `Unweighted P-val` = unw$p,
`Reweighted HR (95% CI)` = wtd$ci, `Reweighted P-val` = wtd$p
)
}))
# Dynamic Sweep Interactive DT Table
datatable(
landmark_dynamic_results,
rownames = FALSE,
options = list(dom = 't', pageLength = 20),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold; color: black;',
'Table 5: Final Side-by-Side Dynamic Landmark Analysis Summary'
)
)cat("Note: 'NE' = Not Estimable due to upper bound infinity.\n",
" 'NC' = Cox model failed (non-convergence or other error).\n",
" 'Underpowered' = below structural safety thresholds for matrix estimation.\n",
" 'Deaths Pre-Landmark' = deaths occurring before the landmark time (excluded).\n",
" 'Censored Pre-Landmark' = administrative censorings before landmark (excluded).\n")## Note: 'NE' = Not Estimable due to upper bound infinity.
## 'NC' = Cox model failed (non-convergence or other error).
## 'Underpowered' = below structural safety thresholds for matrix estimation.
## 'Deaths Pre-Landmark' = deaths occurring before the landmark time (excluded).
## 'Censored Pre-Landmark' = administrative censorings before landmark (excluded).8.1 Focal 9-Month Landmark
Rationale: in Stage IV HR+/HER2- mBC, first-line systemic therapy is initiated within weeks of metastatic diagnosis. Median 1L duration in this cohort was 7.99 months (Group A) and 4.26 months (Group B). By 9 months, the majority of patients in both arms have completed 1L and transitioned to 2L, bracketing the realistic upper bound of the differential classification window during which Group A status was still being established through cumulative agent exposure. A 9-month landmark therefore sits downstream of the full classification window but upstream of substantial treatment-effect accrual on survival. A 12-month landmark would extend into the period during which treatment effects accrue, inflating Type II error without addressing additional immortal time bias.
landmark_time <- 9
eligible_landmarked <- copy(eligible[survival_time >= landmark_time])
if (nrow(eligible_landmarked) == 0)
stop("No patients remain after ", landmark_time, "-month landmark.")
message(nrow(eligible_landmarked), " patients in landmark cohort.")
eligible_landmarked[, survival_time_adj := survival_time - landmark_time]
max_adjusted_time <- max(eligible_landmarked$survival_time_adj, na.rm = TRUE)
plot_breaks_adj <- seq(0, max_adjusted_time, by = 12)
plot_labels_orig <- seq(landmark_time, landmark_time + max_adjusted_time,
by = 12)
eligible_landmarked[, follow_up_status := ifelse(survival_status == 1, 0, 1)]
follow_up_fit_landmarked <- survfit(
Surv(survival_time_adj, follow_up_status) ~ 1,
data = eligible_landmarked)
median_fu_lm <- quantile(follow_up_fit_landmarked, probs = 0.5)
# Interactive DT for Follow-up Stats
fu_summary <- data.table(
Metric = c("Median Follow-up", "95% Lower Bound", "95% Upper Bound"),
Months = c(round(median_fu_lm$quantile, 1),
round(median_fu_lm$lower, 1),
round(median_fu_lm$upper, 1))
)
datatable(
fu_summary,
rownames = FALSE,
options = list(dom = 't'),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold; color: black;',
paste('Table 6: Follow-up Time Summary from', landmark_time, '-Month Landmark')
)
)# Helper: build a landmark ggsurvplot in either unweighted/weighted mode
build_landmark_plot <- function(km_fit, title_txt) {
pp <- ggsurvplot(
km_fit, data = eligible_landmarked, pval = TRUE, pval.method = TRUE,
risk.table = TRUE, risk.table.col = "strata", surv.median.line = "hv",
break.time.by = 12, palette = c("#549ed1", "#f79c59"),
title = title_txt,
legend.labs = c("Group B (Control)", "Group A (Intensive Regimen)"),
ggtheme = theme_pubr())
pp$plot <- pp$plot +
scale_x_continuous(breaks = plot_breaks_adj, labels = plot_labels_orig) +
labs(x = "Time from Initial Diagnosis (Months)",
caption = paste0("Note: Restricted to survivors at ",
landmark_time, " months."))
pp$table <- pp$table +
scale_x_continuous(breaks = plot_breaks_adj, labels = plot_labels_orig) +
labs(x = "Time from Initial Diagnosis (Months)")
pp
}
# Unweighted 9-month landmark (clock reset to 0 at the landmark)
cat("\n--- Unweighted 9-Month Landmarked OS Analysis ---\n")##
## --- Unweighted 9-Month Landmarked OS Analysis ---km_fit_unweighted_landmarked <- survfit2(
Surv(survival_time_adj, survival_status) ~ therapy_flag,
data = eligible_landmarked)
print(km_fit_unweighted_landmarked)## Call: survfit(formula = Surv(survival_time_adj, survival_status) ~
## therapy_flag, data = eligible_landmarked)
##
## n events median 0.95LCL 0.95UCL
## therapy_flag=Group B 45 21 46.9 37.4 NA
## therapy_flag=Group A 30 10 112.6 68.7 NA## Call: survfit(formula = Surv(survival_time_adj, survival_status) ~
## therapy_flag, data = eligible_landmarked)
##
## therapy_flag=Group B
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 45 0 1.000 0.0000 1.000 1.000
## 12 38 3 0.931 0.0386 0.858 1.000
## 24 26 6 0.769 0.0680 0.647 0.915
## 36 18 3 0.670 0.0799 0.530 0.847
## 48 12 5 0.484 0.0913 0.334 0.700
## 60 7 2 0.385 0.0963 0.236 0.628
## 72 5 0 0.385 0.0963 0.236 0.628
##
## therapy_flag=Group A
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 30 0 1.000 0.0000 1.000 1.000
## 12 28 0 1.000 0.0000 1.000 1.000
## 24 24 1 0.960 0.0392 0.886 1.000
## 36 22 1 0.920 0.0543 0.820 1.000
## 48 14 3 0.760 0.0952 0.595 0.972
## 60 12 1 0.702 0.1043 0.524 0.939
## 72 8 1 0.624 0.1183 0.430 0.905##
## Unweighted 9-month landmark median OS (from landmark, months):## strata median lower upper
## 1 therapy_flag=Group B 46.92105 37.38158 NA
## 2 therapy_flag=Group A 112.57895 68.66447 NAcox_model_unweighted_landmarked <- coxph(
Surv(survival_time_adj, survival_status) ~ therapy_flag,
data = eligible_landmarked)
print(summary(cox_model_unweighted_landmarked))## Call:
## coxph(formula = Surv(survival_time_adj, survival_status) ~ therapy_flag,
## data = eligible_landmarked)
##
## n= 75, number of events= 31
##
## coef exp(coef) se(coef) z Pr(>|z|)
## therapy_flagGroup A -1.041 0.353 0.421 -2.473 0.0134 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A 0.353 2.833 0.1546 0.8056
##
## Concordance= 0.637 (se = 0.039 )
## Likelihood ratio test= 6.87 on 1 df, p=0.009
## Wald test = 6.12 on 1 df, p=0.01
## Score (logrank) test = 6.64 on 1 df, p=0.01p_landmarked_km_unweighted <- build_landmark_plot(
km_fit_unweighted_landmarked,
paste0("Unweighted Overall Survival from ", landmark_time, " Months"))
p_landmarked_km_unweighted$table <- p_landmarked_km_unweighted$table +
theme(legend.position = "none")
print(p_landmarked_km_unweighted)
# Dynamic weighted 9-month landmark
cat("\n--- Weighted 9-Month Landmarked OS Analysis (Entropy Balanced) ---\n")##
## --- Weighted 9-Month Landmarked OS Analysis (Entropy Balanced) ---formula_to_use <- if (exists("eb_fit")) eb_fit$formula else landmark_formula
eb_fit_landmarked <- weightit(formula_to_use, data = eligible_landmarked,
method = "ebal", estimand = "ATE")
eligible_landmarked[, ws_landmarked := eb_fit_landmarked$weights]
print(summary(eb_fit_landmarked))## Summary of weights
##
## - Weight ranges:
##
## Min Max
## treated 0.217 |---------------------------| 2.495
## control 0.352 |---------------------| 2.094
##
## - Units with the 5 most extreme weights by group:
##
## 19 22 14 12 3
## treated 1.519 1.519 1.905 2.441 2.495
## 22 33 13 4 9
## control 1.527 1.527 1.752 2.094 2.094
##
## - Weight statistics:
##
## Coef of Var MAD Entropy # Zeros
## treated 0.580 0.452 0.152 0
## control 0.398 0.308 0.072 0
##
## - Effective Sample Sizes:
##
## Control Treated
## Unweighted 45. 30.
## Weighted 38.97 22.63km_fit_weighted_landmarked <- survfit2(
Surv(survival_time_adj, survival_status) ~ therapy_flag,
data = eligible_landmarked, weights = ws_landmarked)
print(km_fit_weighted_landmarked)## Call: survfit(formula = Surv(survival_time_adj, survival_status) ~
## therapy_flag, data = eligible_landmarked, weights = ws_landmarked)
##
## records n events median 0.95LCL 0.95UCL
## therapy_flag=Group B 45 45 20.20 46.9 37.4 NA
## therapy_flag=Group A 30 30 9.51 112.6 58.4 NA## Call: survfit(formula = Surv(survival_time_adj, survival_status) ~
## therapy_flag, data = eligible_landmarked, weights = ws_landmarked)
##
## therapy_flag=Group B
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 45.00 0.00 1.000 0.0000 1.000 1.000
## 12 37.90 2.54 0.942 0.0338 0.878 1.000
## 24 26.45 5.10 0.800 0.0651 0.683 0.939
## 36 17.95 4.14 0.661 0.0895 0.507 0.862
## 48 11.34 5.35 0.464 0.0981 0.307 0.702
## 60 6.83 1.45 0.389 0.1018 0.233 0.650
## 72 4.59 0.00 0.389 0.1018 0.233 0.650
##
## therapy_flag=Group A
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 30.00 0.000 1.000 0.0000 1.000 1.000
## 12 25.99 0.000 1.000 0.0000 1.000 1.000
## 24 23.13 1.346 0.945 0.0534 0.846 1.000
## 36 21.95 0.217 0.936 0.0540 0.836 1.000
## 48 11.75 3.590 0.721 0.1231 0.516 1.000
## 60 10.36 0.599 0.682 0.1247 0.476 0.976
## 72 6.86 0.676 0.620 0.1312 0.410 0.939##
## Weighted 9-month landmark median OS (from landmark, months):## strata median lower upper
## 1 therapy_flag=Group B 46.92105 37.38158 NA
## 2 therapy_flag=Group A 112.57895 58.43421 NAcox_model_weighted_landmarked <- coxph(
Surv(survival_time_adj, survival_status) ~ therapy_flag,
data = eligible_landmarked, weights = ws_landmarked, robust = TRUE)
print(summary(cox_model_weighted_landmarked))## Call:
## coxph(formula = Surv(survival_time_adj, survival_status) ~ therapy_flag,
## data = eligible_landmarked, weights = ws_landmarked, robust = TRUE)
##
## n= 75, number of events= 31
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## therapy_flagGroup A -1.0257 0.3585 0.4376 0.4515 -2.272 0.0231 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A 0.3585 2.789 0.148 0.8687
##
## Concordance= 0.636 (se = 0.044 )
## Likelihood ratio test= 6.21 on 1 df, p=0.01
## Wald test = 5.16 on 1 df, p=0.02
## Score (logrank) test = 5.95 on 1 df, p=0.01, Robust = 5.71 p=0.02
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).## chisq df p
## therapy_flag 0.955 1 0.33
## GLOBAL 0.955 1 0.33p_landmarked_km_weighted <- build_landmark_plot(
km_fit_weighted_landmarked,
paste0("Adjusted Overall Survival (Reweighted Entropy Balanced) from ",
landmark_time, " Months"))
# Replace weighted landmark risk-table counts with unweighted counts.
p_unweighted_landmark_full <- ggsurvplot(
km_fit_unweighted_landmarked, data = eligible_landmarked,
risk.table = TRUE, risk.table.col = "strata", break.time.by = 12,
legend.labs = c("Group B (Control)", "Group A (Intensive Regimen)"),
palette = c("#549ed1", "#f79c59"), ggtheme = theme_pubr())
p_landmarked_km_weighted$table <- p_unweighted_landmark_full$table +
scale_x_continuous(breaks = plot_breaks_adj, labels = plot_labels_orig) +
labs(x = "Time from Initial Diagnosis (Months)") +
theme(legend.position = "none")
print(p_landmarked_km_weighted)
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. Goodness-of-fit is assessed via Hosmer-Lemeshow and AUC, with multicollinearity assessed via VIF.
All AME forest plots in this section use the same point size, errorbar geometry, x-axis breaks, gridline treatment, and caption styling. Limits may be set per-plot or computed across paired panels for shared-axis comparisons.
# ============================================================================
# Shared AME plot styling
# ============================================================================
ame_point_size <- 3.5
ame_errorbar_width <- 0.15
ame_errorbar_lw <- 0.7
ame_vline_color <- "gray60"
ame_break_step <- 0.2
# Compute symmetric, rounded x-axis limits + breaks from an AME data.table.
ame_xaxis <- function(dt, pad = 0.05, step = ame_break_step) {
lo <- floor((min(dt$lower) - pad) * 10) / 10
hi <- ceiling((max(dt$upper) + pad) * 10) / 10
list(limits = c(lo, hi), breaks = seq(lo, hi, by = step))
}
# Build a styled AME forest plot.
build_ame_plot <- function(dt, ax, title_txt, x_label,
subtitle_txt = NULL, caption_txt = NULL) {
ggplot(dt, aes(x = AME, y = reorder(factor, AME))) +
geom_vline(xintercept = 0, linetype = "dashed", color = ame_vline_color) +
geom_errorbar(aes(xmin = lower, xmax = upper),
width = ame_errorbar_width, color = "black",
linewidth = ame_errorbar_lw) +
geom_point(size = ame_point_size, color = "black") +
scale_x_continuous(labels = scales::percent_format(accuracy = 1),
limits = ax$limits, breaks = ax$breaks) +
labs(title = title_txt, subtitle = subtitle_txt,
x = x_label, y = NULL, caption = caption_txt) +
theme_pubr() +
theme(plot.title = element_text(face = "bold", size = 12),
plot.subtitle = element_text(size = 10),
plot.caption = element_text(hjust = 0, size = 9,
face = "italic", color = "grey30"),
axis.text.y = element_text(size = 10))
}
# Helper that fits the GLM, computes diagnostics, and returns the styled AME
# plot + the AME data.table (for downstream shared-axis combined figures).
run_ame_panel <- function(target_group, level_order, label_map,
title_txt, x_lab, ax = NULL, roc_title,
subtitle_txt = NULL, caption_txt = NULL) {
# 1. Outcome
eligible[, therapy_flag_binary := ifelse(therapy_flag == target_group, 1, 0)]
# 2. Apply this panel's reference levels
eligible[, met_age := factor(met_age, levels = level_order$met_age)]
eligible[, tumor_grade_stable := factor(tumor_grade_stable, levels = level_order$tumor_grade_stable)]
eligible[, akt := factor(akt, levels = level_order$akt)]
eligible[, met_pattern_stable := factor(met_pattern_stable, levels = level_order$met_pattern_stable)]
# 3. GLM
fit <- glm(therapy_flag_binary ~ met_age + tumor_grade_stable +
akt + met_pattern_stable,
data = eligible, family = binomial())
print(summary(fit))
# 4. Diagnostics
print(as.data.table(vif(fit)))
hl <- hoslem.test(eligible$therapy_flag_binary, fitted(fit), g = 6)
print(hl)
roc_obj <- roc(response = eligible$therapy_flag_binary,
predictor = fitted(fit),
plot = TRUE, print.auc = TRUE,
col = "black", main = roc_title)
auc_val <- round(auc(roc_obj), 3)
if (is.null(subtitle_txt)) {
subtitle_txt <- paste0("AUC: ", auc_val,
" | HL p-val: ", round(hl$p.value, 3))
}
# 5. AME table
ame_tbl <- as.data.table(summary(margins(fit)))
ame_tbl[, factor := ifelse(factor %in% names(label_map),
label_map[factor], factor)]
# 6. Plot (auto-compute axis if not passed in)
if (is.null(ax)) ax <- ame_xaxis(ame_tbl)
pl <- build_ame_plot(
dt = ame_tbl,
ax = ax,
title_txt = title_txt,
x_label = x_lab,
subtitle_txt = subtitle_txt,
caption_txt = caption_txt)
print(pl)
list(plot = pl, ame = ame_tbl, auc = auc_val, hl = hl)
}Group B (Control) Selection Drivers
panel_b <- run_ame_panel(
target_group = "Group B",
level_order = list(
met_age = c("<=55", ">55"),
tumor_grade_stable = c("Grade 1-2", "Grade 3/Unknown"),
akt = c("Mutant", "Wildtype"),
met_pattern_stable = c("Bone-Only", "Liver Involved",
"Other Sites (Non-Liver)")
),
label_map = 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"
),
title_txt = "Clinical Drivers of Group B (Control) Selection",
x_lab = "Change in Probability of Receiving Group B",
roc_title = "Group B Assignment Prediction (ROC)",
caption_txt = 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.")
)##
## Call:
## glm(formula = therapy_flag_binary ~ met_age + tumor_grade_stable +
## akt + met_pattern_stable, family = binomial(), data = eligible)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1480 0.6792 -0.218 0.8275
## met_age>55 0.9524 0.5390 1.767 0.0772
## tumor_grade_stableGrade 3/Unknown -1.1609 0.5219 -2.224 0.0261
## aktWildtype 1.0201 0.6244 1.634 0.1023
## met_pattern_stableLiver Involved 1.6052 0.8724 1.840 0.0658
## met_pattern_stableOther Sites (Non-Liver) -0.4592 0.5467 -0.840 0.4009
##
## (Intercept)
## met_age>55 .
## tumor_grade_stableGrade 3/Unknown *
## aktWildtype
## met_pattern_stableLiver Involved .
## met_pattern_stableOther Sites (Non-Liver)
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 108.605 on 82 degrees of freedom
## Residual deviance: 93.133 on 77 degrees of freedom
## AIC: 105.13
##
## Number of Fisher Scoring iterations: 4
##
## GVIF Df GVIF^(1/(2*Df))
## <num> <num> <num>
## 1: 1.108242 1 1.052731
## 2: 1.069341 1 1.034089
## 3: 1.084605 1 1.041444
## 4: 1.095890 2 1.023156
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: eligible$therapy_flag_binary, fitted(fit)
## X-squared = 3.1412, df = 4, p-value = 0.5345
datatable(panel_b$ame[, .(Variable = factor,
"Marginal Effect" = round(AME, 3),
Lower = round(lower, 3), Upper = round(upper, 3),
"P-Value" = round(p, 3))],
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'))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.
Note: factor levels are reset here because the Group B panel applied its own reference levels. Each AME panel needs its own contrast reference; releveling here is intentional, not duplicative.
# Apply Standard-Risk Reference Levels for the 1L chemo panel
eligible[, met_age := factor(met_age, levels = c("<=55", ">55"))]
eligible[, tumor_grade_stable := factor(tumor_grade_stable, levels = c("Grade 1-2", "Grade 3/Unknown"))]
eligible[, akt := factor(akt, levels = c("Mutant", "Wildtype"))]
eligible[, met_pattern_stable := factor(met_pattern_stable,
levels = c("Bone-Only",
"Liver Involved",
"Other Sites (Non-Liver)"))]
# Fit the 1L Chemotherapy Selection Model
glm_chemo_1L <- glm(
first_chemo_binary ~ met_age + tumor_grade_stable + akt + met_pattern_stable,
data = eligible,
family = binomial()
)
cat("\n--- 1L Chemotherapy Selection Model Summary ---\n")##
## --- 1L Chemotherapy Selection Model Summary ---##
## Call:
## glm(formula = first_chemo_binary ~ met_age + tumor_grade_stable +
## akt + met_pattern_stable, family = binomial(), data = eligible)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.7876 0.8229 -2.172 0.02984
## met_age>55 -0.5188 0.5671 -0.915 0.36028
## tumor_grade_stableGrade 3/Unknown 0.1872 0.5516 0.339 0.73426
## aktWildtype -0.6495 0.6330 -1.026 0.30487
## met_pattern_stableLiver Involved 2.5301 0.8034 3.149 0.00164
## met_pattern_stableOther Sites (Non-Liver) 1.9682 0.7257 2.712 0.00668
##
## (Intercept) *
## met_age>55
## tumor_grade_stableGrade 3/Unknown
## aktWildtype
## met_pattern_stableLiver Involved **
## met_pattern_stableOther Sites (Non-Liver) **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 97.973 on 82 degrees of freedom
## Residual deviance: 81.180 on 77 degrees of freedom
## AIC: 93.18
##
## Number of Fisher Scoring iterations: 5# Multicollinearity Check (VIF)
vif_check_chemo <- as.data.table(vif(glm_chemo_1L))
cat("\n--- Multicollinearity Check (1L Chemo VIF) ---\n")##
## --- Multicollinearity Check (1L Chemo VIF) ---## GVIF Df GVIF^(1/(2*Df))
## <num> <num> <num>
## 1: 1.062994 1 1.031016
## 2: 1.020218 1 1.010059
## 3: 1.044920 1 1.022213
## 4: 1.061436 2 1.015017# Goodness of Fit & ROC Diagnostics
hoslem_res_chemo <- hoslem.test(eligible$first_chemo_binary, fitted(glm_chemo_1L), g = 6)
cat("\n--- Hosmer-Lemeshow Goodness of Fit ---\n")##
## --- Hosmer-Lemeshow Goodness of Fit ---##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: eligible$first_chemo_binary, fitted(glm_chemo_1L)
## X-squared = 1.2115, df = 4, p-value = 0.8762roc_obj_chemo <- roc(
response = eligible$first_chemo_binary,
predictor = fitted(glm_chemo_1L),
plot = TRUE,
print.auc = TRUE,
col = "black",
main = "1L Chemotherapy Prediction (ROC)"
)
auc_chemo_1L <- round(auc(roc_obj_chemo), 3)
# AME Generation & Standardized Mapping
ame_table_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_table_chemo[, factor := ifelse(factor %in% names(labels_chemo),
labels_chemo[factor], factor)]
datatable(ame_table_chemo[, .(Variable = factor,
"Marginal Effect" = round(AME, 3),
Lower = round(lower, 3), Upper = round(upper, 3),
"P-Value" = round(p, 3))],
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'))# Standalone 1L Chemo plot (per-plot axis)
ax_1L_chemo <- ame_xaxis(ame_table_chemo)
chemo_subtitle <- paste0("AUC: ", auc_chemo_1L,
" | HL p-val: ", round(hoslem_res_chemo$p.value, 3))
chemo_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 reference, holding other factors constant.\n",
"Baseline references reflect: Bone-Only, AKT Mutant, Age <= 55, and Grade 1-2."
)
chemo_plot <- build_ame_plot(
dt = ame_table_chemo,
ax = ax_1L_chemo,
title_txt = "Clinical Drivers of 1L Chemotherapy Selection",
x_label = "Change in Probability of Receiving 1L Chemotherapy",
subtitle_txt = chemo_subtitle,
caption_txt = chemo_caption
)
print(chemo_plot)
Combined Multi-Panel Figures
# --- Survival panel: unweighted + weighted + 9-month landmark unweighted/weighted
extract_and_combine_ggsurvplot_with_table <-
function(ggsurv_obj, label_char = NULL) {
plot_km <- ggsurv_obj$plot
plot_rt <- ggsurv_obj$table
if (!is.null(label_char)) {
current_title <- plot_km$labels$title
plot_km <- plot_km + labs(title = paste0("(", label_char, ") ",
current_title))
}
plot_rt <- plot_rt +
theme(plot.title = element_blank(), axis.title.x = element_blank())
(plot_km / plot_rt) + plot_layout(ncol = 1, heights = c(3, 1))
}
group_a_panel <- extract_and_combine_ggsurvplot_with_table(p1, "a")
group_b_panel <- extract_and_combine_ggsurvplot_with_table(p2, "b")
group_c_panel <- extract_and_combine_ggsurvplot_with_table(
p_landmarked_km_unweighted, "c")
group_d_panel <- extract_and_combine_ggsurvplot_with_table(
p_landmarked_km_weighted, "d")
combined_survival_figure_with_tables <- wrap_plots(
list(group_a_panel, group_b_panel, group_c_panel, group_d_panel),
ncol = 2)
print(combined_survival_figure_with_tables)
ggsave("combined_survival_plots_with_tables.png",
combined_survival_figure_with_tables,
width = 19, height = 15, dpi = 300)
# --- AME combined figure: Group B + 1L Chemo on a shared x-axis ---
# Pool the lower and upper bounds from both AME tables to compute a single
# master axis. Both plots are re-rendered on this shared axis so the
# combined figure has truly comparable scales.
pooled_bounds <- rbind(
panel_b$ame[, .(lower, upper)],
ame_table_chemo[, .(lower, upper)]
)
ax_master <- ame_xaxis(pooled_bounds)
chemo_plot_shared <- build_ame_plot(
dt = ame_table_chemo,
ax = ax_master,
title_txt = "Clinical Drivers of 1L Chemotherapy Selection",
x_label = "Change in Probability of Receiving 1L Chemotherapy"
)
plot_b_shared <- build_ame_plot(
dt = panel_b$ame,
ax = ax_master,
title_txt = "Clinical Drivers of Group B (Control) Selection",
x_label = "Change in Probability of Receiving Group B"
)
x_axis_theme <- theme(axis.title.x = element_text(size = 10, face = "bold"))
combined_ame_figure <- (plot_b_shared + x_axis_theme) +
(chemo_plot_shared + x_axis_theme) +
plot_layout(ncol = 2, guides = "collect") +
plot_annotation(tag_levels = "a")
print(combined_ame_figure)
10. 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.
# NOTE: Please update this file path to match your local environment.
file_path <- "~/Code_AACR/Data/2026_FebIIandMarchOncSummit_DataPullRequest_Combined_vFinal.xlsx"
# February Summit (PME)
dt_pme <- as.data.table(read_xlsx(file_path, sheet = 1, na = " "))[-1, ]
# March Summit (PME)
dt2_pme <- as.data.table(read_xlsx(file_path, sheet = 5, na = " "))[-1, ]
# Combined Summits
dt3_pme <- rbindlist(list(dt_pme, dt2_pme))
# Define Community versus Academic Providers
dt3_pme[, groups := fcase(
grepl("community", P4, ignore.case = TRUE), "community",
default = "academic"
)]
# Load and merge response data
dt3_res <- rbindlist(list(
as.data.table(read_xlsx(file_path, sheet = 3, na = " ")),
as.data.table(read_xlsx(file_path, sheet = 7, na = " "))
))
dt3_pme[, ROID := as.character(ROID)]
dt3_res[, ROID := as.character(ROID)]
combo <- merge(dt3_pme, dt3_res, by = "ROID", all = FALSE)
# --- Dorinda question: standardize choice labels ---
long_q_name <- names(combo)[grep("Dorinda", names(combo))][1]
setnames(combo, old = long_q_name, new = "Dorinda_Choice")
combo[, Dorinda_Choice := fcase(
grepl("Chemotherapy", Dorinda_Choice, ignore.case = TRUE), "Chemotherapy-based regimen",
grepl("Single-agent|Singleagent", Dorinda_Choice, ignore.case = TRUE), "Single-agent endocrine therapy",
default = Dorinda_Choice
)]
# Filter for active breast cancer treaters
dt_plot <- combo[Dorinda_Choice != "I do not manage patients with breast cancer" &
!is.na(Dorinda_Choice)]
# --- Plot 1: Preferred 1L therapy for Dorinda ---
dorinda_total <- dt_plot[, .(N = .N), by = .(choice = Dorinda_Choice)][,
percentage := round((N / sum(N) * 100), 1)]
p1_ars <- ggplot(dorinda_total, aes(x = reorder(choice, -N), y = percentage)) +
geom_col(fill = "#2C5282", color = "black") +
geom_text(aes(label = paste0(N, " (", percentage, "%)")),
vjust = -0.7, size = 3, fontface = "bold") +
scale_y_continuous(expand = expansion(mult = c(0, 0.2)),
labels = function(x) paste0(x, "%")) +
scale_x_discrete(labels = label_wrap(15)) +
theme_classic() +
labs(title = "Dorinda Patient Case: Preferred 1L Therapy (Combined)",
x = "Choice", y = "Percentage (%)",
caption = paste0("Total N = ", sum(dorinda_total$N)))
print(p1_ars)
# # --- Plot 2: Characteristics influencing 1L selection for Dorinda ---
char_q <- "Q5. Which of the following baseline clinical characteristics most influenced your 1L treatment selection for Dorinda? Please select up to 2. (Multiple Choice - Multiple Response)"
dt_filtered <- combo[!is.na(get(char_q)) &
get(char_q) != "" &
!(get(char_q) %like% "I do not manage patients with breast cancer")]
char_long <- dt_filtered[, .(
Category = trimws(unlist(strsplit(as.character(get(char_q)), "\r?\n")))
), by = ROID]
char_long[, Category := fcase(
Category %like% "Age \\(i.e.", "Age (i.e., 55 years old)",
Category %like% "Menopausal status", "Menopausal status (i.e., postmenopausal)",
default = Category
)]
total_respondents <- uniqueN(dt_filtered$ROID)
char_combined <- char_long[, .(N = .N), by = Category]
char_combined[, percentage := (N / total_respondents) * 100]
p4_ars <- ggplot(char_combined, aes(x = reorder(Category, -N), y = percentage)) +
geom_col(fill = "#2C5282", color = "black", width = 0.7) +
geom_text(aes(label = paste0(N, " (", round(percentage, 1), "%)")),
vjust = -0.7, size = 3.5, fontface = "bold") +
theme_classic() +
scale_y_continuous(expand = expansion(mult = c(0, 0.2)),
labels = function(x) paste0(x, "%")) +
scale_x_discrete(labels = label_wrap(15)) +
theme(axis.text.x = element_text(vjust = 1, hjust = 0.5),
plot.title = element_text(face = "bold")) +
labs(title = "Characteristics Influencing 1L Treatment Selection for Dorinda",
subtitle = "Combined Results (Multiple Response: Select up to 2)",
x = "Clinical Characteristic", y = "Percentage (%)",
caption = paste0("Total N (Unique Respondents) = ", total_respondents))
print(p4_ars)
# # --- Plot 3: Reasons for 1L chemotherapy over targeted therapy ---
reason_q <- "Q58. What are the most common reasons why your patients with de novo HR+/HER2– metastatic breast cancer who are eligible for targeted therapy (i.e., endocrine therapy + a CDK4/6 inhibitor) receive a chemotherapy-based regimen in the 1L setting instead? Please select up to 3. (Multiple Choice - Multiple Response)"
dt_filtered_58 <- combo[!is.na(get(reason_q)) &
get(reason_q) != "" &
!(get(reason_q) %like% "I do not manage patients with breast cancer")]
reasons_long <- dt_filtered_58[, .(
Category = trimws(unlist(strsplit(as.character(get(reason_q)), "\r?\n")))
), by = ROID]
reasons_long[, Category := fcase(
Category %like% "Visceral crisis|visceral crisis", "Presence of visceral crisis",
Category %like% "Patient preference|Patient request", "Patient preference",
Category %like% "Symptomatic|symptom burden", "High symptom burden",
Category %like% "Cost|Affordability|insurance", "Financial/Insurance barriers",
Category %like% "Presence of endorgan|Presence of end-organ", "Presence of end-organ dysfunction",
default = Category
)]
total_respondents_58 <- uniqueN(dt_filtered_58$ROID)
reasons_combined <- reasons_long[, .(N = .N), by = Category]
reasons_combined[, percentage := (N / total_respondents_58) * 100]
p5_ars <- ggplot(reasons_combined, aes(x = reorder(Category, -N), y = percentage)) +
geom_col(fill = "#2C5282", color = "black", width = 0.7) +
geom_text(aes(label = paste0(N, " (", round(percentage, 1), "%)")),
vjust = -0.7, size = 3.5, fontface = "bold") +
theme_classic() +
scale_y_continuous(expand = expansion(mult = c(0, 0.2)),
labels = function(x) paste0(x, "%")) +
scale_x_discrete(labels = label_wrap(15)) +
theme(axis.text.x = element_text(vjust = 1, hjust = 0.5),
plot.title = element_text(face = "bold", size = 12),
plot.subtitle = element_text(size = 10)) +
labs(title = "Reasons for Utilizing 1L Chemotherapy over Targeted Therapy",
subtitle = "HR+/HER2– De Novo Patients (Combined Results: Select up to 3)",
x = "Reason for Selecting Chemotherapy", y = "Percentage (%)",
caption = paste0("Total N (Unique Respondents) = ", total_respondents_58))
print(p5_ars)
11. Manuscript Drafts
Statistical Methods
Overall survival (OS) was measured from the date of metastatic diagnosis to death from any cause or last follow-up. Survival distributions were estimated using the Kaplan-Meier method. As the primary analysis, hazard ratios (HR) comparing Group A to Group B were estimated via Cox proportional hazards regression, beginning with an unadjusted univariable model. To address selection bias, Entropy Balancing¹ was utilized to estimate patient-level weights targeting the Average Treatment Effect (ATE). This approach adjusts for sample moments without iterative, misspecification-prone propensity-score modeling², rendering it advantageous for maximizing statistical efficiency and preserving effective sample size. Four baseline covariates were specified a priori: age at metastatic diagnosis (\(\le 55\) vs. \(> 55\)), initial tumor grade (Grade 1–2 vs. Grade 3/Unknown), AKT mutation status (Mutant vs. Wildtype), and metastatic presentation pattern (Bone-only, Liver-involved, Other non-liver visceral).
An absolute standardized mean difference (ASMD) of \(<0.1\) was pre-specified as the threshold for acceptable covariate balance, with the entropy balancing algorithm imposing exact mean balancing constraints during weight estimation. Weighting efficiency was monitored via patient-level weight inspection and effective sample size (ESS) calculation. The primary adjusted HR for the treatment effect was subsequently estimated using a weighted marginal Cox model with Huber-White robust (sandwich) standard errors³ to account for the weighting constraint. The proportional hazards assumption was verified for all Cox models via the Global Schoenfeld residual test. As a secondary analysis to evaluate the independent prognostic impact of baseline clinical characteristics on overall survival under the balancing constraint, an entropy-balanced weighted multivariable Cox model was evaluated. An E-value⁴ was computed to quantify the potential impact of unmeasured confounding on the primary ATE point estimate.
To evaluate clinical drivers of treatment selection, multivariable logistic regression models with Average Marginal Effects (AME) were estimated for cohort assignment (Group B vs. Group A reference groups) and initiation of first-line (1L) chemotherapy. Logistic model discrimination and calibration were evaluated using the Area Under the Receiver Operating Characteristic (ROC) curve and the Hosmer-Lemeshow goodness-of-fit test⁵, respectively. Finally, to contextualize real-world selection patterns, a supplemental cross-sectional analysis of provider preferences was conducted using anonymized audience response system (ARS) data from two recent oncology summits. Descriptive statistics summarized treatment preferences and clinical rationales in response to a standardized high-risk patient vignette. All statistical tests were two-sided with \(p < 0.05\) considered significant. Analyses were performed in R (v. 4.5.2) using the WeightIt and survival packages.
Results
From an initial dataset of 185 patients evaluated, 83 met final eligibility criteria—specifically, being female, aged \(\ge 18\) years, with Stage IV AJCC, HR+/HER2− metastatic breast cancer. This final cohort was divided between two treatment arms: Control Group B (\(N = 53\)) and Intensive Regimen Group A (\(N = 30\)). The median follow-up was 58.9 months (95% CI: 45.6–74.5), with a maximum observation period of 111.7 months in Group B and 166.6 months in Group A; 31 deaths were observed during the study period.
Baseline characteristics were distributed across treatment arms for the majority of sociodemographic and clinical features (Table 1). The cohort was predominantly treated at MSK (51% total; 49% Group B vs. 53% Group A) and was primarily White (77% total; 79% Group B vs. 73% Group A, \(p = 0.7\)). Most patients presented with ductal histology (72% total; 79% Group B vs. 60% Group A, \(p = 0.10\)). Baseline differences were noted in additional sequencing rates, which were performed in 14% of the total cohort, though more frequently in Group A than Group B (27% vs. 7.5%, \(p = 0.024\)).
Per the structural definition of the intensive regimen, 100% of patients in Group A received both an aromatase inhibitor and fulvestrant, compared to 79% (\(p = 0.006\)) and 32% (\(p < 0.001\)) in Group B, respectively. Furthermore, by design, all patients in Group A received at least one targeted therapy. Reflecting this group definition, Group A demonstrated higher rates of exposure to CDK4/6 inhibitors (47% vs. 21%, \(p = 0.026\)) and mTOR inhibitors (70% vs. 17%, \(p < 0.001\)). Group A patients also experienced more total lines of metastatic chemotherapy (median 5.0 vs. 2.0, \(p = 0.004\)), lines of metastatic endocrine therapy (median 4.0 vs. 2.0, \(p < 0.001\)), higher rates of second-line endocrine therapy exposure (83% vs. 55%, \(p = 0.017\)), and a higher total number of systemic metastatic therapies (median 7.5 vs. 4.0, \(p < 0.001\)). Following entropy balancing, covariate balance was achieved, yielding a post-weighting ASMD \(\approx 0\) across all four prespecified factors (age, grade, AKT status, and metastatic site). The weights (range: 0.14–3.04, mean: 1.00) retained an effective sample size (ESS) of 20.1 in the intensive group and 47.3 in the control group.
In the unweighted primary analysis, Group A demonstrated a longer observed median OS than Group B (121.6 months [95% CI: 77.7–NE] vs. 55.9 months [95% CI: 46.4–NE]; HR = 0.35 [95% CI: 0.15–0.81]; \(p = 0.013\); Table 2). The proportional hazards assumption for this unweighted model was satisfied (Global \(p = 0.39\)).
In the primary entropy-balanced model, the survival benefit remained statistically significant. The weighted median OS was 89.9 months (95% CI: 52.3–NE) for Group A versus 55.9 months (95% CI: 46.4–NE) for Group B, yielding a primary adjusted HR of 0.39 (95% CI: 0.16–0.95; \(p = 0.039\)). The proportional hazards assumption for the weighted model was satisfied (Global \(p = 0.32\)). w Quantitative bias analysis yielded an E-value of 3.22 for the primary treatment point estimate, indicating that an unmeasured confounder would require a minimum independent hazard ratio of 3.22 with both treatment allocation and mortality to fully account for the observed therapeutic benefit.
To evaluate the independent prognostic impact of baseline clinical factors, a weighted multivariable Cox model was compared with an unweighted model (Table 3). In the unweighted model, AKT wildtype status was associated with higher risk of mortality (HR = 4.59; 95% CI: 1.07–19.63; \(p = 0.040\)). However, in the entropy-balanced weighted model, AKT wildtype status was no longer significantly associated with survival (weighted HR = 3.07; 95% CI: 0.41–23.26; \(p = 0.278\)). Conversely, baseline liver involvement was identified as an independent predictor of mortality in the weighted model relative to bone-only disease (weighted HR = 2.86; 95% CI: 1.13–7.21; \(p = 0.026\)).
Evaluation of treatment selection drivers indicated allocation patterns based on baseline clinical presentation. Relative to bone-only metastasis, multivariable AME models for Group B control arm assignment (Table 4) demonstrated that patients presenting with liver involvement had a higher absolute probability of being assigned to the control arm (\(+25.2\%\) AME; 95% CI: \(+3.8\%\) to \(+46.5\%\); \(p = 0.021\)). Conversely, presenting with Tumor Grade 3 or Unknown status decreased the probability of assignment to the control arm (\(-23.2\%\) AME; 95% CI: \(-42.9\%\) to \(-3.6\%\); \(p = 0.020\)). This treatment assignment model showed adequate discrimination (ROC = 0.753) and calibration (Hosmer-Lemeshow \(p = 0.534\)).
To investigate clinical mechanisms associated with this therapeutic allocation, drivers of first-line (1L) chemotherapy selection were modeled (Table 5). Baseline liver involvement was associated with frontline chemotherapy choice (ROC = 0.764; Hosmer-Lemeshow \(p = 0.876\)): patients with liver metastases showed a \(43.6\%\) higher absolute probability of receiving 1L chemotherapy compared to those with bone-only disease (95% CI: \(+16.9\%\) to \(+70.3\%\); \(p = 0.001\)). Non-liver visceral site involvement was also associated with a \(30.2\%\) absolute increase in 1L chemotherapy initiation relative to bone-only disease (95% CI: \(+10.5\%\) to \(+49.9\%\); \(p = 0.003\)).
In the cross-sectional provider analysis of ARS data from 112 clinicians across two oncology summits, treatment preferences for a high-risk patient vignette were divided: 47.6% favored an intensive targeted regimen, while 43.7% selected immediate chemotherapy. Half of the participating clinicians (\(N=106\)) identified “visceral crisis involving the liver” as the primary clinical driver of their treatment choice. Among respondents specifically queried about the rationale for prioritizing first-line chemotherapy over targeted options (\(N=112\)), the primary responses were “rapidly progressing disease” (63.4%) and the “presence of symptomatic visceral involvement/crisis” (50.9%).
References
1Hainmueller J. Entropy Balancing for Causal Effects: A Multivariate Reweighting Method to Produce Balanced Samples in Observational Studies. Polit Anal. 2012;20(1):25-46.
2Chen R, Chen G, Yu M. Entropy Balancing for Causal Generalization with Target Sample Summary Information. Biometrics. 2023;79(1):16-28.
3White H. A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity. Econometrica. 1980;48(4):817-838.
4VanderWeele TJ, Ding P. Sensitivity analysis in observational research: introducing the E-value. Ann Intern Med. 2017;167(4):268-274.
5Hosmer DW Jr, Lemeshow S. A goodness-of-fit test for the multiple logistic regression model. Commun Stat Theory Methods. 1980;9(10):1035-1049.