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

這一段的主要目的是「把髒資料洗乾淨」。

# =====================================================================
# 1. 日期修復與讀取資料
# =====================================================================
# 為什麼要寫這個 function?因為癌登資料庫的日期常有缺漏(例如只有 201805),
# 我們把它統一補成 15 號(20180515),這樣 R 才能計算存活時間。
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)))

# 讀取癌登大表
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))

# =====================================================================
# 2. 處理重複的病歷號 (ID) -> 確保一人一筆資料
# =====================================================================
# 針對重複 ID 寫邏輯!如果是雙側就拆開,如果是同側同一年就取最新的一筆紀錄。」
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 討論。

# =====================================================================
# 大池子:抓出所有 HER2 陽性的病人
# =====================================================================
step1 <- final_data %>% filter(ihc %in% c("HR+/HER2+", "HR-/HER2+"))

# =====================================================================
# 主分析世代建立 (Step 7)
# =====================================================================
step7 <- step1 %>%  
  mutate(NAC_binary = case_when(NAC %in% c(10,11,20,21,30,31,40,41,990,999) ~ 1, TRUE ~ 0))

step7_her2 <- step7 %>%
  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)

cat("👉 [主線] 全體 HER2+ 分析人數:", nrow(step7_her2), "人\n")
## 👉 [主線] 全體 HER2+ 分析人數: 1332 人

三、存活時間定義與防呆修正

這一段是統計上的致命關鍵。處理了五年截尾,並且修復了 Recur_Type == 70 的問題。

# =====================================================================
# 1. 存活定義與修正 (先幫所有病人算好存活時間)
# =====================================================================
step7_her2_clean <- step7_her2 %>%
  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 的病人
    rfs_event_raw = case_when(
      Recur_Type == 70 ~ NA_integer_,                          
      Recur_Type %in% c(1:69, 88) & !is.na(Recur_Date2) ~ 1L,  
      TRUE ~ 0L                                                
    ),
    rfs_date_used = if_else(rfs_event_raw == 1L, Recur_Date2, FU_DT),
    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 年截尾
    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)
  )

# 存出主線資料庫
write.csv(step7_her2_clean, "final_her2_data.csv", row.names = FALSE)

# =====================================================================
# 2. 路徑 B:從「已經算好存活時間」的乾淨母體中,抽出極早期小腫瘤!
# =====================================================================
step6_subgroup <- step7_her2_clean %>% 
  filter(`X_1st_Surg` != 0) %>%                                                     # 條件1:有開刀
  filter(NAC == 988) %>%                                                            # 條件2:未做 NAC
  filter(X3.10.Pathologic.T %in% c("1", "1A", "1B", "1C", "1MIC", "1Mic", "1mic")) %>% # 條件3:腫瘤極小 (T1)
  filter(X3.11.Pathologic.N == "0") %>%                                             # 條件4:淋巴結無轉移 (N0)
  filter(X3.12.Pathologic.M != "1")                                                 # 條件5:無遠端轉移 (M0)

# 現在存出來的 csv 裡面就會有 surv_year_rfs_5y 這個欄位了!
write.csv(step6_subgroup, "subgroup_t1n0m0_data.csv", row.names = FALSE)
cat("[支線] 成功篩選出 pT1N0M0 小腫瘤世代人數:", nrow(step6_subgroup), "人\n")
## [支線] 成功篩選出 pT1N0M0 小腫瘤世代人數: 365 人

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

建立四個治療組別的比較表。

# =====================================================================
# 製表與變數轉換
# =====================================================================
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")),
    
    # 防呆修正:沒寫分化的病人改為 "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
  )

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

tbl1
Characteristics Group A
N = 240
1
Group B
N = 30
1
Group C
N = 70
1
Group D
N = 992
1
Overall
N = 1,332
1
p-value2
Age




<0.001
    Mean ± SD; Median [Min, Max] 58 ± 12; 58 [32, 89] 64 ± 16; 68 [32, 86] 55 ± 12; 57 [30, 79] 54 ± 11; 54 [20, 87] 55 ± 12; 55 [20, 89]
age_group




0.008
    < 50 y/o 57 (24%) 7 (23%) 22 (31%) 344 (35%) 430 (32%)
    ≥ 50 y/o 183 (76%) 23 (77%) 48 (69%) 648 (65%) 902 (68%)
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




<0.001
    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.070
    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%)
follow_up_months




0.10
    Mean ± SD; Median [Min, Max] 66 ± 41; 60 [0, 155] 65 ± 39; 58 [9, 147] 77 ± 42; 78 [5, 156] 75 ± 39; 71 [0, 163] 73 ± 40; 67 [0, 163]
final_stage




<0.001
    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 Pearson’s Chi-squared test
tbl1 %>% as_flex_table() %>% autofit() %>% save_as_docx(path = "Table1_Clinicopathologic_Characteristics.docx")

五、存活分析與 Cox 比例風險模型

展示各期別存活曲線,並以迴歸校正干擾因子。

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

# 1. 繪製全期別的 KM 存活圖
for (stg in stages) {
  df_stage <- df_analysis %>% filter(final_stage == stg)
  if (nrow(df_stage) > 0 && n_distinct(df_stage$group) >= 2) {
    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), title = paste0("5-Year RFS - Stage ", stg), risk.table.height = 0.25, tables.theme = theme_cleantable())
    
    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), title = paste0("5-Year OS - Stage ", stg), risk.table.height = 0.25, tables.theme = theme_cleantable())
    
    grid.newpage()
    grid.arrange(p_rfs$plot, p_os$plot, p_rfs$table, p_os$table, ncol = 2, heights = c(4, 1.2))
  }
}

# 2. 跑 Cox 模型 
library(survival)
library(dplyr)
library(broom)
library(flextable)
library(officer)

# === 定義生存物件 (改用我們整理好的、最乾淨完整的 df_analysis) ===
# OS (使用完整的追蹤年份與事件)
surv_os <- Surv(time = df_analysis$surv_year_os, event = df_analysis$Vita_m == 1)

# RFS (使用完整的追蹤年份與事件)
surv_rfs <- Surv(time = df_analysis$surv_year_rfs, event = df_analysis$rfs_event_raw == 1)

# === 要跑的變數 ===
vars <- c("group", "Age", "hr_status", "grade_f", "final_stage")

# --- function: 單因子 + 多因子 Cox ---
run_cox <- function(surv_obj, data, vars) {
  
  # 單因子 Cox
  uni_results <- lapply(vars, function(v) {
    f <- as.formula(paste("surv_obj ~", v))
    coxph(f, data = data)
  })
  
  uni_table <- lapply(uni_results, tidy, 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)
  
  # 多因子 Cox
  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)
  
  # 合併單因子與多因子
  cox_table <- 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 = c("term")
    )
  
  # ✨ 智慧攔截:處理 Inf 或因為 0 事件導致的荒謬超大數值 (>100) → 轉為 NA
  cox_table <- cox_table %>%
    mutate(across(c(Uni_HR, Uni_LCL, Uni_UCL, Multi_HR, Multi_LCL, Multi_UCL),
                  ~ ifelse(. > 100 | is.infinite(.), NA, .)))
  
  return(cox_table)
}

# === 格式化成表格(乾淨 Level + 顯示 ref group) ===
format_cox <- function(cox_table, title) {
  
  cox_table_fmt <- cox_table %>%
    mutate(
      # 乾淨的 Level 名稱
      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
      ),
      # 數值格式化 (遇到 NA 就顯示 Not estimable)
      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)
  
  # 各變數的 reference group
  ref_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)",      "",   "",     "",     ""
  )
  
  # 插入 ref 並讓它排在各變數最前面
  cox_table_fmt <- bind_rows(cox_table_fmt, ref_rows) %>%
    mutate(Level = trimws(Level)) %>%
    arrange(factor(Variable, levels = c("group","Age","hr_status","grade_f","final_stage")),
            case_when(
              grepl("ref", Level) ~ 0,  # ref 永遠排最上面
              TRUE ~ 1
            ),
            Level)
  
  flextable(cox_table_fmt) %>%
    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)
}

# === 跑 OS & RFS ===
cox_os <- run_cox(surv_os, df_analysis, vars)
cox_rfs <- run_cox(surv_rfs, df_analysis, vars)

# === 轉換成表格 (直接在 RMarkdown 顯示) ===
tbl_os <- format_cox(cox_os, "Table 3a. Cox Regression for Overall Survival (OS)")
tbl_rfs <- format_cox(cox_rfs, "Table 3b. Cox Regression for Recurrence-Free Survival (RFS)")

tbl_os

Table 3a. Cox Regression for Overall Survival (OS)

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

tbl_rfs

Table 3b. Cox Regression for Recurrence-Free Survival (RFS)

Variable

Level

Univariate HR (95% CI)

Univariate p

Multivariate HR (95% CI)

Multivariate p

group

Group A (ref)

group

Group B

2.07 (0.84–5.09)

0.112

0.68 (0.26–1.82)

0.448

group

Group C

1.08 (0.50–2.34)

0.837

0.53 (0.23–1.22)

0.136

group

Group D

0.66 (0.41–1.05)

0.080

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.363

hr_status

HR negative (ref)

hr_status

HR positive

0.72 (0.49–1.04)

0.078

0.76 (0.52–1.12)

0.168

grade_f

Grade I (ref)

grade_f

II

Not estimable

0.993

Not estimable

0.995

grade_f

III

Not estimable

0.993

Not estimable

0.995

grade_f

Unknown

Not estimable

0.993

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.39 (1.32–4.32)

0.004

3.93 (2.03–7.60)

0.000

final_stage

2B

1.90 (0.92–3.91)

0.082

3.87 (1.72–8.70)

0.001

final_stage

3A

2.08 (0.95–4.55)

0.065

4.41 (1.84–10.55)

0.001

final_stage

3B

7.40 (2.73–20.07)

0.000

10.42 (3.63–29.89)

0.000

final_stage

3C

7.59 (4.26–13.51)

0.000

15.02 (7.40–30.48)

0.000

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

這就是你原本想留著的趨勢圖!用來看醫師給藥習慣的演變。

# =====================================================================
# NAC 趨勢折線圖
# =====================================================================
# 將 size 改為新版的 linewidth,解決警告問題
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) +
  facet_wrap(~ final_stage, ncol = 3) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(title = "NAC Adoption Rate Over Time by Stage and HR Status", x = "Year of Diagnosis", y = "NAC Rate", color = "HR Status") +
  theme_minimal(base_size = 14) + theme(legend.position = "bottom")

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

Step 2 到 Step 6,就是在為這一步鋪陳!針對這些早期患者,真的有必要打那麼多藥嗎?我們特別把這群人獨立出來看他們的存活率。

# =====================================================================
# 極早期患者降階治療探討
# =====================================================================
# 讀取第二節 (Step 6) 存出來的極早期病人資料 (366人)
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) {
  # 我們特別幫這些小腫瘤的病人重新定義存活計算,畫出專屬他們的圖
  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), title = "Discussion: 5-Year RFS for pT1N0M0 (No NAC)", risk.table.height = 0.25)
  print(p_sub_rfs)
}