Kaplan-Meier trong Survival Analysis bằng R
Phần 2
Liên hệ:
Facebook: R stats
Email: research@raisinghopevn.com
Điện thoại: 0965276046
Xem lại phần 1 tại đây nhé!
Thư viện cần thiết
Nếu chưa cài package pacman thì sử dụng lệnh install.packages(‘pacman’)
pacman::p_load(survival, ggplot2, survminer, ggthemes)Trong đó:
survival: phân tích survival analysis (Kaplan-Meier, Cox model,…)
ggplot2, survminer: vẽ biểu đồ đường cong survival curve
ggthemes: Các theme khác nhau cho biểu đồ đẹp hơn
Data
Số liệu veteran về ung thư phổi có sẵn trong package survival
df = survival::veteran
head(df)## trt celltype time status karno diagtime age prior
## 1 1 squamous 72 1 60 7 69 0
## 2 1 squamous 411 1 70 5 64 10
## 3 1 squamous 228 1 60 3 38 0
## 4 1 squamous 126 1 60 9 63 10
## 5 1 squamous 118 1 70 11 65 10
## 6 1 squamous 10 1 20 5 49 0Với:
celltype: 1=vảy, 2=tb nhỏ, 3=tuyến, 4=lớn
time: survival time
status: censoring status
karno: Karnofsky performance score (100=good)
diagtime: thời gian (tháng) từ lúc chẩn đoán đến lúc thử nghiệm
age: tuổi (năm)
prior: phương pháp điểu trị trước đó, 0=kkhông, 10=có
Tạo đối tượng survival object bằng Surv()
surv_object = with(df, Surv(time, status))
surv_object## [1] 72 411 228 126 118 10 82 110 314 100+ 42 8 144 25+ 11
## [16] 30 384 4 54 13 123+ 97+ 153 59 117 16 151 22 56 21
## [31] 18 139 20 31 52 287 18 51 122 27 54 7 63 392 10
## [46] 8 92 35 117 132 12 162 3 95 177 162 216 553 278 12
## [61] 260 200 156 182+ 143 105 103 250 100 999 112 87+ 231+ 242 991
## [76] 111 1 587 389 33 25 357 467 201 1 30 44 283 15 25
## [91] 103+ 21 13 87 2 20 7 24 99 8 99 61 25 95 80
## [106] 51 29 24 18 83+ 31 51 90 52 73 8 36 48 7 140
## [121] 186 84 19 45 80 52 164 19 53 15 43 340 133 111 231
## [136] 378 49Lúc này những đối tượng censored sẽ có thêm dấu +
Bắt đầu phân tích
Hàm: survfit(formula, data)
Formula của ta lúc này là Surv(time, status)~1
1 có nghĩa là tất cả, không chia nhóm.
surv_fit = survfit(Surv(time, status)~1, data = df)
summary(surv_fit, time = c(1, 300, 600, 900))## Call: survfit(formula = Surv(time, status) ~ 1, data = df)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 137 2 0.985 0.0102 0.96552 1.0000
## 300 13 113 0.117 0.0295 0.07147 0.1917
## 600 2 11 0.018 0.0126 0.00459 0.0707
## 900 2 0 0.018 0.0126 0.00459 0.07071, 300, 600, 900 là những khoảng thời gian bạn muốn báo cáo.
Ví dụ: Sau 300 ngày theo dõi còn lại 13 đối tượng. Tỷ lệ sống là 11.7% (95%CI 7,1%-19,2%)
Biểu đồ đường cong survival curve
ggsurvplot(surv_fit)
Tính trung vị
Trung vị trong survival analysis
Trung vị trong survival analysis là thời gian mà 100% số đối tượng tham gia nghiên cứu ban đầu còn lại 50%.
Bạn có thể liên tưởng định nghĩa này với khái niệm thời gian bán hủy của thuốc (thời gian mà nồng độ của thuốc trong cơ thể còn lại một nửa).
surv_fit## Call: survfit(formula = Surv(time, status) ~ 1, data = df)
##
## n events median 0.95LCL 0.95UCL
## [1,] 137 128 80 52 105Ta có trung vị là 80 ngày. Kết quả này có nghĩa là sau 80 ngày, chỉ còn lại 1/2 số bệnh nhân còn sống.
Vẽ trung vị lên biểu đồ
ggsurvplot(surv_fit, surv.median.line = "hv",
ggtheme = theme_grey())
Trong đó:
So sánh 2 nhóm điều trị với nhau (biến trt)
Lúc này ta thay số 1 trong formula thành trt là được
surv_fit = survfit(Surv(time, status)~trt, data = df)
ggsurvplot(surv_fit, pval = T, pval.method = T,
risk.table = T,
ggtheme = theme_grey())
Mô hình Cox (Cox proportional hazards model)
Ta tiến hành đưa thêm biến vào formula bằng dấu +
cox = coxph(Surv(time, status) ~ trt + celltype +
karno + diagtime +
age + prior ,
data = df)
summary(cox)## Call:
## coxph(formula = Surv(time, status) ~ trt + celltype + karno +
## diagtime + age + prior, data = df)
##
## n= 137, number of events= 128
##
## coef exp(coef) se(coef) z Pr(>|z|)
## trt 2.946e-01 1.343e+00 2.075e-01 1.419 0.15577
## celltypesmallcell 8.616e-01 2.367e+00 2.753e-01 3.130 0.00175 **
## celltypeadeno 1.196e+00 3.307e+00 3.009e-01 3.975 7.05e-05 ***
## celltypelarge 4.013e-01 1.494e+00 2.827e-01 1.420 0.15574
## karno -3.282e-02 9.677e-01 5.508e-03 -5.958 2.55e-09 ***
## diagtime 8.132e-05 1.000e+00 9.136e-03 0.009 0.99290
## age -8.706e-03 9.913e-01 9.300e-03 -0.936 0.34920
## prior 7.159e-03 1.007e+00 2.323e-02 0.308 0.75794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## trt 1.3426 0.7448 0.8939 2.0166
## celltypesmallcell 2.3669 0.4225 1.3799 4.0597
## celltypeadeno 3.3071 0.3024 1.8336 5.9647
## celltypelarge 1.4938 0.6695 0.8583 2.5996
## karno 0.9677 1.0334 0.9573 0.9782
## diagtime 1.0001 0.9999 0.9823 1.0182
## age 0.9913 1.0087 0.9734 1.0096
## prior 1.0072 0.9929 0.9624 1.0541
##
## Concordance= 0.736 (se = 0.021 )
## Likelihood ratio test= 62.1 on 8 df, p=2e-10
## Wald test = 62.37 on 8 df, p=2e-10
## Score (logrank) test = 66.74 on 8 df, p=2e-11Nhận xét:
Ung thư phổi loại tế bào nhỏ (p<0.01) và tuyến (p<0.001) làm tăng nguy cơ tử vong (HR lần lượt 2.37, 3.31).
KPS cao thì giảm nhẹ nguy cơ tử vong (HR=0.97, p <0.001).
4 dòng cuối cùng là chỉ số/kiểm định độ phù hợp của mô hình.
Hết phần 2!