AKT1-Mutated Breast Cancer: Master Analysis
Samuel Baird
2026-06-10
1. Environment Setup & Data Ingestion
This analysis evaluates the real-world overall survival (OS) of patients with AKT1-mutated metastatic breast cancer. This section establishes the R environment, loads necessary statistical packages, 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(EValue)
library(margins)
library(car)
library(ResourceSelection)
library(pROC)
library(ggpubr)
library(patchwork)
library(gt)
library(gtsummary)
library(broom)
library(dplyr)
library(vctrs)
library(DT)
library(scales)
# NOTE: Please update this file path to match your local environment.
tryCatch({
dt <- as.data.table(read_excel(
path = "~/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'.")
})
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")]
ineligible <- dt[eligibility_flag == FALSE]
eligible <- dt[eligibility_flag == TRUE]
cat("Eligible patients:", eligible[, uniqueN(patient_id)], "\n")## Eligible patients: 83# 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)]
# Create first_chemo and first_chemo_binary early so both §9 (AME model)
# and §10 (therapy association tests) can reference them without errors.
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]
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 repetitive per-variable blocks.
Each call:
(a) recodes "NA" strings to "Unknown" when
requested,
(b) prints group distribution, and
(c) runs Fisher’s exact or chi-square association
vs. therapy_flag.
eligible[, met_age_cont := as.numeric(met_age)]
eligible[, seq_int_cont := as.numeric(seq_int)]
eligible[, TMB_cont := as.numeric(TMB)]
eligible[, frac_alt_cont := as.numeric(frac_alt)]
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, B = 2000)
else fisher.test(eligible[[var]], eligible$therapy_flag),
chisq = chisq.test(eligible[[var]], eligible$therapy_flag)
)
print(res)
invisible(res)
}assoc_test("met_biopsy", "fisher", as_factor = TRUE)
assoc_test("HR_status", "fisher", na_to_unknown = FALSE)
assoc_test("bone_mets", "fisher")
assoc_test("lymph_mets", "chisq", as_factor = TRUE)
assoc_test("multiple_mets", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
assoc_test("other_mets", "chisq")
assoc_test("soft_mets", "fisher", as_factor = TRUE)
assoc_test("visceral_mets", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
assoc_test("seq_meth", "fisher", na_to_unknown = FALSE, simulate = TRUE,
seed = 100, as_factor = TRUE)
assoc_test("discordant", "fisher", simulate = TRUE, seed = 101, as_factor = TRUE)
# TMB: continuous -> categorical
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)
# Fraction genome altered: 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)
# Race recoding
eligible[race %in% c("Black", "Other", "NA"), race := "Non-White"]
assoc_test("race", "chisq", na_to_unknown = FALSE)
# Cancer type recoding
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)
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)
# Age at metastatic diagnosis: continuous -> categorical
eligible[, met_age := cut(as.numeric(met_age), c(-Inf, 55, Inf),
labels = c("<=55", ">55"), right = TRUE)]
assoc_test("met_age", "chisq", na_to_unknown = FALSE)
# 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)
# Sequencing interval
eligible[, seq_int := cut(as.numeric(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)
assoc_test("add_seq", "fisher", na_to_unknown = TRUE, as_factor = TRUE)
eligible[, samp_type := as.factor(`Sample Type`)]
assoc_test("samp_type", "fisher", na_to_unknown = FALSE)
assoc_test("boneonly_mets", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
assoc_test("akt", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
eligible[HR_change == "NA" | is.na(HR_change) | HR_change == "Unk", HR_change := "Unknown"]
assoc_test("HR_change", "fisher", na_to_unknown = FALSE, simulate = TRUE,
seed = 106, as_factor = TRUE)
assoc_test("clinical_trial", "chisq", na_to_unknown = FALSE, as_factor = TRUE)
assoc_test("center", "fisher", na_to_unknown = FALSE, simulate = TRUE,
seed = 105, as_factor = TRUE)# Composite stabilized variables for the PS model
eligible[, met_pattern_stable := factor(fcase(
liver_mets == "Yes", "Liver Involved",
boneonly_mets == "Yes", "Bone-Only",
default = "Other Sites (Non-Liver)"),
levels = c("Bone-Only", "Other Sites (Non-Liver)", "Liver Involved"))]
print(eligible[, .N, .(therapy_flag, met_pattern_stable)][
, prop := N / sum(N), by = therapy_flag])## therapy_flag met_pattern_stable N prop
## <fctr> <fctr> <int> <num>
## 1: Group B Liver Involved 13 0.24528302
## 2: Group B Other Sites (Non-Liver) 16 0.30188679
## 3: Group B Bone-Only 24 0.45283019
## 4: Group A Other Sites (Non-Liver) 15 0.50000000
## 5: Group A Bone-Only 13 0.43333333
## 6: Group A Liver Involved 2 0.06666667##
## Fisher's Exact Test for Count Data
##
## data: eligible$met_pattern_stable and eligible$therapy_flag
## p-value = 0.06186
## alternative hypothesis: two.sidedeligible[, tumor_grade_stable := factor(
ifelse(tumor_grade == "Grade 1-2", "Grade 1-2", "Grade 3/Unknown"),
levels = c("Grade 1-2", "Grade 3/Unknown"))]
print(eligible[, .N, .(therapy_flag, tumor_grade_stable)][
, prop := N / sum(N), by = therapy_flag])## therapy_flag tumor_grade_stable N prop
## <fctr> <fctr> <int> <num>
## 1: Group B Grade 1-2 35 0.6603774
## 2: Group B Grade 3/Unknown 18 0.3396226
## 3: Group A Grade 1-2 13 0.4333333
## 4: Group A Grade 3/Unknown 17 0.5666667##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: eligible$tumor_grade_stable and eligible$therapy_flag
## X-squared = 3.1719, df = 1, p-value = 0.074924. 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.
surv_object <- Surv(time = eligible[, survival_time],
event = eligible[, survival_status])
# Overall KM
km_fit_overall <- survfit2(surv_object ~ 1, data = eligible)
print(km_fit_overall)## Call: survfit(formula = surv_object ~ 1, data = eligible)
##
## n events median 0.95LCL 0.95UCL
## [1,] 83 31 77.7 55.9 NA## Call: survfit(formula = surv_object ~ 1, data = eligible)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 12 73 0 1.000 0.0000 1.000 1.000
## 24 62 5 0.929 0.0307 0.871 0.991
## 36 46 7 0.812 0.0494 0.721 0.915
## 48 35 4 0.733 0.0582 0.628 0.857
## 60 24 6 0.601 0.0684 0.480 0.751
## 72 16 3 0.522 0.0731 0.397 0.687## strata median lower upper
## 1 All 77.66447 55.92105 NA## Call: survfit(formula = surv_object ~ therapy_flag, data = eligible)
##
## n events median 0.95LCL 0.95UCL
## therapy_flag=Group B 53 21 55.9 46.4 NA
## therapy_flag=Group A 30 10 121.6 77.7 NA## Call: survfit(formula = surv_object ~ therapy_flag, data = eligible)
##
## therapy_flag=Group B
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 12 44 0 1.000 0.0000 1.000 1.000
## 24 34 5 0.880 0.0503 0.787 0.985
## 36 23 5 0.737 0.0724 0.608 0.894
## 48 16 4 0.596 0.0867 0.448 0.792
## 60 11 3 0.484 0.0913 0.334 0.700
## 72 6 2 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
## 12 29 0 1.000 0.0000 1.000 1.000
## 24 28 0 1.000 0.0000 1.000 1.000
## 36 23 2 0.920 0.0543 0.820 1.000
## 48 19 0 0.920 0.0543 0.820 1.000
## 60 13 3 0.760 0.0952 0.595 0.972
## 72 10 1 0.702 0.1043 0.524 0.939## strata median lower upper
## 1 therapy_flag=Group B 55.92105 46.38158 NA
## 2 therapy_flag=Group A 121.57895 77.66447 NA# Unweighted Cox model and PH test
cox_model_unweighted <- coxph(surv_object ~ therapy_flag, data = eligible)
summary(cox_model_unweighted)## 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
# Overall KM plot
ggsurvplot(km_fit_overall, data = eligible, conf.int = TRUE, risk.table = TRUE,
risk.table.col = "strata", surv.median.line = "hv",
linetype = "solid", break.time.by = 12, xscale = "m_y",
title = "Unweighted Overall OS",
xlab = "Time from Initial Diagnosis (Years)",
ggtheme = theme_pubr())
# Stratified KM 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, xscale = "m_y",
title = "Unweighted Overall Survival",
xlab = "Time from Initial Diagnosis (Years)",
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
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)$quantile
cat("The median follow-up time is:", round(median_follow_up, 1), "months\n")## The median follow-up time is: 58.9 months5. 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", moments = 1, tols = 0.01
)
summary(weightit_model)## Summary of weights
##
## - Weight ranges:
##
## Min Max
## treated 0.48 |---------------------------| 7.975
## control 0.72 |-----| 2.639
##
## - Units with the 5 most extreme weights by group:
##
## 60 56 77 72 33
## treated 3.986 3.986 6.196 7.074 7.975
## 19 31 17 14 8
## control 2.291 2.498 2.639 2.639 2.639
##
## - Weight statistics:
##
## Coef of Var MAD Entropy # Zeros
## treated 0.662 0.502 0.193 0
## control 0.315 0.246 0.047 0
##
## - Effective Sample Sizes:
##
## Control Treated
## Unweighted 53. 30.
## Weighted 48.3 21.07final_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 Measures
## Type Diff.Un Diff.Adj
## met_age_>55 Binary 0.1239 0.01
## tumor_grade_stable_Grade 3/Unknown Binary 0.2270 0.01
## akt_Wildtype Binary 0.1302 0.01
## met_pattern_stable_Bone-Only Binary 0.0195 0.01
## met_pattern_stable_Other Sites (Non-Liver) Binary 0.1981 0.01
## met_pattern_stable_Liver Involved Binary 0.1786 0.02
## 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.02 Balanced, <0.1
##
## Effective sample sizes
## Group B Group A
## Unadjusted 53. 30.
## Adjusted 48.3 21.07eligible[, ws := weightit_model$weights]
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
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 166 59.8 81.3 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 145.0 0.00 1.000 0.0000 1.000 1.000
## 24 124.8 6.91 0.950 0.0226 0.907 0.996
## 36 99.4 12.81 0.846 0.0484 0.756 0.946
## 48 76.9 7.35 0.776 0.0578 0.671 0.898
## 60 45.8 15.75 0.595 0.0828 0.453 0.782
## 72 33.0 4.15 0.539 0.0843 0.396 0.732## strata median lower upper
## 1 All 81.28289 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 83 32.9 55.9 46.4 NA
## therapy_flag=Group A 30 83 26.9 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 83.00 0.00 1.000 0.0000 1.000 1.000
## 12 69.94 0.00 1.000 0.0000 1.000 1.000
## 24 53.74 6.91 0.896 0.0454 0.811 0.989
## 36 36.48 8.35 0.742 0.0768 0.605 0.908
## 48 25.38 7.35 0.580 0.0930 0.423 0.794
## 60 17.04 5.04 0.464 0.0957 0.310 0.696
## 72 9.59 2.53 0.384 0.0995 0.231 0.638
##
## therapy_flag=Group A
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 83.0 0.00 1.000 0.0000 1.000 1.000
## 12 75.0 0.00 1.000 0.0000 1.000 1.000
## 24 71.0 0.00 1.000 0.0000 1.000 1.000
## 36 62.9 4.46 0.934 0.0575 0.828 1.000
## 48 51.6 0.00 0.934 0.0575 0.828 1.000
## 60 28.8 10.70 0.698 0.1342 0.478 1.000
## 72 23.5 1.62 0.658 0.1345 0.441 0.982## 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.9781 0.3760 0.2828 0.4554 -2.148 0.0317 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A 0.376 2.659 0.154 0.918
##
## Concordance= 0.647 (se = 0.054 )
## Likelihood ratio test= 12.37 on 1 df, p=4e-04
## Wald test = 4.61 on 1 df, p=0.03
## Score (logrank) test = 12.87 on 1 df, p=3e-04, Robust = 5.7 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 2.53 1 0.11
## GLOBAL 2.53 1 0.11# Overall weighted KM plot
ggsurvplot(km_fit_weighted_overall, data = eligible, conf.int = TRUE,
risk.table = TRUE, risk.table.col = "strata",
surv.median.line = "hv", linetype = "solid",
break.time.by = 12, xscale = "m_y",
title = "Weighted OS",
xlab = "Time from Initial Diagnosis (Years)",
ggtheme = theme_pubr())
# Stratified weighted KM plot with unweighted risk table
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, xscale = "m_y",
title = "Adjusted Overall Survival (Entropy Balanced)",
xlab = "Time from Initial Diagnosis (Years)",
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
# Proportional hazards assumption: numerical + graphical
cox.zph(cox_model_weighted, transform = "km", terms = TRUE)## chisq df p
## therapy_flag 2.53 1 0.11
## GLOBAL 2.53 1 0.11ggcoxzph(cox.zph(cox_model_weighted), 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")
## point lower upper
## RR 0.5150192 0.2826446 0.9438471
## E-values 3.2938690 NA 1.31055777. 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.
covariates_filtered <- c("therapy_flag", "met_age", "tumor_grade_stable",
"akt", "met_pattern_stable")
full_label_map <- c(
"therapy_flagGroup A" = "Intensive Regimen (Group A)",
"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)"
)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 2: 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 3: 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 4: 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 5: 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.
8.1 Dynamic Landmark Exploratory Sweep (3–48 Months)
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)
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]
n_per_arm <- d_land[, .N, keyby = therapy_flag]$N
events_per_arm <- d_land[, sum(survival_status), keyby = therapy_flag]$V1
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)
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
)
}))
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 6: 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.2 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.
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) / 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)
fu_summary <- data.table(
Metric = c("Median Follow-up", "95% Lower Bound", "95% Upper Bound"),
`Months (from Landmark)` = 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 7: Follow-up Time Summary from', landmark_time, '-Month Landmark')
)
)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 (Years)",
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 (Years)")
pp
}
# Unweighted 9-month 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)
# 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.543 |---------------------------| 6.237
## control 0.587 |-------------| 3.491
##
## - Units with the 5 most extreme weights by group:
##
## 52 48 69 64 26
## treated 3.798 3.798 4.763 6.103 6.237
## 41 13 24 9 4
## control 2.545 2.545 2.92 3.491 3.491
##
## - Weight statistics:
##
## Coef of Var MAD Entropy # Zeros
## treated 0.58 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 75 33.7 46.9 37.4 NA
## therapy_flag=Group A 30 75 23.8 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 75.00 0.00 1.000 0.0000 1.000 1.000
## 12 63.16 4.23 0.942 0.0338 0.878 1.000
## 24 44.09 8.50 0.800 0.0651 0.683 0.939
## 36 29.91 6.91 0.661 0.0895 0.507 0.862
## 48 18.89 8.91 0.464 0.0981 0.307 0.702
## 60 11.38 2.42 0.389 0.1018 0.233 0.650
## 72 7.66 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 75.0 0.000 1.000 0.0000 1.000 1.000
## 12 65.0 0.000 1.000 0.0000 1.000 1.000
## 24 57.8 3.366 0.945 0.0534 0.846 1.000
## 36 54.9 0.543 0.936 0.0540 0.836 1.000
## 48 29.4 8.976 0.721 0.1231 0.516 1.000
## 60 25.9 1.498 0.682 0.1247 0.476 0.976
## 72 17.1 1.690 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.0319 0.3563 0.2949 0.4539 -2.273 0.023 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A 0.3563 2.806 0.1464 0.8674
##
## Concordance= 0.65 (se = 0.051 )
## Likelihood ratio test= 13 on 1 df, p=3e-04
## Wald test = 5.17 on 1 df, p=0.02
## Score (logrank) test = 13.28 on 1 df, p=3e-04, Robust = 6.01 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).## chisq df p
## therapy_flag 2.29 1 0.13
## GLOBAL 2.29 1 0.13p_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 (Years)") +
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.
ame_point_size <- 3.5
ame_errorbar_width <- 0.15
ame_errorbar_lw <- 0.7
ame_vline_color <- "gray60"
ame_break_step <- 0.2
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_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))
}
run_ame_panel <- function(target_group, level_order, label_map,
title_txt, x_lab, ax = NULL, roc_title,
subtitle_txt = NULL, caption_txt = NULL) {
eligible[, therapy_flag_binary := ifelse(therapy_flag == target_group, 1, 0)]
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)]
fit <- glm(therapy_flag_binary ~ met_age + tumor_grade_stable +
akt + met_pattern_stable,
data = eligible, family = binomial())
print(summary(fit))
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))
}
ame_tbl <- as.data.table(summary(margins(fit)))
ame_tbl[, factor := ifelse(factor %in% names(label_map),
label_map[factor], factor)]
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
ame_table_b <- panel_b$ame
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 8: Average Marginal Effects for Group B Assignment'))1L Chemotherapy Selection Drivers
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.
# Apply 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)"))]
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: 5vif_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.015017hoslem_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_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 9: Average Marginal Effects for 1L Chemotherapy Selection'))10. Therapy Usage Association Tests
This section evaluates associations between treatment group assignment and detailed therapy usage patterns across chemotherapy lines, endocrine therapy sequences, targeted agents, and treatment durations.
# --- Chemotherapy ---
# Note: first_chemo and first_chemo_binary were created in §2 (cohort_definition)
# to ensure availability for the §9 AME model. Association test runs directly here.
first_chemo_dis <- eligible[, .N, .(therapy_flag, first_chemo)][
, prop := N/sum(N), therapy_flag]
print(first_chemo_dis)
chisq.test(eligible[, first_chemo], eligible[, therapy_flag])
second_chemo_dis <- eligible[, .N, .(therapy_flag, `Chemotherapy in second-line treatment`)][
, prop := N/sum(N), therapy_flag]
print(second_chemo_dis)
eligible[, second_chemo := as.factor(`Chemotherapy in second-line treatment`)]
eligible[, `Chemotherapy in second-line treatment` := NULL]
fisher.test(eligible[, second_chemo], eligible[, therapy_flag])
eligible[, chemo_overall := as.factor(`Chemotherapy Overall`)]
eligible[, `Chemotherapy Overall` := NULL]
chisq.test(eligible[, chemo_overall], eligible[, therapy_flag])
eligible[, met_chemo := as.numeric(as.factor(`Chemotherapy lines received in Metastatic Treatment`))]
eligible[, `Chemotherapy lines received in Metastatic Treatment` := NULL]
eligible[therapy_flag == "Group A", summary(met_chemo)]
eligible[therapy_flag == "Group B", summary(met_chemo)]
wilcox.test(eligible[, met_chemo] ~ eligible[, therapy_flag], conf.int = TRUE)
# --- Endocrine Therapy ---
eligible[, prior_endo := as.factor(`Total Number of Prior Endocrine Therapies`)]
eligible[, `Total Number of Prior Endocrine Therapies` := NULL]
set.seed(105)
fisher.test(eligible[, prior_endo], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[, met_endo := as.numeric(`Endocrine Therapies received in Metastatic Treatment`)]
eligible[, `Endocrine Therapies received in Metastatic Treatment` := NULL]
wilcox.test(eligible[, met_endo] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[, met_endo_cat := cut(met_endo, c(-Inf, 1, 2, 3, Inf),
labels = c("0", "1", "2", "3+"), right = FALSE)]
set.seed(106)
fisher.test(eligible[, met_endo_cat], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[`Overall Endocrine Therapy Sensitivity` == "NA",
`Overall Endocrine Therapy Sensitivity` := "Unknown"]
eligible[, endo_sens := as.factor(`Overall Endocrine Therapy Sensitivity`)]
eligible[, `Overall Endocrine Therapy Sensitivity` := NULL]
set.seed(107)
fisher.test(eligible[, endo_sens], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[, first_endo := as.factor(`Endocrine Therapy in first-line treatment?`)]
eligible[, `Endocrine Therapy in first-line treatment?` := NULL]
chisq.test(eligible[, first_endo], eligible[, therapy_flag])
eligible[, second_endo := as.factor(`Endocrine Therapy in second-line treatment?`)]
eligible[, `Endocrine Therapy in second-line treatment?` := NULL]
fisher.test(eligible[, second_endo], eligible[, therapy_flag])
# --- Fulvestrant ---
eligible[, fulvestrant := as.factor(fulvestrant)]
set.seed(108)
fisher.test(eligible[, fulvestrant], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[, first_fulv := as.factor(`Fulvestrant in first-line treatment`)]
eligible[, `Fulvestrant in first-line treatment` := NULL]
set.seed(109)
fisher.test(eligible[, first_fulv], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[, second_fulv := as.factor(`Fulvestrant in second-line treatment`)]
eligible[, `Fulvestrant in second-line treatment` := NULL]
fisher.test(eligible[, second_fulv], eligible[, therapy_flag])
eligible[, dur_fulv := as.numeric(`Overall Duration of Fulvestrant Treatment (any line) (Months)`)]
eligible[, `Overall Duration of Fulvestrant Treatment (any line) (Months)` := NULL]
wilcox.test(eligible[, dur_fulv] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[`Overall Fulvestrant Censoring status (any line)` == "NA",
`Overall Fulvestrant Censoring status (any line)` := "Unknown"]
eligible[, cen_fulv := as.factor(`Overall Fulvestrant Censoring status (any line)`)]
eligible[, `Overall Fulvestrant Censoring status (any line)` := NULL]
set.seed(110)
fisher.test(eligible[, cen_fulv], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[, tamoxifen := as.factor(`Tamoxifen Overall`)]
eligible[, `Tamoxifen Overall` := NULL]
chisq.test(eligible[, tamoxifen], eligible[, therapy_flag])
eligible[, aromatase_inhibitor := as.factor(aromatase_inhibitor)]
set.seed(111)
fisher.test(eligible[, aromatase_inhibitor], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[`Ovarian function suppression Overall` == "NA",
`Ovarian function suppression Overall` := "Unknown"]
eligible[, ovarian_sup := as.factor(`Ovarian function suppression Overall`)]
eligible[, `Ovarian function suppression Overall` := NULL]
set.seed(112)
fisher.test(eligible[, ovarian_sup], eligible[, therapy_flag], simulate.p.value = TRUE)
# --- Targeted Therapies ---
eligible[, AKT_inhibitor := as.factor(AKT_inhibitor)]
fisher.test(eligible[, AKT_inhibitor], eligible[, therapy_flag])
eligible[, cdk46_inhibitor := as.factor(`CDK4/6`)]
eligible[, `CDK4/6` := NULL]
chisq.test(eligible[, cdk46_inhibitor], eligible[, therapy_flag])
eligible[, cdk46_dur := as.numeric(`Duration of CDK4/6 treatment (any line) (Months)`)]
eligible[, `Duration of CDK4/6 treatment (any line) (Months)` := NULL]
wilcox.test(eligible[, cdk46_dur] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[`CDK4/6 Censoring Status (any line)` == "NA",
`CDK4/6 Censoring Status (any line)` := "Unknown"]
eligible[, cdk46_censoring := as.factor(`CDK4/6 Censoring Status (any line)`)]
eligible[, `CDK4/6 Censoring Status (any line)` := NULL]
set.seed(113)
fisher.test(eligible[, cdk46_censoring], eligible[, therapy_flag], simulate.p.value = TRUE)
chisq.test(eligible[, mTOR], eligible[, therapy_flag])
eligible[, mTOR_dur := as.numeric(`mTOR Duration of Treatment (any line)`)]
eligible[, `mTOR Duration of Treatment (any line)` := NULL]
wilcox.test(eligible[, mTOR_dur] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[`Overall mTOR Censoring Status (any line) (Months)` == "NA",
`Overall mTOR Censoring Status (any line) (Months)` := "Unknown"]
eligible[, mTOR_censoring := as.factor(`Overall mTOR Censoring Status (any line) (Months)`)]
eligible[, `Overall mTOR Censoring Status (any line) (Months)` := NULL]
set.seed(114)
fisher.test(eligible[, mTOR_censoring], eligible[, therapy_flag], simulate.p.value = TRUE)
# --- Treatment Lines Summary ---
eligible[, met_therapies := as.numeric(`Total number of therapies received in metastatic disease treatment`)]
eligible[, `Total number of therapies received in metastatic disease treatment` := NULL]
wilcox.test(eligible[, met_therapies] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[, first_dur := as.numeric(`Duration of first-line treatment (Months)`)]
eligible[, `Duration of first-line treatment (Months)` := NULL]
wilcox.test(eligible[, first_dur] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[`First-line treatment Censoring Status` == "NA",
`First-line treatment Censoring Status` := "Unknown"]
eligible[, first_cen := as.factor(`First-line treatment Censoring Status`)]
eligible[, `First-line treatment Censoring Status` := NULL]
set.seed(115)
fisher.test(eligible[, first_cen], eligible[, therapy_flag], simulate.p.value = TRUE)
eligible[, second_dur := as.numeric(`Duration of second-line treatment (Months)`)]
eligible[, `Duration of second-line treatment (Months)` := NULL]
wilcox.test(eligible[, second_dur] ~ eligible[, therapy_flag], conf.int = TRUE)
eligible[`Second-line treatment Censoring Status` == "NA",
`Second-line treatment Censoring Status` := "Unknown"]
eligible[, second_cen := as.factor(`Second-line treatment Censoring Status`)]
eligible[, `Second-line treatment Censoring Status` := NULL]
set.seed(116)
fisher.test(eligible[, second_cen], eligible[, therapy_flag], simulate.p.value = TRUE)
# 1L chemo driver validation (first_chemo_binary already created in §2)
eligible[, first_chemo_num := ifelse(first_chemo == "Yes", 1, 0)]
fisher.test(table(eligible$liver_mets, eligible$first_chemo_num))11. Publication Tables
This section generates all publication-ready formatted tables and
exports them as .docx files. Tables cover patient
demographics, univariate and multivariate Cox results, primary outcomes,
and Average Marginal Effects.
# --- Table 1: Demographics ---
table_labels <- list(
met_age_cont ~ "Age at Metastatic Diagnosis (Continuous, Years)",
met_age ~ "Age at Metastatic Diagnosis (Categorical)",
race ~ "Racial Cohort",
cancer_type ~ "Histological Cancer Type",
tumor_grade ~ "Tumor Grade (Initial)",
HR_status ~ "Hormone Receptor (HR) Status",
akt ~ "AKT Mutation Status",
TMB_cont ~ "Tumor Mutational Burden (Continuous)",
TMB ~ "Tumor Mutational Burden (Categorical)",
frac_alt_cont ~ "Fraction of Genome Altered (Continuous)",
frac_alt ~ "Fraction of Genome Altered (Categorical)",
samp_type ~ "Biopsy Sample Location Type",
met_biopsy ~ "Metastatic Biopsy Performed",
discordant ~ "Biopsy Receptor Discordance",
seq_meth ~ "Next-Gen Sequencing Method",
seq_int_cont ~ "Interval to Sequencing (Continuous, Months)",
seq_int ~ "Interval to Sequencing (Categorical)",
clinical_trial ~ "Enrolled in Clinical Trial",
bone_mets ~ "Bone Metastases",
boneonly_mets ~ "Bone-Only Metastatic Disease",
liver_mets ~ "Liver Metastases",
lung_mets ~ "Lung Metastases",
brain_mets ~ "Brain Metastases",
lymph_mets ~ "Lymph Node Metastases",
soft_mets ~ "Soft Tissue Metastases",
visceral_mets ~ "Visceral Metastasis Involvement",
multiple_mets ~ "Multiple Metastatic Sites Present",
met_pattern_stable ~ "Metastatic Pattern Profile (PS Model)",
tumor_grade_stable ~ "Tumor Grade Risk Status (PS Model)",
first_chemo ~ "Chemotherapy in First-Line (1L) Regimen",
second_chemo ~ "Chemotherapy in Second-Line (2L) Regimen",
chemo_overall ~ "Ever Received Chemotherapy",
met_chemo ~ "Total Lines of Metastatic Chemotherapy",
first_endo ~ "Endocrine Therapy in First-Line (1L) Regimen",
second_endo ~ "Endocrine Therapy in Second-Line (2L) Regimen",
prior_endo ~ "Number of Prior Endocrine Therapy Lines",
met_endo ~ "Metastatic Endocrine Therapies Received",
met_endo_cat ~ "Metastatic Endocrine Lines Grouped",
endo_sens ~ "Overall Endocrine Therapy Sensitivity",
tamoxifen ~ "Ever Received Tamoxifen Overall",
aromatase_inhibitor~ "Ever Received Aromatase Inhibitor",
fulvestrant ~ "Ever Received Fulvestrant Overall",
first_fulv ~ "Fulvestrant in First-Line (1L) Regimen",
second_fulv ~ "Fulvestrant in Second-Line (2L) Regimen",
dur_fulv ~ "Duration of Fulvestrant Therapy (Months)",
cen_fulv ~ "Fulvestrant Treatment Censoring Status",
ovarian_sup ~ "Ovarian Function Suppression (OFS) Received",
AKT_inhibitor ~ "Ever Received AKT Inhibitor",
cdk46_inhibitor ~ "Ever Received CDK4/6 Inhibitor",
cdk46_dur ~ "Duration of CDK4/6 Inhibitor Therapy (Months)",
cdk46_censoring ~ "CDK4/6 Inhibitor Censoring Status",
mTOR ~ "Ever Received mTOR Inhibitor",
mTOR_dur ~ "Duration of mTOR Inhibitor Therapy (Months)",
mTOR_censoring ~ "mTOR Inhibitor Censoring Status",
met_therapies ~ "Total Lines of Systemic Metastatic Therapies",
first_dur ~ "Duration of First-Line (1L) Treatment (Months)",
first_cen ~ "First-Line (1L) Treatment Censoring Status",
second_dur ~ "Duration of Second-Line (2L) Treatment (Months)",
second_cen ~ "Second-Line (2L) Treatment Censoring Status"
)
fisher_vars <- c("center", "HR_change", "add_seq", "met_biopsy", "HR_status",
"bone_mets", "soft_mets", "brain_mets", "liver_mets",
"tumor_grade", "prior_endo", "met_endo_cat", "endo_sens",
"second_chemo", "fulvestrant", "first_fulv", "second_fulv",
"cen_fulv", "aromatase_inhibitor", "ovarian_sup",
"AKT_inhibitor", "cdk46_censoring", "mTOR_censoring",
"first_cen", "second_cen")
all_vars <- c("center", "met_age_cont", "met_age", "race", "cancer_type",
"tumor_grade", "HR_status", "HR_change", "add_seq", "akt",
"TMB_cont", "TMB", "frac_alt_cont", "frac_alt", "samp_type",
"met_biopsy", "discordant", "seq_meth", "seq_int_cont",
"seq_int", "clinical_trial", "bone_mets", "boneonly_mets",
"liver_mets", "lung_mets", "brain_mets", "lymph_mets",
"soft_mets", "visceral_mets", "multiple_mets",
"met_pattern_stable", "tumor_grade_stable", "first_chemo",
"second_chemo", "chemo_overall", "met_chemo", "first_endo",
"second_endo", "prior_endo", "met_endo", "met_endo_cat",
"endo_sens", "tamoxifen", "aromatase_inhibitor", "fulvestrant",
"first_fulv", "second_fulv", "dur_fulv", "cen_fulv",
"ovarian_sup", "AKT_inhibitor", "cdk46_inhibitor", "cdk46_dur",
"cdk46_censoring", "mTOR", "mTOR_dur", "mTOR_censoring",
"met_therapies", "first_dur", "first_cen", "second_dur",
"second_cen")
continuous_vars <- c("met_age_cont", "TMB_cont", "frac_alt_cont",
"seq_int_cont", "met_chemo", "met_endo", "dur_fulv",
"cdk46_dur", "mTOR_dur", "met_therapies",
"first_dur", "second_dur")
categorical_vars <- setdiff(all_vars, continuous_vars)
chisq_vars <- setdiff(categorical_vars, fisher_vars)
pretty_demographics <- eligible |>
tbl_summary(
by = therapy_flag, include = all_of(all_vars),
label = c(list(center = "Center", HR_change = "HR Change",
add_seq = "Additional Sequencing"), table_labels),
type = list(all_of(continuous_vars) ~ "continuous"),
statistic = list(all_categorical() ~ "{n} ({p}%)",
all_continuous() ~ "{median} ({p25}, {p75})"),
digits = all_continuous() ~ 1, missing = "no"
) |>
add_overall(col_label = "**Total Cohort**<br>N = {N}") |>
add_p(test = list(all_of(continuous_vars) ~ "wilcox.test",
all_of(fisher_vars) ~ "fisher.test",
all_of(chisq_vars) ~ "chisq.test"),
test.args = list(
all_tests("fisher.test") ~ list(simulate.p.value = TRUE, B = 2000))) |>
bold_labels() |> as_gt() |>
tab_header(title = "Table 1. Demographic and Clinical Characteristics of the Patient Cohort by Treatment Allocation") |>
cols_align(align = "left", columns = label) |>
cols_align(align = "center", columns = c(starts_with("stat_"), "p.value")) |>
tab_options(
table.border.top.color = "black", table.border.top.width = px(2),
table.border.bottom.color = "black", table.border.bottom.width = px(2),
column_labels.border.bottom.color = "black",
column_labels.border.bottom.width = px(1.5),
table_body.border.bottom.color = "black",
table_body.border.bottom.width = px(1.5),
row_group.border.bottom.color = "gray80",
data_row.padding = px(4), table.font.size = px(11))
pretty_demographics| Table 1. Demographic and Clinical Characteristics of the Patient Cohort by Treatment Allocation | ||||
| Characteristic | Total Cohort N = 831 |
Group B N = 531 |
Group A N = 301 |
p-value2 |
|---|---|---|---|---|
| Center | 0.7 | |||
| Â Â Â Â DFCI | 17 (20%) | 10 (19%) | 7 (23%) | |
| Â Â Â Â GRCC | 2 (2.4%) | 2 (3.8%) | 0 (0%) | |
| Â Â Â Â MDA | 13 (16%) | 9 (17%) | 4 (13%) | |
| Â Â Â Â MSK | 42 (51%) | 26 (49%) | 16 (53%) | |
| Â Â Â Â UHN | 1 (1.2%) | 0 (0%) | 1 (3.3%) | |
| Â Â Â Â VICC | 8 (9.6%) | 6 (11%) | 2 (6.7%) | |
| Age at Metastatic Diagnosis (Continuous, Years) | 55.0 (48.0, 62.0) | 55.0 (49.0, 63.0) | 55.0 (48.0, 59.0) | 0.3 |
| Age at Metastatic Diagnosis (Categorical) | 0.4 | |||
| Â Â Â Â <=55 | 46 (55%) | 27 (51%) | 19 (63%) | |
| Â Â Â Â >55 | 37 (45%) | 26 (49%) | 11 (37%) | |
| Racial Cohort | 0.7 | |||
| Â Â Â Â Non-White | 19 (23%) | 11 (21%) | 8 (27%) | |
| Â Â Â Â White | 64 (77%) | 42 (79%) | 22 (73%) | |
| Histological Cancer Type | 0.10 | |||
| Â Â Â Â Ductal | 60 (72%) | 42 (79%) | 18 (60%) | |
| Â Â Â Â Non-Ductal | 23 (28%) | 11 (21%) | 12 (40%) | |
| Tumor Grade (Initial) | 0.12 | |||
| Â Â Â Â Grade 1-2 | 48 (58%) | 35 (66%) | 13 (43%) | |
| Â Â Â Â Grade 3 | 26 (31%) | 14 (26%) | 12 (40%) | |
| Â Â Â Â Unknown | 9 (11%) | 4 (7.5%) | 5 (17%) | |
| Hormone Receptor (HR) Status | >0.9 | |||
| Â Â Â Â ER+ & PR- | 6 (7.2%) | 4 (7.5%) | 2 (6.7%) | |
| Â Â Â Â ER+ & PR+ | 77 (93%) | 49 (92%) | 28 (93%) | |
| HR Change | 0.7 | |||
| Â Â Â Â ER+ & PR- | 5 (6.0%) | 3 (5.7%) | 2 (6.7%) | |
| Â Â Â Â ER+ & PR+ | 6 (7.2%) | 5 (9.4%) | 1 (3.3%) | |
| Â Â Â Â Unknown | 72 (87%) | 45 (85%) | 27 (90%) | |
| Additional Sequencing | 12 (14%) | 4 (7.5%) | 8 (27%) | 0.024 |
| AKT Mutation Status | 0.3 | |||
| Â Â Â Â Mutant | 18 (22%) | 9 (17%) | 9 (30%) | |
| Â Â Â Â Wildtype | 65 (78%) | 44 (83%) | 21 (70%) | |
| Tumor Mutational Burden (Continuous) | 4.0 (2.0, 6.0) | 4.0 (1.0, 6.0) | 4.0 (2.0, 6.0) | 0.5 |
| Tumor Mutational Burden (Categorical) | >0.9 | |||
| Â Â Â Â low | 75 (90%) | 48 (91%) | 27 (90%) | |
| Â Â Â Â high | 3 (3.6%) | 2 (3.8%) | 1 (3.3%) | |
| Â Â Â Â Unknown | 5 (6.0%) | 3 (5.7%) | 2 (6.7%) | |
| Fraction of Genome Altered (Continuous) | 0.2 (0.1, 0.3) | 0.3 (0.2, 0.3) | 0.2 (0.1, 0.3) | 0.073 |
| Fraction of Genome Altered (Categorical) | 0.4 | |||
| Â Â Â Â low fra | 36 (43%) | 21 (40%) | 15 (50%) | |
| Â Â Â Â high fra | 7 (8.4%) | 6 (11%) | 1 (3.3%) | |
| Â Â Â Â Unknown | 40 (48%) | 26 (49%) | 14 (47%) | |
| Biopsy Sample Location Type | 0.7 | |||
| Â Â Â Â LRR | 1 (1.2%) | 1 (1.9%) | 0 (0%) | |
| Â Â Â Â Metastatic | 32 (39%) | 19 (36%) | 13 (43%) | |
| Â Â Â Â Primary | 49 (59%) | 32 (60%) | 17 (57%) | |
| Â Â Â Â Unk | 1 (1.2%) | 1 (1.9%) | 0 (0%) | |
| Metastatic Biopsy Performed | 0.9 | |||
| Â Â Â Â No | 26 (31%) | 16 (30%) | 10 (33%) | |
| Â Â Â Â Unknown | 5 (6.0%) | 4 (7.5%) | 1 (3.3%) | |
| Â Â Â Â Yes | 52 (63%) | 33 (62%) | 19 (63%) | |
| Biopsy Receptor Discordance | >0.9 | |||
| Â Â Â Â No | 38 (46%) | 24 (45%) | 14 (47%) | |
| Â Â Â Â Unknown | 32 (39%) | 21 (40%) | 11 (37%) | |
| Â Â Â Â Yes | 13 (16%) | 8 (15%) | 5 (17%) | |
| Next-Gen Sequencing Method | 0.7 | |||
| Â Â Â Â Hiseq | 40 (48%) | 25 (47%) | 15 (50%) | |
| Â Â Â Â ion-proton/PGM/torrent | 15 (18%) | 11 (21%) | 4 (13%) | |
| Â Â Â Â Miseq | 1 (1.2%) | 0 (0%) | 1 (3.3%) | |
| Â Â Â Â MSK-impact | 22 (27%) | 14 (26%) | 8 (27%) | |
| Â Â Â Â Sequenom | 5 (6.0%) | 3 (5.7%) | 2 (6.7%) | |
| Interval to Sequencing (Continuous, Months) | 3.0 (1.0, 6.0) | 2.0 (1.0, 4.0) | 5.0 (3.0, 7.0) | 0.006 |
| Interval to Sequencing (Categorical) | 0.084 | |||
| Â Â Â Â seq int 1-3 months | 45 (54%) | 33 (62%) | 12 (40%) | |
| Â Â Â Â seq int 4+ months | 38 (46%) | 20 (38%) | 18 (60%) | |
| Enrolled in Clinical Trial | 25 (30%) | 14 (26%) | 11 (37%) | 0.5 |
| Bone Metastases | 0.12 | |||
| Â Â Â Â No | 7 (8.4%) | 4 (7.5%) | 3 (10%) | |
| Â Â Â Â Unknown | 9 (11%) | 3 (5.7%) | 6 (20%) | |
| Â Â Â Â Yes | 67 (81%) | 46 (87%) | 21 (70%) | |
| Bone-Only Metastatic Disease | 37 (45%) | 24 (45%) | 13 (43%) | >0.9 |
| Liver Metastases | 0.11 | |||
| Â Â Â Â No | 49 (59%) | 30 (57%) | 19 (63%) | |
| Â Â Â Â Unknown | 19 (23%) | 10 (19%) | 9 (30%) | |
| Â Â Â Â Yes | 15 (18%) | 13 (25%) | 2 (6.7%) | |
| Lung Metastases | 0.2 | |||
| Â Â Â Â No | 52 (63%) | 36 (68%) | 16 (53%) | |
| Â Â Â Â Unknown | 17 (20%) | 11 (21%) | 6 (20%) | |
| Â Â Â Â Yes | 14 (17%) | 6 (11%) | 8 (27%) | |
| Brain Metastases | 0.8 | |||
| Â Â Â Â No | 56 (67%) | 37 (70%) | 19 (63%) | |
| Â Â Â Â Unknown | 24 (29%) | 14 (26%) | 10 (33%) | |
| Â Â Â Â Yes | 3 (3.6%) | 2 (3.8%) | 1 (3.3%) | |
| Lymph Node Metastases | 0.7 | |||
| Â Â Â Â No | 42 (51%) | 27 (51%) | 15 (50%) | |
| Â Â Â Â Unknown | 16 (19%) | 9 (17%) | 7 (23%) | |
| Â Â Â Â Yes | 25 (30%) | 17 (32%) | 8 (27%) | |
| Soft Tissue Metastases | 0.4 | |||
| Â Â Â Â No | 54 (65%) | 37 (70%) | 17 (57%) | |
| Â Â Â Â Unknown | 22 (27%) | 12 (23%) | 10 (33%) | |
| Â Â Â Â Yes | 7 (8.4%) | 4 (7.5%) | 3 (10%) | |
| Visceral Metastasis Involvement | 25 (30%) | 17 (32%) | 8 (27%) | 0.8 |
| Multiple Metastatic Sites Present | 34 (41%) | 23 (43%) | 11 (37%) | 0.7 |
| Metastatic Pattern Profile (PS Model) | 0.067 | |||
| Â Â Â Â Bone-Only | 37 (45%) | 24 (45%) | 13 (43%) | |
| Â Â Â Â Liver Involved | 15 (18%) | 13 (25%) | 2 (6.7%) | |
| Â Â Â Â Other Sites (Non-Liver) | 31 (37%) | 16 (30%) | 15 (50%) | |
| Tumor Grade Risk Status (PS Model) | 0.075 | |||
| Â Â Â Â Grade 1-2 | 48 (58%) | 35 (66%) | 13 (43%) | |
| Â Â Â Â Grade 3/Unknown | 35 (42%) | 18 (34%) | 17 (57%) | |
| Chemotherapy in First-Line (1L) Regimen | 23 (28%) | 15 (28%) | 8 (27%) | >0.9 |
| Chemotherapy in Second-Line (2L) Regimen | 16 (19%) | 11 (21%) | 5 (17%) | 0.8 |
| Ever Received Chemotherapy | 59 (71%) | 35 (66%) | 24 (80%) | 0.3 |
| Total Lines of Metastatic Chemotherapy | 3.0 (1.0, 5.0) | 2.0 (1.0, 4.0) | 5.0 (2.0, 6.0) | 0.004 |
| Endocrine Therapy in First-Line (1L) Regimen | 59 (71%) | 38 (72%) | 21 (70%) | >0.9 |
| Endocrine Therapy in Second-Line (2L) Regimen | 54 (65%) | 29 (55%) | 25 (83%) | 0.017 |
| Number of Prior Endocrine Therapy Lines | <0.001 | |||
| Â Â Â Â >=3 | 48 (58%) | 21 (40%) | 27 (90%) | |
| Â Â Â Â 1 | 18 (22%) | 18 (34%) | 0 (0%) | |
| Â Â Â Â 2 | 14 (17%) | 11 (21%) | 3 (10%) | |
| Â Â Â Â No | 3 (3.6%) | 3 (5.7%) | 0 (0%) | |
| Metastatic Endocrine Therapies Received | 3.0 (1.0, 4.0) | 2.0 (1.0, 3.0) | 4.0 (3.0, 5.0) | <0.001 |
| Metastatic Endocrine Lines Grouped | <0.001 | |||
| Â Â Â Â 0 | 3 (3.6%) | 3 (5.7%) | 0 (0%) | |
| Â Â Â Â 1 | 18 (22%) | 18 (34%) | 0 (0%) | |
| Â Â Â Â 2 | 14 (17%) | 11 (21%) | 3 (10%) | |
| Â Â Â Â 3+ | 48 (58%) | 21 (40%) | 27 (90%) | |
| Overall Endocrine Therapy Sensitivity | 0.026 | |||
| Â Â Â Â No | 10 (12%) | 8 (15%) | 2 (6.7%) | |
| Â Â Â Â Unknown | 8 (9.6%) | 8 (15%) | 0 (0%) | |
| Â Â Â Â Yes | 65 (78%) | 37 (70%) | 28 (93%) | |
| Ever Received Tamoxifen Overall | 33 (40%) | 21 (40%) | 12 (40%) | >0.9 |
| Ever Received Aromatase Inhibitor | 72 (87%) | 42 (79%) | 30 (100%) | 0.006 |
| Ever Received Fulvestrant Overall | 47 (57%) | 17 (32%) | 30 (100%) | <0.001 |
| Fulvestrant in First-Line (1L) Regimen | 6 (7.2%) | 2 (3.8%) | 4 (13%) | 0.2 |
| Fulvestrant in Second-Line (2L) Regimen | 11 (13%) | 5 (9.4%) | 6 (20%) | 0.2 |
| Duration of Fulvestrant Therapy (Months) | 4.6 (2.8, 8.8) | 3.9 (2.8, 8.8) | 4.9 (2.5, 11.8) | 0.6 |
| Fulvestrant Treatment Censoring Status | <0.001 | |||
| Â Â Â Â censoring | 5 (6.0%) | 2 (3.8%) | 3 (10%) | |
| Â Â Â Â event | 42 (51%) | 15 (28%) | 27 (90%) | |
| Â Â Â Â Unknown | 36 (43%) | 36 (68%) | 0 (0%) | |
| Ovarian Function Suppression (OFS) Received | 0.5 | |||
| Â Â Â Â No | 11 (13%) | 6 (11%) | 5 (17%) | |
| Â Â Â Â Unknown | 50 (60%) | 31 (58%) | 19 (63%) | |
| Â Â Â Â Yes | 22 (27%) | 16 (30%) | 6 (20%) | |
| Ever Received AKT Inhibitor | 5 (6.0%) | 2 (3.8%) | 3 (10%) | 0.3 |
| Ever Received CDK4/6 Inhibitor | 25 (30%) | 11 (21%) | 14 (47%) | 0.026 |
| Duration of CDK4/6 Inhibitor Therapy (Months) | 4.6 (2.4, 9.6) | 8.1 (3.7, 15.7) | 3.2 (2.3, 5.0) | 0.2 |
| CDK4/6 Inhibitor Censoring Status | 0.044 | |||
| Â Â Â Â censoring | 3 (3.6%) | 1 (1.9%) | 2 (6.7%) | |
| Â Â Â Â event | 22 (27%) | 10 (19%) | 12 (40%) | |
| Â Â Â Â Unknown | 58 (70%) | 42 (79%) | 16 (53%) | |
| Ever Received mTOR Inhibitor | 30 (36%) | 9 (17%) | 21 (70%) | <0.001 |
| Duration of mTOR Inhibitor Therapy (Months) | 4.5 (1.9, 10.1) | 6.2 (2.6, 6.7) | 3.6 (1.9, 11.6) | >0.9 |
| mTOR Inhibitor Censoring Status | <0.001 | |||
| Â Â Â Â event | 29 (35%) | 9 (17%) | 20 (67%) | |
| Â Â Â Â Unknown | 54 (65%) | 44 (83%) | 10 (33%) | |
| Total Lines of Systemic Metastatic Therapies | 5.0 (2.0, 8.0) | 4.0 (2.0, 5.0) | 7.5 (4.0, 11.0) | <0.001 |
| Duration of First-Line (1L) Treatment (Months) | 6.0 (2.8, 14.9) | 4.3 (2.8, 15.4) | 8.0 (4.2, 14.9) | 0.2 |
| First-Line (1L) Treatment Censoring Status | 0.4 | |||
| Â Â Â Â censoring | 3 (3.6%) | 3 (5.7%) | 0 (0%) | |
| Â Â Â Â event | 79 (95%) | 49 (92%) | 30 (100%) | |
| Â Â Â Â Unknown | 1 (1.2%) | 1 (1.9%) | 0 (0%) | |
| Duration of Second-Line (2L) Treatment (Months) | 4.2 (2.3, 14.0) | 3.5 (1.7, 12.3) | 6.9 (3.0, 15.8) | 0.2 |
| Second-Line (2L) Treatment Censoring Status | 0.006 | |||
| Â Â Â Â censoring | 3 (3.6%) | 2 (3.8%) | 1 (3.3%) | |
| Â Â Â Â event | 68 (82%) | 39 (74%) | 29 (97%) | |
| Â Â Â Â Unknown | 12 (14%) | 12 (23%) | 0 (0%) | |
| 1 n (%); Median (Q1, Q3) | ||||
| 2 Fisher’s Exact Test for Count Data with simulated p-value (based on 2000 replicates); Wilcoxon rank sum test; Pearson’s Chi-squared test; Fisher’s exact test; Wilcoxon rank sum exact test | ||||
# --- Table 2: Univariate Cox ---
merged_tab_uni <- merge(
uni_unweighted_final[, .(Variable, HR_un = HR, L_un = Lower, U_un = Upper, P_un = P)],
uni_weighted_final[, .(Variable, HR_w = HR, L_w = Lower, U_w = Upper, P_w = P)],
by = "Variable", all = TRUE)
merged_tab_uni[, `:=`(
`Unweighted HR (95% CI)` = sprintf("%.2f (%.2f\u2013%.2f)", HR_un, L_un, U_un),
`Weighted HR (95% CI)` = sprintf("%.2f (%.2f\u2013%.2f)", HR_w, L_w, U_w))]
gt_table_2 <- merged_tab_uni[, .(Variable,
`Unweighted HR (95% CI)`, `P-value (Unweighted)` = P_un,
`Weighted HR (95% CI)`, `P-value (Weighted)` = P_w)] |>
gt() |>
tab_header(title = "Table 2. Univariate Cox Proportional Hazards Analysis for Overall Survival") |>
cols_align(align = "left", columns = Variable) |>
cols_align(align = "center", columns = -Variable) |>
tab_options(
table.border.top.color = "black", table.border.top.width = px(2),
table.border.bottom.color = "black", table.border.bottom.width = px(2),
column_labels.border.bottom.color = "black",
column_labels.border.bottom.width = px(1.5),
table_body.border.bottom.color = "black",
table_body.border.bottom.width = px(1.5),
data_row.padding = px(6))
gt_table_2| Table 2. Univariate Cox Proportional Hazards Analysis for Overall Survival | ||||
| Variable | Unweighted HR (95% CI) | P-value (Unweighted) | Weighted HR (95% CI) | P-value (Weighted) |
|---|---|---|---|---|
| AKT Wildtype | 4.45 (1.06–18.67) | 0.041 | 3.66 (0.61–21.98) | 0.156 |
| Age > 55 at Diagnosis | 0.90 (0.42–1.93) | 0.782 | 0.72 (0.32–1.62) | 0.430 |
| Intensive Regimen (Group A) | 0.35 (0.15–0.81) | 0.013 | 0.38 (0.15–0.92) | 0.032 |
| Liver Involved | 2.23 (0.85–5.84) | 0.103 | 3.39 (1.41–8.14) | 0.006 |
| Other Sites (Non-Liver) | 1.04 (0.46–2.34) | 0.933 | 1.30 (0.54–3.16) | 0.563 |
| Tumor Grade: 3 or Unknown | 1.52 (0.75–3.09) | 0.243 | 1.94 (0.92–4.09) | 0.083 |
# --- Table 3: Multivariate Cox ---
mv_merged <- merge(
mv_unweighted_final[, .(Variable, HR_un = HR, L_un = Lower, U_un = Upper, P_un = P)],
mv_weighted_final[, .(Variable, HR_w = HR, L_w = Lower, U_w = Upper, P_w = P)],
by = "Variable", all = TRUE)
mv_merged[, `:=`(
`Unweighted HR (95% CI)` = sprintf("%.2f (%.2f\u2013%.2f)", HR_un, L_un, U_un),
`Weighted HR (95% CI)` = sprintf("%.2f (%.2f\u2013%.2f)", HR_w, L_w, U_w))]
gt_table_3 <- mv_merged[, .(Variable,
`Unweighted HR (95% CI)`, `P-value (Unweighted)` = P_un,
`Weighted HR (95% CI)`, `P-value (Weighted)` = P_w)] |>
gt() |>
tab_header(
title = "Table 3. Multivariate Cox Proportional Hazards Analysis for Overall Survival",
subtitle = "Mutually Adjusted Models (Unweighted vs. Weighted)") |>
cols_align(align = "left", columns = Variable) |>
cols_align(align = "center", columns = -Variable) |>
tab_options(
table.border.top.color = "black", table.border.top.width = px(2),
table.border.bottom.color = "black", table.border.bottom.width = px(2),
column_labels.border.bottom.color = "black",
column_labels.border.bottom.width = px(1.5),
table_body.border.bottom.color = "black",
table_body.border.bottom.width = px(1.5),
data_row.padding = px(6))
gt_table_3| Table 3. Multivariate Cox Proportional Hazards Analysis for Overall Survival | ||||
| Mutually Adjusted Models (Unweighted vs. Weighted) | ||||
| Variable | Unweighted HR (95% CI) | P-value (Unweighted) | Weighted HR (95% CI) | P-value (Weighted) |
|---|---|---|---|---|
| AKT Wildtype | 3.68 (0.84–16.13) | 0.084 | 3.64 (0.54–24.41) | 0.184 |
| Age > 55 at Diagnosis | 0.83 (0.38–1.84) | 0.647 | 0.79 (0.33–1.92) | 0.610 |
| Intensive Regimen (Group A) | 0.32 (0.13–0.80) | 0.016 | 0.32 (0.14–0.74) | 0.008 |
| Liver Involved | 2.12 (0.78–5.73) | 0.140 | 3.42 (1.43–8.18) | 0.006 |
| Other Sites (Non-Liver) | 1.18 (0.51–2.72) | 0.703 | 1.47 (0.59–3.67) | 0.406 |
| Tumor Grade: 3 or Unknown | 1.95 (0.90–4.22) | 0.091 | 1.87 (0.78–4.45) | 0.160 |
# --- Table 4: Primary Outcomes ---
get_median_string_by_name <- function(fit, group_name = NULL) {
sum_fit <- summary(fit)
if (is.null(group_name)) {
med <- sum_fit$table["median"]
lower <- sum_fit$table["0.95LCL"]
upper <- sum_fit$table["0.95UCL"]
} else {
row_names <- rownames(sum_fit$table)
row_index <- grep(group_name, row_names)
med <- sum_fit$table[row_index, "median"]
lower <- sum_fit$table[row_index, "0.95LCL"]
upper <- sum_fit$table[row_index, "0.95UCL"]
}
paste0(
if (is.na(med) | is.infinite(med)) "NE" else round(med, 1), " (",
if (is.na(lower) | is.infinite(lower)) "NE" else round(lower, 1), "\u2013",
if (is.na(upper) | is.infinite(upper)) "NE" else round(upper, 1), ")")
}
get_cox_metrics <- function(model) {
sum_mod <- summary(model)
hr <- round(sum_mod$conf.int[1, "exp(coef)"], 2)
lower <- round(sum_mod$conf.int[1, "lower .95"], 2)
upper <- round(sum_mod$conf.int[1, "upper .95"], 2)
p_val <- sum_mod$coefficients[1, ncol(sum_mod$coefficients)]
list(
hr_str = paste0(hr, " (",
if (is.na(lower) | is.infinite(lower)) "NE" else lower, "\u2013",
if (is.na(upper) | is.infinite(upper)) "NE" else upper, ")"),
p_str = ifelse(p_val < 0.001, "<0.001", sprintf("%.3f", p_val)))
}
cox_unweighted_m <- get_cox_metrics(cox_model_unweighted)
cox_weighted_m <- get_cox_metrics(cox_model_weighted)
stable_data <- data.table(
Section = c("Total Study Population",
rep("Unweighted Primary Analysis", 2),
rep("Entropy-Balanced Weighted Analysis", 2)),
Group = c("Overall Cohort",
"Group B (Reference)", "Group A",
"Group B (Reference)", "Group A"),
Median_OS = c(
get_median_string_by_name(km_fit_overall),
get_median_string_by_name(km_fit, "Group B"),
get_median_string_by_name(km_fit, "Group A"),
get_median_string_by_name(km_fit_weighted, "Group B"),
get_median_string_by_name(km_fit_weighted, "Group A")),
HR = c("NE", "1.00 (Reference)", cox_unweighted_m$hr_str,
"1.00 (Reference)", cox_weighted_m$hr_str),
P_Val = c("NE", "NE", cox_unweighted_m$p_str, "NE", cox_weighted_m$p_str)
)
gt_table_4 <- stable_data |>
gt(groupname_col = "Section") |>
tab_header(title = "Table 4. Clinical Outcomes and Hazard Ratios for Overall Survival Across Analysis Cohorts") |>
cols_label(Group = "Analysis Cohort & Group",
Median_OS = "Median OS, Months (95% CI)",
HR = "Hazard Ratio (95% CI)", P_Val = "P-value") |>
cols_align(align = "left", columns = Group) |>
cols_align(align = "center", columns = c(Median_OS, HR, P_Val)) |>
tab_options(
table.border.top.color = "black", table.border.top.width = px(2),
table.border.bottom.color = "black", table.border.bottom.width = px(2),
column_labels.border.bottom.color = "black",
column_labels.border.bottom.width = px(1.5),
table_body.border.bottom.color = "black",
table_body.border.bottom.width = px(1.5),
row_group.font.weight = "bold", row_group.background.color = "gray95",
data_row.padding = px(6), table.font.size = px(12))
gt_table_4| Table 4. Clinical Outcomes and Hazard Ratios for Overall Survival Across Analysis Cohorts | |||
| Analysis Cohort & Group | Median OS, Months (95% CI) | Hazard Ratio (95% CI) | P-value |
|---|---|---|---|
| Total Study Population | |||
| Overall Cohort | 77.7 (55.9–NE) | NE | NE |
| Unweighted Primary Analysis | |||
| Group B (Reference) | 55.9 (46.4–NE) | 1.00 (Reference) | NE |
| Group A | 121.6 (77.7–NE) | 0.35 (0.15–0.81) | 0.013 |
| Entropy-Balanced Weighted Analysis | |||
| Group B (Reference) | 55.9 (46.4–NE) | 1.00 (Reference) | NE |
| Group A | 89.9 (52.3–NE) | 0.38 (0.15–0.92) | 0.032 |
# --- Tables 5 & 6: AMEs ---
format_ame_tbl <- function(raw_dt, title_text) {
dt <- copy(raw_dt)
dt[, `:=`(ame_pct = AME * 100,
lower_pct = lower * 100,
upper_pct = upper * 100)]
dt[, `:=`(
`AME (95% CI)` = sprintf("%+.1f%% (%+.1f%% to %+.1f%%)",
ame_pct, lower_pct, upper_pct),
`P-value` = ifelse(p < 0.001, "<0.001", sprintf("%.3f", p)))]
disp <- dt[, .(Factor = factor, `AME (95% CI)`, `P-value`)]
setorder(disp, -`AME (95% CI)`)
disp |> gt() |>
tab_header(title = title_text) |>
cols_align(align = "left", columns = Factor) |>
cols_align(align = "center", columns = -Factor) |>
tab_options(
table.border.top.color = "black", table.border.top.width = px(2),
table.border.bottom.color = "black", table.border.bottom.width = px(2),
column_labels.border.bottom.color = "black",
column_labels.border.bottom.width = px(1.5),
table_body.border.bottom.color = "black",
table_body.border.bottom.width = px(1.5),
data_row.padding = px(7))
}
gt_table_5 <- format_ame_tbl(ame_table_b, "Table 5. Average Marginal Effects for Group B Selection")
gt_table_6 <- format_ame_tbl(ame_table_chemo, "Table 6. Average Marginal Effects for 1L Chemotherapy Selection")
gt_table_5| Table 5. Average Marginal Effects for Group B Selection | ||
| Factor | AME (95% CI) | P-value |
|---|---|---|
| Other Sites (Non-Liver) | -9.8% (-32.6% to +13.1%) | 0.401 |
| Tumor Grade: 3 or Unknown | -23.2% (-42.9% to -3.6%) | 0.020 |
| Liver Involved | +25.2% (+3.8% to +46.5%) | 0.021 |
| AKT Wildtype | +20.5% (-3.8% to +44.9%) | 0.099 |
| Age > 55 at Diagnosis | +17.9% (-0.7% to +36.4%) | 0.059 |
| Table 6. Average Marginal Effects for 1L Chemotherapy Selection | ||
| Factor | AME (95% CI) | P-value |
|---|---|---|
| Age > 55 at Diagnosis | -8.4% (-26.1% to +9.3%) | 0.352 |
| AKT Wildtype | -11.1% (-32.9% to +10.7%) | 0.318 |
| Liver Involved | +43.6% (+16.9% to +70.3%) | 0.001 |
| Other Sites (Non-Liver) | +30.2% (+10.5% to +49.9%) | 0.003 |
| Tumor Grade: 3 or Unknown | +3.1% (-14.6% to +20.7%) | 0.735 |
gtsave(pretty_demographics, "Table1_Patient_Demographics.docx")
gtsave(gt_table_2, "Table2_Univariate_OS.docx")
gtsave(gt_table_3, "Table3_Multivariate_OS.docx")
gtsave(gt_table_4, "Table4_Primary_Outcomes.docx")
gtsave(gt_table_5, "Table5_AME_GroupB.docx")
gtsave(gt_table_6, "Table6_AME_1L_Chemo.docx")
cat("All tables exported successfully.\n")## All tables exported successfully.12. Publication Figures
This section assembles and exports the two primary composite figures: a four-panel survival figure and a two-panel AME forest plot figure.
# --- Four-panel Survival Figure ---
extract_and_combine <- function(ggsurv_obj, label_char) {
plot_km <- ggsurv_obj$plot +
labs(title = paste0("(", label_char, ") ", ggsurv_obj$plot$labels$title))
plot_rt <- ggsurv_obj$table +
labs(y = "Cohort", fill = "Cohort", color = "Cohort") +
theme(plot.title = element_blank(), axis.title.x = element_blank())
(plot_km / plot_rt) + plot_layout(ncol = 1, heights = c(3, 1))
}
combined_survival_figure <- wrap_plots(
list(extract_and_combine(p1, "a"),
extract_and_combine(p2, "b")),
ncol = 2)
print(combined_survival_figure)
ggsave("combined_survival_plots.png", combined_survival_figure,
width = 19, height = 8, dpi = 300)
# --- Two-panel AME Figure (shared x-axis) ---
pooled_bounds <- rbind(ame_table_b[, .(lower, upper)],
ame_table_chemo[, .(lower, upper)])
ax_master <- ame_xaxis(pooled_bounds)
plot_b_shared <- build_ame_plot(
ame_table_b, ax_master,
"Clinical Drivers of Group B (Control) Selection",
"Change in Probability of Receiving Group B")
plot_chemo_shared <- build_ame_plot(
ame_table_chemo, ax_master,
"Clinical Drivers of 1L Chemotherapy Selection",
"Change in Probability of Receiving 1L Chemotherapy")
combined_ame_figure <- (plot_b_shared + plot_chemo_shared) +
plot_layout(ncol = 2, guides = "collect") &
theme(axis.title.x = element_text(size = 10, face = "bold")) &
plot_annotation(tag_levels = "a")
print(combined_ame_figure)
Supplementary Analysis: Weight Truncation Sensitivity
To evaluate the influence of extreme entropy balancing weights, we performed a sensitivity analysis by truncating weights at a fixed cap of 4. This threshold was selected to limit the influence of highly weighted observations while preserving the majority of the weighted sample.
# ---- 1. Apply weight truncation ----
cap_val <- 4
eligible[, ws_trim := pmin(ws, cap_val)]
cat("Weight cap applied:", cap_val, "\n")## Weight cap applied: 4##
## --- Weight Distribution ---## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4801 1.3257 1.5995 2.0000 2.2915 7.9747## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4801 1.3257 1.5995 1.8886 2.2915 4.0000# ---- 3. Effective Sample Size (ESS) ----
calc_ess <- function(w) (sum(w)^2) / sum(w^2)
ess_original <- calc_ess(eligible$ws)
ess_trimmed <- calc_ess(eligible$ws_trim)
cat("\n--- Effective Sample Size ---\n")##
## --- Effective Sample Size ---## Original ESS: 58.7## Trimmed ESS: 67.2# ---- 4. Refit Cox model with truncated weights ----
cox_model_trimmed <- coxph(
Surv(survival_time, survival_status) ~ therapy_flag,
data = eligible,
weights = ws_trim,
robust = TRUE
)
cat("\n--- Trimmed-weight Cox model ---\n")##
## --- Trimmed-weight Cox model ---## Call:
## coxph(formula = Surv(survival_time, survival_status) ~ therapy_flag,
## data = eligible, weights = ws_trim, robust = TRUE)
##
## n= 83, number of events= 31
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## therapy_flagGroup A -1.0384 0.3540 0.2933 0.4482 -2.317 0.0205 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## therapy_flagGroup A 0.354 2.825 0.1471 0.8522
##
## Concordance= 0.648 (se = 0.052 )
## Likelihood ratio test= 13.19 on 1 df, p=3e-04
## Wald test = 5.37 on 1 df, p=0.02
## Score (logrank) test = 13.6 on 1 df, p=2e-04, Robust = 6.33 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).# ---- 5. Extract results for comparison table ----
extract_hr <- function(model) {
s <- summary(model)
data.table(
HR = round(s$conf.int[1, "exp(coef)"], 3),
Lower = round(s$conf.int[1, "lower .95"], 3),
Upper = round(s$conf.int[1, "upper .95"], 3),
P = format_p(s$coefficients[1, "Pr(>|z|)"])
)
}
res_original <- extract_hr(cox_model_weighted)
res_trimmed <- extract_hr(cox_model_trimmed)
comparison_table <- rbindlist(list(
cbind(Model = "Primary (No Truncation)", res_original),
cbind(Model = "Truncated (Cap = 4)", res_trimmed)
))
datatable(
comparison_table,
rownames = FALSE,
options = list(dom = 't'),
caption = htmltools::tags$caption(
style = 'caption-side: top; text-align: left; font-weight: bold;',
'Supplement Table: Sensitivity Analysis with Weight Truncation'
)
)