Publication-ready KM plot in R

pub_ready_KM-plot.R

# 安裝並載入所需的套件，使用pacman套件管理器
# 如果pacman未安裝，則先安裝pacman
if (!require("pacman")) install.packages("pacman")
# 使用pacman一次性載入多個套件：ggsurvfit(繪製生存曲線)、survival(生存分析)、ggplot2(繪圖系統)
pacman::p_load(ggsurvfit, survival, ggplot2, showtext) 
# Find first TTF font in current directory and load it as "default"
ttf_files <- list.files(pattern = "\\.ttf$", ignore.case = TRUE)
if (length(ttf_files) > 0) {
  font_add("default", ttf_files[1])
  cat("Loaded font:", ttf_files[1], "\n")
} else {
  warning("No TTF fonts found in current directory")
}
showtext_auto()

# 設定全局字體大小變量
base_font_size <- 16

# 創建生存分析的樣本數據集
# 設定隨機種子，確保結果可重現
set.seed(123)
# 設定樣本數量為200
n <- 200
# 生成隨機的生存時間，使用指數分佈，參數rate=0.01
time <- rexp(n, rate = 0.01)
# 生成隨機的狀態指標(0=審查，1=事件發生)，設定發生事件的概率為70%
status <- sample(0:1, n, replace = TRUE, prob = c(0.3, 0.7))
# 創建分組因子，並設定描述性標籤
# 隨機分配患者到兩個藥物組，各50%的概率
group <- factor(sample(1:2, n, replace = TRUE, prob = c(0.5, 0.5)),
                # 設定因子的水平為1和2
                levels = c(1, 2),
                # 為因子水平設定有意義的標籤：藥物A和藥物B
                labels = c("Expensive Drug", "Cheap Drug"))
# 將生成的數據組合成一個數據框
df_sample <- data.frame(
  # 生存時間
  time = time,
  # 狀態指標(0=審查，1=事件發生)
  status = status,
  # 治療組別
  group = group
)

# 創建生存曲線圖
p <- survfit2(Surv(time, status) ~ group, data = df_sample) |>
  # 使用ggsurvfit繪製生存曲線，設定線寬為1，不顯示圖例
  ggsurvfit(linewidth = 1, show.legend = FALSE) +
  # 添加審查標記，不顯示圖例
  add_censor_mark(show.legend = FALSE) +
  # 添加分位數線
  add_quantile(y_value = 0.5, color = "gray50", linewidth = 0.5, show.legend = FALSE) +
  # 添加生存曲線的顏色
  scale_ggsurvfit() +
  # 使用經典的主題 - 這只會影響KM曲線圖部分
  scale_color_manual(values = c('#FF3300', '#007AAB')) +
  theme_classic(base_size = base_font_size) +
  # 設定y軸的標題和標籤
  scale_y_continuous(
    # 設定y軸的標題
    name = "Survival Probability",
    # 設定y軸的標籤格式
    labels = scales::number_format(accuracy = 0.1),
    # 設定y軸的刻度
    breaks = seq(0, 1, by = 0.2),
    # 設定y軸的範圍
    limits = c(0, 1),  
    # 設定y軸的填充
    expand = c(0, 0)   
  ) +
  # 設定x軸的填充
  scale_x_continuous(
    # 設定x軸的填充
    name = "Months",
    expand = c(0.04, 0),
    # Convert days to months (30.44 days per month)
    breaks = seq(0, max(df_sample$time), by = 30.44 * 3),
    labels = function(x) round(x / 30.44)
  ) +
  # 設定KM曲線圖的主題 - 這些設定只會影響KM曲線圖部分，不會影響風險表
  theme(
    # 設定KM曲線圖軸線的顏色
    axis.text = element_text(color = "black"),
    # 設定KM曲線圖軸線的樣式
    axis.line = element_line(color = "black"),
    # 設定KM曲線圖刻度線的長度
    axis.ticks.length = unit(0.3, "cm"),
    # 設定KM曲線圖刻度線的樣式
    axis.ticks = element_line(color = "black", linewidth = 0.8),
    # 隱藏KM曲線圖的圖例
    legend.position = "none",  
    # 設定KM曲線圖標題的樣式
    plot.title = element_text(face = "bold", size = base_font_size + 2, hjust = 0, margin = margin(0, 0, 30, -180)),
    # 設定KM曲線圖y軸標題的樣式
    axis.title.y = element_text(vjust = 8, face = "bold", size = base_font_size),
    axis.title.x = element_text(vjust = -5, face = "bold", size = base_font_size),
  )

# 添加風險表 - 這是一個單獨的組件，有自己的主題設定
p <- p + add_risktable(
  # 顯示風險表中的風險人數
  risktable_stats = "n.risk",
  # 設定風險表中的標籤
  stats_label = list(n.risk = "No. at Risk"),
  # 設定風險表中的字體大小
  size = 5,  
  # 設定風險表的高度
  risktable_height = 0.16,
  # 設定風險表的主題 - 這些設定只會影響風險表，不會影響KM曲線圖
  theme = list(
    # 使用預設的風險表主題
    theme_risktable_default(
      # 設定標籤的字體大小
      axis.text.y.size = base_font_size,  
      # 設定標題的字體大小
      plot.title.size = base_font_size    
    ),
    # 添加自定義的風險表主題 - 這些設定只會影響風險表，不會影響KM曲線圖
    theme(
      # 設定風險表標題的樣式
      plot.title = element_text(face = "bold", hjust = 0, margin = margin(0, 0, 0, -150)),  
      # 設定風險表標籤的樣式
      axis.text.y = element_text(hjust = 0, margin = margin(r = 30))  
    )
  )
)

# 添加整體的主題設定 - 這會影響整個圖表，包括KM曲線圖和風險表
p <- p + theme(
  # 設定整個圖表的邊距 - 上、右、下、左
  plot.margin = margin(10, 20, 40, 30)
)

# 顯示基本的生存曲線圖
print(p)

# 創建自定義的生存曲線圖
p_custom <- p +
  # 限制圖表的時間範圍為0-200單位
  coord_cartesian(
    # 設定x軸的範圍 (in days, but will display as months)
    xlim = c(0, 200), 
    # 設定y軸的範圍
    ylim = c(0, 1),
    # 不使用預設的座標系統
    default = FALSE,
    # 啟用填充
    expand = TRUE,
    # 啟用裁剪
    clip = "on"
  ) +
  # Update x-axis for custom plot to show months
  scale_x_continuous(
    name = "Months",
    breaks = seq(0, 200, by = 30.44 * 1),  # Every 2 months within the 200-day range
    labels = function(x) round(x / 30.44),
    expand = c(0.04, 0)
  ) +
  # 更新圖表的標題和標籤
  labs(
    # 設定圖表的標題
    title = "A. Survival by Treatment Group"
  )

# 顯示自定義的生存曲線圖
print(p_custom)

# 儲存圖表
ggsave("survival_plot_basic.png", plot = p, width = 10, height = 6)
ggsave("survival_plot_custom.png", plot = p_custom, width = 10, height = 6)