Kaplan–Meier 生存曲线

# Kaplan–Meier 生存曲线：两组生存曲线及风险人数表
# 项目原创脚本，MIT 许可；模拟数据 CC0-1.0；图形 CC BY 4.0。
# 在全新 R 会话运行。每次运行会更新下方 output_dir 中的本示例文件。
# 图形标签使用英文以避免不同电脑缺少中文字体；可在参数和 labs() 中修改。

# 1. 依赖：首次使用时手动运行以下安装命令，脚本不会自动安装。
# install.packages(c("ggplot2", "jsonlite", "survival", "patchwork"))
required_packages <- c("ggplot2", "jsonlite", "survival", "patchwork")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages)) stop("Please install: ", paste(missing_packages, collapse = ", "))
library(ggplot2)

# 2. 可修改参数：颜色、字体、导出尺寸和路径。输出路径相对于当前工作目录。
set.seed(20261014)
colors <- c("#007F96", "#B66A46")
font_family <- "sans"
figure_width <- 8
figure_height <- 7
output_dir <- file.path(getwd(), "rfig-survival-output")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
theme_set(theme_minimal(base_size = 12, base_family = font_family) +
 theme(plot.background = element_rect(fill = "white", color = NA), panel.grid.minor = element_blank(),
       plot.title = element_text(face = "bold", size = 15), legend.position = "bottom",
       plot.margin = margin(14, 18, 14, 14)))

# 3. 示例数据：每行一位受试者，记录随访时间与事件状态。 替换 time 和 event，确认 1 为事件、0 为删失，时间单位一致，每人一行。
# 使用自己的 CSV 时，将以下数据生成段替换为：
# data <- read.csv("your-data.csv", stringsAsFactors = FALSE)
# 再按本例要求设置分类顺序、检查字段、单位及缺失值。
data <- data.frame(id = 1:160, group = factor(rep(c("Control", "Treatment"), each = 80), levels = c("Control", "Treatment")))
event_time <- rexp(160, rate = ifelse(data$group == "Control", .055, .035))
censor_time <- runif(160, 6, 30)
data$time <- pmin(event_time, censor_time); data$event <- as.integer(event_time <= censor_time)
stopifnot(!anyNA(data))

# 4. 统计计算与绘图。使用 Kaplan–Meier 乘积极限估计，十字标记为删失。风险人数为各指定时间点之前仍处于风险中的人数。未展示置信带或组间检验。演示采用独立右删失；真实数据需明确终点和删失机制。
fit <- survival::survfit(survival::Surv(time, event) ~ group, data = data)
curve <- data.frame(time = fit$time, survival = fit$surv, censor = fit$n.censor, group = rep(sub("group=", "", names(fit$strata)), fit$strata))
curve <- rbind(data.frame(time = 0, survival = 1, censor = 0, group = levels(data$group)), curve)
curve$group <- factor(curve$group, levels = levels(data$group)); curve <- curve[order(curve$group, curve$time), ]
times <- c(0, 6, 12, 18, 24, 30)
risk <- summary(fit, times = times, extend = TRUE)
risk_data <- data.frame(time = risk$time, n = risk$n.risk, group = sub("group=", "", as.character(risk$strata)))
curve_plot <- ggplot(curve, aes(time, survival, color = group)) + geom_step(aes(linetype = group), linewidth = .9) +
 geom_point(data = subset(curve, censor > 0), shape = 3, size = 2) + scale_color_manual(values = colors) +
 scale_x_continuous(breaks = times, limits = c(0, 30)) + scale_y_continuous(limits = c(0,1), labels = function(x) paste0(x*100, "%")) +
 labs(x = "Follow-up (months)", y = "Survival probability", color = NULL, linetype = NULL, title = "Kaplan-Meier survival estimates")
risk_plot <- ggplot(risk_data, aes(time, group, label = n, color = group)) + geom_text(size = 3.5) +
 scale_color_manual(values = colors) + scale_x_continuous(breaks = times, limits = c(0,30)) + labs(x = "Follow-up (months)", y = NULL, title = "Number at risk") +
 theme(legend.position = "none", panel.grid = element_blank(), plot.title = element_text(size = 11))
p <- patchwork::wrap_plots(curve_plot, risk_plot, heights = c(3.5, 1.2))
at12 <- summary(fit, times = 12)
metrics <- list(survival12 = at12$surv, risk12 = at12$n.risk)
print(p)

# 5. 导出：CSV 与绘图所用 data 完全一致，preview 与科研图片均来自同一个 p。
write.csv(data, file.path(output_dir, "data.csv"), row.names = FALSE, fileEncoding = "UTF-8")
ggsave(file.path(output_dir, "preview.png"), p, width = figure_width, height = figure_height, dpi = 160, bg = "white")
ggsave(file.path(output_dir, "figure.png"), p, width = figure_width, height = figure_height, dpi = 300, bg = "white")
ggsave(file.path(output_dir, "figure.pdf"), p, width = figure_width, height = figure_height, device = grDevices::pdf, bg = "white")
ggsave(file.path(output_dir, "figure.tiff"), p, width = figure_width, height = figure_height, dpi = 300, compression = "lzw", bg = "white")
jsonlite::write_json(metrics, file.path(output_dir, "metrics.json"), auto_unbox = TRUE, digits = 16, pretty = TRUE)