Vẽ đường Kaplan-Meier bằng ggplot R là một trong những việc mình làm nhiều nhất khi phân tích dữ liệu theo dõi dọc. Gói survminer kết hợp với ggsurvplot cho ra biểu đồ đẹp hơn RevMan, tùy chỉnh được màu sắc và nhãn, và thêm bảng rủi ro ngay dưới đường cong. Template mình chia sẻ ở đây là file mình dùng trong thực tế, paste vào là chạy, chỉ cần đổi tên biến.
Trước khi đọc bài này, chắc chắn bạn đã có dữ liệu sống còn với ít nhất 2 cột: thời gian theo dõi và biến cố. Nếu cần ôn lại cách làm sạch và chuẩn bị dữ liệu, xem quy trình R cho bài báo lâm sàng.
Cài gói và chuẩn bị dữ liệu
install.packages(c("survival", "survminer"))
library(survival)
library(survminer)
# Giả sử df có cột:
# - thoi_gian: số tháng theo dõi
# - bien_co: 1 = xảy ra, 0 = kiểm duyệt
# - nhom: "Phau thuat A" hoặc "Phau thuat B"
# Tạo đối tượng Surv
surv_obj <- Surv(time = df$thoi_gian, event = df$bien_co)
# Fit Kaplan-Meier theo nhóm
km_fit <- survfit(surv_obj ~ nhom, data = df)
Vẽ đường cong cơ bản
ggsurvplot(
fit = km_fit,
data = df,
pval = TRUE, # giá trị p log-rank
conf.int = TRUE, # khoảng tin cậy 95%
risk.table = TRUE, # bảng số bệnh nhân còn lại theo thời gian
xlab = "Thoi gian theo doi (thang)",
ylab = "Ti le khong bien co",
legend.title = "Nhom",
legend.labs = c("Phau thuat A", "Phau thuat B"),
palette = c("#E41A1C", "#377EB8")
)
Đây là phiên bản cơ bản. Biểu đồ ra có đường cong, khoảng tin cậy mờ xung quanh, giá trị p log-rank ở góc trên trái, và bảng rủi ro ngay dưới.
Chỉnh định dạng cho đăng báo
Biểu đồ mặc định chưa đẹp bằng yêu cầu journal. Mình thêm một vài tùy chọn:
ggsurvplot(
fit = km_fit,
data = df,
pval = TRUE,
pval.method = TRUE, # hiện "Log-rank" trên biểu đồ
conf.int = TRUE,
risk.table = TRUE,
risk.table.col = "strata",
xlab = "Thoi gian theo doi (thang)",
ylab = "Ti le song con tich luy",
legend.title = "Nhom dieu tri",
legend.labs = c("Phau thuat A", "Phau thuat B"),
palette = c("#E41A1C", "#377EB8"),
ggtheme = theme_bw(), # nền trắng, trông chuyên nghiệp hơn
font.main = c(14, "plain"),
font.x = c(12, "plain"),
font.y = c(12, "plain"),
font.tickslab = c(10, "plain"),
tables.height = 0.3, # tỷ lệ chiều cao bảng rủi ro so với biểu đồ
tables.theme = theme_cleantable()
)
Trong phân tích sống còn một nghiên cứu thuần tập của mình, mình dùng template gần như y chang, chỉ đổi tên nhóm và chỉnh cột thời gian từ ngày sang tháng. Reviewer không yêu cầu chỉnh thêm, biểu đồ pass ngay lần đầu.
Xuất file ảnh chất lượng cao
duong_cong <- ggsurvplot(
fit = km_fit,
data = df,
pval = TRUE,
conf.int = TRUE,
risk.table = TRUE,
palette = c("#E41A1C", "#377EB8"),
ggtheme = theme_bw()
)
# Xuất PNG 300 dpi cho nộp bài
ggsave(
"figure_km_curve.png",
print(duong_cong),
dpi = 300,
width = 8,
height = 6,
units = "in"
)
Journal thường yêu cầu 300 dpi minimum. Code trên xuất PNG đúng chuẩn đó. Nếu cần TIFF, đổi tên file thành .tiff và thêm device = "tiff".
Kiểm tra giả định log-rank
Trước khi báo giá trị p log-rank trong bài, kiểm tra giả định hazards tỷ lệ (proportional hazards):
# Kiểm tra proportional hazards bằng Schoenfeld residuals
kiem_tra_ph <- cox.zph(coxph(surv_obj ~ nhom, data = df))
print(kiem_tra_ph)
ggcoxzph(kiem_tra_ph)
Nếu p nhỏ (thường dưới 0.05 theo Schoenfeld test), giả định bị vi phạm, không dùng log-rank được nữa mà cần phân tích stratified hoặc weighted Kaplan-Meier. Reviewer từ tạp chí lớn thường hỏi điều này.
Median survival và khoảng tin cậy
# Xem median và 95% CI cho từng nhóm
summary(km_fit)$table
# Hoặc trực tiếp
km_fit
Đây là số liệu mình đưa vào bảng kết quả bên cạnh biểu đồ. Reviewer thường muốn thấy median survival cụ thể, không chỉ nhìn vào đường cong.
Để xem thêm về cách tổ chức quy trình phân tích R cho toàn bài báo, từ giai đoạn làm sạch đến xuất bảng và biểu đồ, đọc thêm ở lộ trình tự học R cho bác sĩ.
Nếu bạn muốn học phân tích sống còn cùng toàn bộ quy trình R cho paper lâm sàng một cách có hệ thống, khóa R cho nghiên cứu y khoa bao gồm module riêng cho survival analysis.