Trastuzumab-Based Therapy Improves Survival in HER2-Positive Breast Cancer: A Real-World Cohort Study from Taipei Veterans General Hospital

一、資料前處理與病歷整理 (Data Curation)

本區塊之學術目的在於解決台灣癌症登記資料庫（癌登）中常見的日期缺漏與病患重複登錄問題，確保世代研究中「一人僅能擁有一筆最精準之原發癌紀錄」。

# 1. 模糊日期標準化函數 (月/日缺失者，統一補正為 15 號)
fix_date_string <- function(x) {
  x <- str_replace_all(as.character(x), "[^0-9]", "")
  x[nchar(x) != 8] <- NA
  x_valid <- !is.na(x)
  fixed <- x
  if (any(x_valid)) {
    month_part <- as.integer(substr(x[x_valid], 5, 6))
    day_part   <- as.integer(substr(x[x_valid], 7, 8))
    month_fixed <- ifelse(month_part >= 1 & month_part <= 12, sprintf("%02d", month_part), "06")
    day_fixed <- ifelse(day_part >= 1 & day_part <= 31, sprintf("%02d", day_part), "15")
    fixed[x_valid] <- paste0(substr(x[x_valid], 1, 4), month_fixed, day_fixed)
  }
  return(fixed)
}
safe_ymd <- function(x) suppressWarnings(ymd(fix_date_string(x)))

# 2. 讀取北榮癌登大表
data0 <- read.csv("C:/Users/user/OneDrive/文件/公衛實習/癌登/concise2_1140803.csv", fileEncoding = "UTF-8") %>%
  mutate(ID = as.character(ID), DIAG_DT = safe_ymd(X2.5.Date.of.Initial.Diagnosis), FU_DT = safe_ymd(FU_Date))

# 3. 處理重複病歷號 (確保一人一筆資料)
id_freq <- table(data0$ID)
data1 <- data0 %>% filter(ID %in% names(id_freq[id_freq == 1]))
data2 <- data0 %>% filter(ID %in% names(id_freq[id_freq == 2])) %>% arrange(ID)
data3 <- data0 %>% filter(ID %in% names(id_freq[id_freq >= 3]))

split_data <- lapply(unique(data2$ID), function(i) {
  tmp <- data2 %>% filter(ID == i) %>% mutate(year = year(DIAG_DT), LATERAL = X2.7.Laterality, CASITE = ICD, HIST = Hist, CLASS = X2.3.Class.of.Case)
  if (length(unique(tmp$LATERAL)) > 1) { tmp$ID <- paste0(tmp$ID, "_l", tmp$LATERAL); return(tmp) }
  if (tmp$year[1] == tmp$year[2] && tmp$CASITE[1] == tmp$CASITE[2] && tmp$HIST[1] == tmp$HIST[2] && abs(as.numeric(diff(tmp$DIAG_DT))) <= 30) {
    class12 <- tmp %>% filter(CLASS %in% c(1, 2))
    return(if (nrow(class12) > 0) class12[1, , drop = FALSE] else tmp %>% slice_max(FU_DT, n = 1))
  }
  tmp$ID <- paste0(tmp$ID, "_m", 1:2)
  return(tmp)
})
data2_processed <- bind_rows(split_data)
data3_processed <- data3 %>% group_by(ID) %>% slice_max(DIAG_DT, n = 1) %>% ungroup()

all_cols <- unique(c(names(data1), names(data2_processed), names(data3_processed)))
fill_missing <- function(df, cols) { for (c in setdiff(cols, names(df))) df[[c]] <- NA; df[, cols] }

# 建立最終無重複、清洗完成之大世代資料庫
final_data <- bind_rows(fill_missing(data1, all_cols), fill_missing(data2_processed, all_cols), fill_missing(data3_processed, all_cols))

二、存活變數計算與兩大世代建立

本區進行核心資料流處理：首先計算好所有病患的存活天數與五年截尾變數，接著像漏斗一樣分流出「主分析大世代」與「次族群小世代」，確保兩份資料表皆具備完整的存活變數。 研究邏輯！從「所有 HER2+病人」這個大池子（Step1），分出兩條路徑： * 路徑 A（主線 Step 7）：抓取所有期別的病人，用來看全體存活率與跑 Cox 模型。 * 路徑 B（支線 Step 2~6）：一步步嚴格篩選出「有開刀、沒做化療、腫瘤極小的 pT1N0M0 病人」，留到報告最後的 Discussion 討論。

# =====================================================================
# 1. 抓出 HER2 陽性母體，並「預先算好所有存活變數」
# =====================================================================
step1 <- final_data %>% filter(ihc %in% c("HR+/HER2+", "HR-/HER2+"))

step1_clean <- step1 %>%
  mutate(
    DIAG_DT = as.Date(DIAG_DT), FU_DT = as.Date(FU_DT),
    Recur_Date2 = suppressWarnings(ymd(ifelse(is.na(Recur_Date) | Recur_Date %in% c(0, 99999999), NA, as.character(Recur_Date)))),
    
    # 修正 RFS 事件：排除從未 disease-free 的病人 (設為 NA)
    rfs_event_raw = case_when(
  Recur_Type == 70 ~ NA_integer_,
  Recur_Type == 99 ~ NA_integer_,

  Recur_Type %in% c(
    4, 6, 10, 13, 14, 15, 16, 17,
    20, 21, 22, 25, 26, 27, 30, 36,
    40, 46, 51, 52, 53, 54, 55, 56,
    57, 58, 59, 60, 62, 88
  ) & !is.na(Recur_Date2) ~ 1L,

  Recur_Type == 0 ~ 0L,
  TRUE ~ NA_integer_
),

    rfs_date_used = if_else(rfs_event_raw == 1L, Recur_Date2, FU_DT),

    
    # 【完整追蹤期變數】(供 Cox 迴歸使用)
    surv_year_rfs = if_else(is.na(rfs_event_raw), NA_real_, as.numeric(pmax(rfs_date_used - DIAG_DT, 0)) / 365.25),
    surv_year_os  = as.numeric(survival) / 12,
    Vita_m        = if_else(!is.na(Vital) & Vital == 0, 1, 0),
    
    # 【5 年截尾存活變數】(供 KM Plot 視覺化使用)
    relapse_5y       = case_when(is.na(rfs_event_raw) ~ NA_real_, rfs_event_raw == 1L & surv_year_rfs <= censor_years ~ 1, TRUE ~ 0),
    surv_year_rfs_5y = if_else(is.na(surv_year_rfs), NA_real_, pmin(surv_year_rfs, censor_years)),
    Vita_m_5y        = if_else(Vita_m == 1 & surv_year_os <= censor_years, 1, 0),
    surv_year_os_5y  = pmin(surv_year_os, censor_years)
  )

# =====================================================================
# 2. 🚀 路徑 A：主分析世代建立 (全期別)
# =====================================================================
step7_her2_clean <- step1_clean %>%  
  mutate(NAC_binary = case_when(NAC %in% c(10,11,20,21,30,31,40,41,990,999) ~ 1, TRUE ~ 0)) %>%
  mutate(
    final_stage = factor(case_when(
      NAC_binary == 1 ~ X3.7.Clinical.Stage.Group,
      NAC_binary == 0 ~ X3.13.Pathologic.Stage.Group,
      TRUE ~ NA_character_
    ), levels = stages)
  ) %>% filter(!is.na(final_stage) & final_stage %in% stages)

# 寫出主分析資料庫
write.csv(step7_her2_clean, "final_her2_data.csv", row.names = FALSE)
cat("👉 [主線] 全體 HER2+ 且擁有明確期別之主要分析人數：", nrow(step7_her2_clean), "人\n")

## 👉 [主線] 全體 HER2+ 且擁有明確期別之主要分析人數： 1332 人

# =====================================================================
# 3. 🔎 路徑 B：極早期小腫瘤漏斗篩選 (Patient Flowchart)
# =====================================================================
cat("\n=== 🔎 臨床次族群篩選：Patient Flowchart ===\n")

## 
## === 🔎 臨床次族群篩選：Patient Flowchart ===

step2 <- step1_clean %>% filter(X_1st_Surg != 0)
cat("👉 Step 2 (接受第一線手術開刀者):", nrow(step2), "人\n")

## 👉 Step 2 (接受第一線手術開刀者): 1652 人

step3 <- step2 %>% filter(NAC == 988)
cat("👉 Step 3 (排除新輔助治療，確定為純開刀首選者):", nrow(step3), "人\n")

## 👉 Step 3 (排除新輔助治療，確定為純開刀首選者): 1333 人

step4 <- step3 %>% filter(X3.10.Pathologic.T %in% c("1", "1A", "1B", "1C", "1MIC", "1Mic", "1mic"))
cat("👉 Step 4 (病理診斷確定為腫瘤小於2公分之 T1 族群):", nrow(step4), "人\n")

## 👉 Step 4 (病理診斷確定為腫瘤小於2公分之 T1 族群): 483 人

step5 <- step4 %>% filter(X3.11.Pathologic.N == "0")
cat("👉 Step 5 (病理檢查確定淋巴結無轉移 N0 族群):", nrow(step5), "人\n")

## 👉 Step 5 (病理檢查確定淋巴結無轉移 N0 族群): 367 人

# 這裡修正了名稱，統一叫 step6_subgroup
step6_subgroup <- step5 %>% filter(X3.12.Pathologic.M != "1")
cat("👉 Step 6 (最終排除遠端轉移 M1，確定符合 pT1N0M0 之臨床總人數):", nrow(step6_subgroup), "人\n")

## 👉 Step 6 (最終排除遠端轉移 M1，確定符合 pT1N0M0 之臨床總人數): 366 人

# 寫出次族群資料庫 (此表已自動繼承 surv_year_rfs_5y)
write.csv(step6_subgroup, "subgroup_t1n0m0_data.csv", row.names = FALSE)

# 檢查原始資料 ID 出現次數
id_freq <- table(data0$ID)

cat("Original unique IDs:", length(id_freq), "\n")

## Original unique IDs: 9770

cat("IDs with 1 record:", sum(id_freq == 1), "\n")

## IDs with 1 record: 9226

cat("IDs with 2 records:", sum(id_freq == 2), "\n")

## IDs with 2 records: 539

cat("IDs with 3+ records:", sum(id_freq >= 3), "\n")

## IDs with 3+ records: 5

# 檢查最終 HER2+ 世代的原始 ID 是否有多筆紀錄
step7_id_check <- step7_her2_clean %>%
  mutate(original_ID = sub("_[lm][0-9]+$", "", ID)) %>%
  count(original_ID, name = "n_records") %>%
  filter(n_records > 1)

cat("Original IDs represented by multiple analysis rows:",
    nrow(step7_id_check), "\n")

## Original IDs represented by multiple analysis rows: 14

print(head(step7_id_check, 20))

##    original_ID n_records
## 1     22238894         2
## 2     22616409         2
## 3     22709069         2
## 4     26232704         2
## 5     30076161         2
## 6     36579438         2
## 7     38162526         2
## 8     38208257         2
## 9     39458064         2
## 10    42131327         2
## 11    42339064         2
## 12    43765474         2
## 13    44577769         2
## 14     4524938         2

step7_her2_clean %>%
  count(CT, TT, .drop = FALSE) %>%
  arrange(CT, TT) %>%
  print()

##   CT TT   n
## 1  0  0 240
## 2  0  1  30
## 3  1  0  70
## 4  1  1 992

step7_her2_clean %>%
  summarise(
    n_rows = n(),
    n_unique_ID = n_distinct(ID),
    n_missing_CT = sum(is.na(CT)),
    n_missing_TT = sum(is.na(TT))
  ) %>%
  print()

##   n_rows n_unique_ID n_missing_CT n_missing_TT
## 1   1332        1332            0            0

三、患者基線特徵 (Table 1)

建立四大治療組別的基線特徵比較表，並補正 Grade 缺失值以保留檢定力。

# 讀取具有乾淨存活欄位的主線資料庫
df_analysis <- read.csv("final_her2_data.csv") %>%
  mutate(
    # 定義治療組 A 到 D
    group = factor(case_when(
      CT == 0 & TT == 0 ~ "Group A", 
      CT == 0 & TT == 1 ~ "Group B",
      CT == 1 & TT == 0 ~ "Group C", 
      CT == 1 & TT == 1 ~ "Group D",
      TRUE ~ NA_character_
    ), levels = c("Group A", "Group B", "Group C", "Group D")),
    
    age_group = ifelse(Age < 50, "< 50 y/o", "≥ 50 y/o"),
    hr_status = factor(case_when(grepl("HR\\+", ihc) ~ "HR positive", grepl("HR\\-", ihc) ~ "HR negative", TRUE ~ NA_character_), levels = c("HR negative", "HR positive")),
    
    # 防呆修正：將 Grade 缺失值轉換為 "Unknown" 類別，避免被 R 自動刪除
    grade_f = factor(ifelse(is.na(Grade) | !(Grade %in% c(1, 2, 3)), "Unknown", as.character(Grade)), levels = c("1", "2", "3", "Unknown"), labels = c("I", "II", "III", "Unknown")),
    histology = factor(case_when(Hist == 8500 ~ "IDC", Hist == 8520 ~ "ILC", TRUE ~ "Others"), levels = c("IDC", "ILC", "Others")),
    endocrine = factor(ifelse(HT == 1, "Yes", "No"), levels = c("No", "Yes")),
    rt_value = suppressWarnings(as.numeric(X4.2.1.1.RT.Target.Summary)),
    rt_done  = factor(ifelse(!is.na(rt_value) & rt_value > 0, "Done", "None"), levels = c("None", "Done")),
    follow_up_months = surv_year_os * 12
  )


# ============================================================
# Table 1. Baseline Characteristics by Treatment Group
# ============================================================

tbl1 <- df_analysis %>%
  select(
    group,
    Age,
    hr_status,
    grade_f,
    histology,
    endocrine,
    rt_done,
    final_stage
  ) %>%
  tbl_summary(
    by = group,
    type = list(Age ~ "continuous2"),
    statistic = list(
      Age ~ "{mean} ± {sd}; {median} [{min}, {max}]",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = list(Age ~ 1),
    missing = "ifany"
  ) %>%
  add_overall(last = TRUE) %>%
  add_p(
    test = list(
      Age ~ "kruskal.test"
    )
  ) %>%
  modify_header(
    label = "**Characteristics**"
  ) %>%
  bold_labels()

# 顯示表格
tbl1

Characteristics	Group A N = 2401	Group B N = 301	Group C N = 701	Group D N = 9921	Overall N = 1,3321	p-value2
Age						<0.001
Mean ± SD; Median [Min, Max]	58.1 ± 11.9; 57.5 [32.0, 89.0]	63.8 ± 16.2; 67.5 [32.0, 86.0]	54.9 ± 12.1; 56.5 [30.0, 79.0]	54.3 ± 11.3; 54.0 [20.0, 87.0]	55.3 ± 11.7; 55.0 [20.0, 89.0]
hr_status						<0.001
HR negative	148 (62%)	9 (30%)	23 (33%)	428 (43%)	608 (46%)
HR positive	92 (38%)	21 (70%)	47 (67%)	564 (57%)	724 (54%)
grade_f
I	7 (2.9%)	0 (0%)	0 (0%)	4 (0.4%)	11 (0.8%)
II	76 (32%)	7 (23%)	27 (39%)	270 (27%)	380 (29%)
III	44 (18%)	9 (30%)	24 (34%)	276 (28%)	353 (27%)
Unknown	113 (47%)	14 (47%)	19 (27%)	442 (45%)	588 (44%)
histology						0.14
IDC	224 (93%)	27 (90%)	67 (96%)	946 (95%)	1,264 (95%)
ILC	2 (0.8%)	2 (6.7%)	0 (0%)	12 (1.2%)	16 (1.2%)
Others	14 (5.8%)	1 (3.3%)	3 (4.3%)	34 (3.4%)	52 (3.9%)
endocrine	89 (37%)	21 (70%)	43 (61%)	536 (54%)	689 (52%)	<0.001
rt_done						<0.001
None	168 (70%)	22 (73%)	49 (70%)	489 (49%)	728 (55%)
Done	72 (30%)	8 (27%)	21 (30%)	503 (51%)	604 (45%)
final_stage
1A	191 (80%)	12 (40%)	30 (43%)	251 (25%)	484 (36%)
1B	0 (0%)	0 (0%)	0 (0%)	7 (0.7%)	7 (0.5%)
2A	35 (15%)	5 (17%)	28 (40%)	287 (29%)	355 (27%)
2B	6 (2.5%)	6 (20%)	3 (4.3%)	175 (18%)	190 (14%)
3A	2 (0.8%)	2 (6.7%)	4 (5.7%)	130 (13%)	138 (10%)
3B	4 (1.7%)	0 (0%)	0 (0%)	17 (1.7%)	21 (1.6%)
3C	2 (0.8%)	5 (17%)	5 (7.1%)	125 (13%)	137 (10%)
1 n (%)
2 Kruskal-Wallis rank sum test; Pearson’s Chi-squared test; NA; Fisher’s exact test

# 匯出 Word
tbl1 %>%
  as_flex_table() %>%
  autofit() %>%
  save_as_docx(
    path = "Table1_Clinicopathologic_Characteristics.docx"
  )

四、5年存活率分析 (KM Plots)

使用 5 年截尾時間變數，並套用指定的經典四色色票（A:紅、B:藍、C:綠、D:紫）進行全期別分層展示。

df_analysis <- df_analysis %>% 
  filter(!is.na(group)) %>% 
  mutate(group = relevel(factor(group), ref = "Group A"), final_stage = relevel(factor(final_stage), ref = "1A"))

for (stg in stages) {
  df_stage <- df_analysis %>% filter(final_stage == stg)
  if (nrow(df_stage) > 0 && n_distinct(df_stage$group) >= 2) {
    
    # 5年無復發存活率 (RFS) Plot
    fit_rfs <- survfit(Surv(surv_year_rfs_5y, relapse_5y) ~ group, data = df_stage)
    p_rfs <- ggsurvplot(fit_rfs, data = df_stage, pval = TRUE, risk.table = TRUE, xlim = c(0, 5), palette = km_colors, title = paste0("5-Year RFS - Stage ", stg), risk.table.height = 0.25, tables.theme = theme_cleantable())
    
    # 5年整體存活率 (OS) Plot
    fit_os <- survfit(Surv(surv_year_os_5y, Vita_m_5y) ~ group, data = df_stage)
    p_os <- ggsurvplot(fit_os, data = df_stage, pval = TRUE, risk.table = TRUE, xlim = c(0, 5), palette = km_colors, title = paste0("5-Year OS - Stage ", stg), risk.table.height = 0.25, tables.theme = theme_cleantable())
    
    # 使用 grid 引擎左右合體排列
    grid.newpage()
    grid.arrange(p_rfs$plot, p_os$plot, p_rfs$table, p_os$table, ncol = 2, heights = c(4, 1.2))
  }
}

# ============================================================
# Log-rank tests by stage: 5-year OS and RFS
# Place this after the KM plotting loop
# ============================================================

get_logrank_by_stage <- function(data, time_var, event_var, outcome) {

  bind_rows(lapply(stages, function(stg) {

    dat <- data %>%
      filter(final_stage == stg) %>%
      select(all_of(c("group", time_var, event_var))) %>%
      filter(
        !is.na(group),
        !is.na(.data[[time_var]]),
        !is.na(.data[[event_var]])
      ) %>%
      mutate(group = droplevels(factor(group)))

    n_total <- nrow(dat)
    n_events <- sum(dat[[event_var]] == 1)
    n_groups <- n_distinct(dat$group)

    # Return NA if there are insufficient groups or no events
    if (n_groups < 2 || n_events == 0) {
      return(tibble(
        Outcome = outcome,
        Stage = stg,
        N = n_total,
        Events = n_events,
        Groups = n_groups,
        Chi_square = NA_real_,
        df = NA_integer_,
        p_value = NA_real_
      ))
    }

    # Overall log-rank test across treatment groups
    fit <- survdiff(
      as.formula(
        paste0("Surv(", time_var, ", ", event_var, ") ~ group")
      ),
      data = dat
    )

    tibble(
      Outcome = outcome,
      Stage = stg,
      N = n_total,
      Events = n_events,
      Groups = n_groups,
      Chi_square = unname(fit$chisq),
      df = length(fit$n) - 1L,
      p_value = pchisq(
        fit$chisq,
        df = length(fit$n) - 1L,
        lower.tail = FALSE
      )
    )
  }))
}

# Overall log-rank p-value for each stage
logrank_os <- get_logrank_by_stage(
  df_analysis,
  time_var = "surv_year_os_5y",
  event_var = "Vita_m_5y",
  outcome = "OS"
)

logrank_rfs <- get_logrank_by_stage(
  df_analysis,
  time_var = "surv_year_rfs_5y",
  event_var = "relapse_5y",
  outcome = "RFS"
)

logrank_results <- bind_rows(logrank_os, logrank_rfs) %>%
  mutate(
    p_display = case_when(
      is.na(p_value) ~ "Not estimable",
      p_value < 0.001 ~ "<0.001",
      TRUE ~ sprintf("%.3f", p_value)
    )
  ) %>%
  select(
    Outcome, Stage, N, Events, Groups,
    Chi_square, df, p_value, p_display
  )

print(logrank_results, n = Inf)

## # A tibble: 14 × 9
##    Outcome Stage     N Events Groups Chi_square    df   p_value p_display    
##    <chr>   <chr> <int>  <int>  <int>      <dbl> <int>     <dbl> <chr>        
##  1 OS      1A      484      9      4       1.90     3  5.93e- 1 0.593        
##  2 OS      1B        7      0      1      NA       NA NA        Not estimable
##  3 OS      2A      355     16      4      26.4      3  8.04e- 6 <0.001       
##  4 OS      2B      190     15      4      31.3      3  7.42e- 7 <0.001       
##  5 OS      3A      138     12      4      14.7      3  2.14e- 3 0.002        
##  6 OS      3B       21      4      2       3.45     1  6.34e- 2 0.063        
##  7 OS      3C      137     22      4       2.62     3  4.54e- 1 0.454        
##  8 RFS     1A      481     17      4       9.17     3  2.71e- 2 0.027        
##  9 RFS     1B        7      0      1      NA       NA NA        Not estimable
## 10 RFS     2A      353     29      4      31.2      3  7.57e- 7 <0.001       
## 11 RFS     2B      189     13      4       4.93     3  1.77e- 1 0.177        
## 12 RFS     3A      136     10      4      24.4      3  2.08e- 5 <0.001       
## 13 RFS     3B       21      4      2       3.88     1  4.90e- 2 0.049        
## 14 RFS     3C      135     36      4      57.7      3  1.81e-12 <0.001

五、Cox 比例風險模型 (Cox Regression)

評估各項臨床預後因子對於無復發存活（RFS）與整體存活（OS）之多變項影響（此處採用完整追蹤期變數以獲取最高之統計檢定力與長期 Hazard Ratio 估計）。

# 基於統計效力考量，採用完整追蹤天數
surv_os_full  <- Surv(time = df_analysis$surv_year_os,  event = df_analysis$Vita_m == 1)
surv_rfs_full <- Surv(time = df_analysis$surv_year_rfs, event = df_analysis$rfs_event_raw == 1)
vars <- c("group", "Age", "hr_status", "grade_f", "final_stage")

# 單多因子合併函數
run_cox <- function(surv_obj, data, vars) {
  uni_table <- lapply(vars, function(v) tidy(coxph(as.formula(paste("surv_obj ~", v)), data = data), exponentiate = TRUE, conf.int = TRUE)) %>% bind_rows(.id = "Variable") %>% mutate(Variable = vars[as.numeric(Variable)]) %>% select(Variable, term, estimate, conf.low, conf.high, p.value)
  
  multi_model <- coxph(surv_obj ~ group + Age + hr_status + grade_f + final_stage, data = data)
  multi_table <- tidy(multi_model, exponentiate = TRUE, conf.int = TRUE) %>% select(term, estimate, conf.low, conf.high, p.value)
  
  uni_table %>% rename(Uni_HR = estimate, Uni_LCL = conf.low, Uni_UCL = conf.high, Uni_p = p.value) %>% left_join(multi_table %>% rename(Multi_HR = estimate, Multi_LCL = conf.low, Multi_UCL = conf.high, Multi_p = p.value), by = "term") %>% mutate(across(c(Uni_HR, Uni_LCL, Uni_UCL, Multi_HR, Multi_LCL, Multi_UCL), ~ ifelse(. > 100 | is.infinite(.), NA, .)))
}

# 格式化輸出函數
format_cox <- function(cox_table, title) {
  cox_table %>% mutate(Level = case_when(grepl("group", term) ~ sub("group", "", term), grepl("hr_status", term) ~ sub("hr_status", "", term), grepl("grade_f", term) ~ sub("grade_f", "", term), grepl("final_stage", term) ~ sub("final_stage", "", term), term == "Age" ~ "Age (per year increase)", TRUE ~ term), Uni = ifelse(is.na(Uni_HR), "Not estimable", sprintf("%.2f (%.2f–%.2f)", Uni_HR, Uni_LCL, Uni_UCL)), Multi = ifelse(is.na(Multi_HR), "Not estimable", sprintf("%.2f (%.2f–%.2f)", Multi_HR, Multi_LCL, Multi_UCL)), Uni_p = ifelse(is.na(Uni_p), "", sprintf("%.3f", Uni_p)), Multi_p = ifelse(is.na(Multi_p), "", sprintf("%.3f", Multi_p))) %>% select(Variable, Level, Uni, Uni_p, Multi, Multi_p) %>% bind_rows(tribble(~Variable, ~Level, ~Uni, ~Uni_p, ~Multi, ~Multi_p, "group", "Group A (ref)", "", "", "", "", "hr_status", "HR negative (ref)", "", "", "", "", "grade_f", "Grade I (ref)", "", "", "", "", "final_stage", "Stage 1A (ref)", "", "", "", "")) %>% arrange(factor(Variable, levels = c("group", "Age", "hr_status", "grade_f", "final_stage")), case_when(grepl("ref", Level) ~ 0, TRUE ~ 1), Level) %>% flextable() %>% set_header_labels(Variable="Variable", Level="Level", Uni="Univariate HR (95% CI)", Uni_p="Univariate p", Multi="Multivariate HR (95% CI)", Multi_p="Multivariate p") %>% add_header_row(values = title, colwidths = 6) %>% autofit() %>% bold(i = ~ as.numeric(Uni_p) < 0.05 | as.numeric(Multi_p) < 0.05, bold = TRUE)
}

# 印出標準發表格式表格
format_cox(run_cox(surv_os_full, df_analysis, vars), "Table 2a. Cox Regression for Overall Survival (Full Follow-up)")

Table 2a. Cox Regression for Overall Survival (Full Follow-up)
Variable	Level	Univariate HR (95% CI)	Univariate p	Multivariate HR (95% CI)	Multivariate p
group	Group A (ref)
group	Group B	2.33 (0.94–5.76)	0.068	0.71 (0.27–1.90)	0.500
group	Group C	1.29 (0.59–2.81)	0.524	0.97 (0.42–2.23)	0.935
group	Group D	0.59 (0.36–0.97)	0.039	0.35 (0.19–0.64)	0.001
Age	Age (per year increase)	1.07 (1.05–1.09)	0.000	1.05 (1.03–1.07)	0.000
hr_status	HR negative (ref)
hr_status	HR positive	0.81 (0.54–1.22)	0.320	0.80 (0.52–1.24)	0.318
grade_f	Grade I (ref)
grade_f	II	0.49 (0.12–2.02)	0.324	0.19 (0.04–0.88)	0.034
grade_f	III	0.53 (0.13–2.20)	0.385	0.12 (0.03–0.59)	0.009
grade_f	Unknown	0.10 (0.02–0.49)	0.004	0.04 (0.01–0.19)	0.000
final_stage	Stage 1A (ref)
final_stage	1B	0.00 (0.00–NA)	0.994	0.00 (0.00–NA)	0.995
final_stage	2A	2.53 (1.26–5.09)	0.009	3.25 (1.49–7.09)	0.003
final_stage	2B	3.75 (1.81–7.78)	0.000	7.53 (3.20–17.69)	0.000
final_stage	3A	3.72 (1.70–8.16)	0.001	7.51 (3.06–18.44)	0.000
final_stage	3B	9.27 (2.99–28.75)	0.000	7.70 (2.30–25.76)	0.001
final_stage	3C	7.23 (3.62–14.46)	0.000	12.28 (5.39–27.99)	0.000

format_cox(run_cox(surv_rfs_full, df_analysis, vars), "Table 2b. Cox Regression for Recurrence-Free Survival (Full Follow-up)")

Table 2b. Cox Regression for Recurrence-Free Survival (Full Follow-up)
Variable	Level	Univariate HR (95% CI)	Univariate p	Multivariate HR (95% CI)	Multivariate p
group	Group A (ref)
group	Group B	2.06 (0.84–5.05)	0.116	0.68 (0.26–1.82)	0.445
group	Group C	1.10 (0.51–2.37)	0.811	0.54 (0.24–1.24)	0.148
group	Group D	0.65 (0.41–1.04)	0.075	0.25 (0.14–0.46)	0.000
Age	Age (per year increase)	1.03 (1.01–1.04)	0.002	1.01 (0.99–1.02)	0.373
hr_status	HR negative (ref)
hr_status	HR positive	0.72 (0.50–1.04)	0.082	0.76 (0.52–1.12)	0.170
grade_f	Grade I (ref)
grade_f	II	Not estimable	0.994	Not estimable	0.995
grade_f	III	Not estimable	0.994	Not estimable	0.996
grade_f	Unknown	Not estimable	0.994	Not estimable	0.995
final_stage	Stage 1A (ref)
final_stage	1B	0.00 (0.00–NA)	0.994	0.00 (0.00–NA)	0.996
final_stage	2A	2.38 (1.32–4.31)	0.004	3.93 (2.03–7.61)	0.000
final_stage	2B	1.91 (0.93–3.93)	0.080	3.94 (1.75–8.87)	0.001
final_stage	3A	2.08 (0.95–4.54)	0.067	4.42 (1.85–10.58)	0.001
final_stage	3B	7.38 (2.72–20.01)	0.000	10.48 (3.65–30.07)	0.000
final_stage	3C	7.56 (4.25–13.47)	0.000	15.06 (7.43–30.56)	0.000

六、新輔助治療 (NAC) 時代演變趨勢

觀察不同期別與荷爾蒙狀態下，歷史採用新輔助化療（NAC）的趨勢。折線上已加入各年份實際的「總人數（N數）」，以防小樣本比例失真。

summary_nac <- df_analysis %>%
  mutate(year = year(DIAG_DT), NAC = as.integer(NAC_binary == 1)) %>%
  filter(!is.na(year), !is.na(hr_status)) %>%
  group_by(year, final_stage, hr_status) %>%
  summarise(total = n(), nac_rate = sum(NAC) / total, .groups = "drop") %>%
  filter(final_stage %in% c("1A", "2A", "2B", "3A", "3B", "3C"))

ggplot(summary_nac, aes(x = year, y = nac_rate, color = hr_status)) +
  geom_line(linewidth = 1) + 
  geom_point(size = 2) +
  # 數據點上直接標記樣本數 (N)
  geom_text(aes(label = total), vjust = -1, size = 3, show.legend = FALSE) +
  facet_wrap(~ final_stage, ncol = 3) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1), expand = expansion(mult = c(0.1, 0.2))) +
  labs(title = "NAC Adoption Rate Over Time by Stage and HR Status with Sample Counts", x = "Year of Diagnosis", y = "NAC Rate", color = "HR Status") +
  theme_minimal(base_size = 14) + theme(legend.position = "bottom")

七、討論區：pT1N0M0 極早期小腫瘤之療效分析 (Subgroup Discussion)

專門探討針對第二章經由漏斗篩選出的極早期小腫瘤患者，實施「降階治療」的療效。本圖已套用專屬四色色票（km_colors）。

# 讀取在第二章時已經算好存活變數的極早期小腫瘤檔案
df_sub <- read.csv("subgroup_t1n0m0_data.csv") %>%
  mutate(group = factor(case_when(
      CT == 0 & TT == 0 ~ "Group A (None)", 
      CT == 0 & TT == 1 ~ "Group B (TT only)",
      CT == 1 & TT == 0 ~ "Group C (CT only)", 
      CT == 1 & TT == 1 ~ "Group D (CT+TT)",
      TRUE ~ NA_character_
    ), levels = c("Group A (None)", "Group B (TT only)", "Group C (CT only)", "Group D (CT+TT)"))) %>% 
  filter(!is.na(group))

if (n_distinct(df_sub$group) >= 2) {
  # 繪製 pT1N0M0 專屬 5 年 RFS
  fit_sub_rfs <- survfit(Surv(surv_year_rfs_5y, relapse_5y) ~ group, data = df_sub)
  p_sub_rfs <- ggsurvplot(fit_sub_rfs, data = df_sub, pval = TRUE, risk.table = TRUE, xlim = c(0, 5), palette = km_colors, title = "5-Year RFS for pT1N0M0", risk.table.height = 0.25, tables.theme = theme_cleantable())
  
  # 繪製 pT1N0M0 專屬 5 年 OS
  fit_sub_os <- survfit(Surv(surv_year_os_5y, Vita_m_5y) ~ group, data = df_sub)
  p_sub_os <- ggsurvplot(fit_sub_os, data = df_sub, pval = TRUE, risk.table = TRUE, xlim = c(0, 5), palette = km_colors, title = "5-Year OS for pT1N0M0", risk.table.height = 0.25, tables.theme = theme_cleantable())
  
  # 雙圖輸出
  grid.newpage()
  grid.arrange(p_sub_rfs$plot, p_sub_os$plot, p_sub_rfs$table, p_sub_os$table, ncol = 2, heights = c(4, 1.2))
}