KAPTwithluv
Khánh An
2023-07-28
THỰC HÀNH PHÂN TÍCH SỐNG CÒN VỚI PACKAGE SURVIVAL
Cài đặt và sử dụng package survival
Để sử dụng package survival trên
RStudio, ta thực hiện các bước sau:
Bước 1: Cài đặt package survival
Bạn có thể cài đặt package survival bằng cách chạy
lệnh sau trong RStudio:
install.packages("survival")Nếu package đã được cài đặt trước đó, bạn có thể bỏ qua bước này.
Bước 2: Tải package survival
Sau khi cài đặt, bạn cần tải package survival bằng
lệnh sau:
library(survival) Bước 3: Sử dụng các hàm của package
survival
Sau khi đã tải package survival, bạn có thể sử dụng
các hàm của package này để phân tích dữ liệu thời gian đến sự kiện. Các
hàm này bao gồm:
survfit(): Tính toán các ước lượng sống sót bằng phương pháp Kaplan-Meier hoặc phương pháp Nelson-Aalen.coxph(): Phân tích mô hình hồi quy Cox để tìm các yếu tố ảnh hưởng đến thời gian đến sự kiện.survreg(): Phân tích mô hình hồi quy Cox để tìm các yếu tố ảnh hưởng đến thang đo thời gian đến sự kiện không phải là dạng số nguyên (ví dụ: thời gian sống sót được đo bằng giờ, ngày, tháng,…).cox.zph(): Kiểm định giả thuyết về sự đồng nhất của các hệ số hồi quy Cox theo thời gian.
Dataset lung - dữ liệu về ung thư phổi
Trong suốt phần này, chúng tôi sẽ sử dụng dataset
lung từ package survival làm dữ liệu phân
tích.
Dataset lung là một trong những tập dữ liệu tiêu
biểu của package survival được sử dụng trong nhiều ví dụ và
hướng dẫn về phân tích dữ liệu. Dataset này được tạo ra từ nghiên
cứu lâm sàng về ung thư phổi và bao gồm thông tin về 228
bệnh nhân ung thư phổi với các biến như tuổi, giới tính, mức độ
suy giảm chức năng phổi,…
Dataset lung được sử dụng để minh họa cho nhiều khía cạnh khác nhau của phân tích dữ liệu, bao gồm mô hình hồi quy tuyến tính, mô hình hồi quy logistic, phân tích phân nhóm và đặc biệt là phân tích sống sót.
Để truy cập dataset lung từ package survival trong R,
ta sử dụng lệnh sau:
# gọi dataset về dữ liệu bệnh nhân của cách bệnh ung thư -> chọn dataset lung
data(cancer)Dữ liệu ung thư phổi - lung có cấu trúc như sau:
str(lung)## 'data.frame': 228 obs. of 10 variables:
## $ inst : num 3 3 3 5 1 12 7 11 1 7 ...
## $ time : num 306 455 1010 210 883 ...
## $ status : num 2 2 1 2 2 1 2 2 2 2 ...
## $ age : num 74 68 56 57 60 74 68 71 53 61 ...
## $ sex : num 1 1 1 1 1 1 2 2 1 1 ...
## $ ph.ecog : num 1 0 0 1 0 1 2 2 1 2 ...
## $ ph.karno : num 90 90 90 90 100 50 70 60 70 70 ...
## $ pat.karno: num 100 90 90 60 90 80 60 80 80 70 ...
## $ meal.cal : num 1175 1225 NA 1150 NA ...
## $ wt.loss : num NA 15 15 11 0 0 10 1 16 34 ...Mô tả các biến của dữ liệu ung thư phổi - lung
inst: Mã số của bệnh nhân.
time: Thời gian sống sót hoặc thời gian theo dõi (ngày).
status: Tình trạng sống sau thời gian theo dõi (1 = Sống sót, 2 = Qua đời).
age: Tuổi của bệnh nhân.
sex: Giới tính của bệnh nhân (1 = Nam, 2 = Nữ).
ph.ecog: Mức độ suy giảm phổi của bệnh nhân (0 = Tốt, 1 = Tàn phế nhẹ, 2 = Tàn phế trung bình, 3 = Tàn phế nặng).
ph.karno: Điểm Karnofsky của bệnh nhân, đo lường sức khỏe tổng thể (100 điểm là sức khỏe tốt nhất, 0 điểm là bệnh. nhân đã qua đời)
pat.karno: Điểm Karnofsky của bệnh nhân dựa trên người chăm sóc (nếu có).
meal.cal: Lượng calo tiêu thụ hàng ngày của bệnh nhân.
wt.loss: Mức độ giảm cân của bệnh nhân trong 6 tháng trước đó.
Các biến time và status được sử dụng trong phân tích sống còn, trong khi các biến còn lại được sử dụng để đánh giá các yếu tố ảnh hưởng đến sống sót của bệnh nhân. Các biến ph.ecog và ph.karno cung cấp thông tin về sức khỏe tổng thể của bệnh nhân, trong khi các biến meal.cal và wt.loss cung cấp thông tin về chế độ ăn uống và mức độ giảm cân của bệnh nhân.
Xử lý dữ liệu ung thư phổi - lung
Lưu ý rằng trạng thái (status) được mã hóa theo cách không chuẩn trong bộ dữ liệu này. Thông thường, dữ liệu của biến status gồm 2 giá trị: 0= còn sống, 1= đã chết.
Để thay đổi giá trị dữ liệu của biến về giá trị 0
và 1, chúng ta sẽ dùng hàm mutate() trong
package dplyr. Hàm mutate() trong package
dplyr được sử dụng để thay đổi giá trị của các biến trong
một tập dữ liệu. Để thay đổi giá trị của một biến bằng
mutate(), bạn cần xác định tên biến và cung cấp cho nó một
biểu thức hoặc hàm để tính toán giá trị mới. Quá trình thay đổi giá trị
dữ liệu của biến được thực hiện như sau:
lung <- lung %>%
mutate(status = recode(status, '1' = 0, '2' = 1))Ngoài ra, chúng ta sẽ tuổi của bệnh nhân thành 3 nhóm cho tiện quá trình phân tích:
lung <- lung %>%
mutate(age_group = case_when(
age < 50 ~ "dưới 50 tuổi",
age >= 50 & age <= 65 ~ "50 đến 65 tuổi",
age > 65 ~ "trên 65 tuổi",
TRUE ~ "không rõ"
))Bây giờ chúng ta một bộ dữ liệu mới với:
time: Thời gian sống sót hoặc thời gian theo dõi (ngày)
status: Tình trạng sống sau thời gian theo dõi (0 = Sống sót, 1 = Đã chết)
age_group: Tuổi của bệnh nhân theo nhóm tuổi (dưới 50, 50 đến 65, trên 65)
head(lung[, c("time", "status", "sex", "age_group")])## time status sex age_group
## 1 306 1 1 trên 65 tuổi
## 2 455 1 1 trên 65 tuổi
## 3 1010 0 1 50 đến 65 tuổi
## 4 210 1 1 50 đến 65 tuổi
## 5 883 1 1 50 đến 65 tuổi
## 6 1022 0 1 trên 65 tuổi Lưu ý: hàm Surv() trong package survival
mặc định chấp nhận TRUE/FALSE, trong đó TRUE là sự kiện và FALSE là kiểm
duyệt; 1/0 trong đó 1 là sự kiện và 0 là kiểm duyệt; hoặc 2/1
trong đó 2 là sự kiện và 1 là kiểm duyệt. Hãy cẩn thận để đảm bảo dữ
liệu được định dạng đúng.
Tạo đối tượng sinh tồn
Trong phân tích sống còn, để thực hiện các phương pháp phân tích dữ liệu sống sót, ta cần tạo đối tượng sống còn (survival object) dựa trên các thông tin về thời gian sống sót và trạng thái sự kiện của các đối tượng trong một nghiên cứu. Sau đó, ta có thể sử dụng các hàm để vẽ đường sống sót và ước lượng các mô hình sống sót.
Trong R, ta có thể tạo đối tượng sống còn bằng hàm
Surv(). Cú pháp của hàm này như sau:
Surv(time, event, type = "right")Trong đó:
time: Là vector chứa thông tin về thời gian sống sót hoặc thời gian theo dõi của các đối tượng. Đây có thể là thời gian từ lần khám đầu tiên đến lần khám cuối cùng, thời gian từ điều trị đến sự kiện xảy ra, hoặc thời gian từ lần phát hiện căn bệnh đến sự kiện xảy ra.
event: Là vector chứa thông tin về trạng thái sự kiện của các đối tượng. Giá trị của event có thể là 0 (còn sống) hoặc 1 (chết hoặc xảy ra sự kiện quan tâm).
type: Là tham số xác định cách tính toán thời gian sống sót. Giá trị mặc định của type là “right”, có nghĩa là thời gian sống sót được tính từ thời điểm khám cuối cùng trừ đi thời điểm xảy ra sự kiện. Nếu giá trị của type là “left”, thời gian sống sót sẽ được tính từ thời điểm khám đầu tiên trừ đi thời điểm xảy ra sự kiện.
Ví dụ, ta có dữ liệu về thời gian sống sót và trạng thái sự kiện của 10 bệnh nhân đầu tiên như sau:
Surv(lung$time, lung$status)[1:10]## [1] 306 455 1010+ 210 883 1022+ 310 361 218 166Ta thấy đối tượng 1 có sự kiện tại thời điểm 306 ngày, đối tượng 2 có sự kiện tại thời điểm 455 ngày, đối tượng 3 bị kiểm duyệt tại thời điểm 1010 ngày, …
Tính toán đường cong sinh tồn
Hàm survfit() là một hàm dùng để tính toán và vẽ
đường sống sót ước tính (survival curves) dựa trên dữ
liệu về thời gian sống sót và trạng thái sự kiện của các đối tượng trong
một nghiên cứu. Hàm survfit() là một trong những hàm cơ bản
của package survival được dùng để phân tích sống còn trên
R.
Cú pháp của hàm survfit() như sau:
survfit(formula, data, weights, ...)Trong đó:
formula: là công thức mô tả mối quan hệ giữa đối tượng sống còn và các biến giải thích. Công thức này có dạng
Surv(time, event) ~ x1 + x2 + …data: là data frame chứa các biến trong công thức.
weights: là vector chứa các trọng số cho mỗi quan sát trong mô hình.
Tôi sẽ dùng lại đối tượng sinh tồn
Surv(lung$time, lung$status) để tạo các đường cong sinh tồn
bằng phương pháp Kaplan-Meier. Tôi tạo đường cong sinh
tồn tổng thể với tham số “~1” chỉ định rằng không có
biến độc lập nào được sử dụng để phân nhóm, gán nó cho đối tượng
survfit_lung và xem xét cấu trúc bằng hàm
str():
survfit_lung <- survfit(Surv(time, status) ~ 1, data = lung)
str(survfit_lung)## List of 16
## $ n : int 228
## $ time : num [1:186] 5 11 12 13 15 26 30 31 53 54 ...
## $ n.risk : num [1:186] 228 227 224 223 221 220 219 218 217 215 ...
## $ n.event : num [1:186] 1 3 1 2 1 1 1 1 2 1 ...
## $ n.censor : num [1:186] 0 0 0 0 0 0 0 0 0 0 ...
## $ surv : num [1:186] 0.996 0.982 0.978 0.969 0.965 ...
## $ std.err : num [1:186] 0.0044 0.00885 0.00992 0.01179 0.01263 ...
## $ cumhaz : num [1:186] 0.00439 0.0176 0.02207 0.03103 0.03556 ...
## $ std.chaz : num [1:186] 0.00439 0.0088 0.00987 0.01173 0.01257 ...
## $ type : chr "right"
## $ logse : logi TRUE
## $ conf.int : num 0.95
## $ conf.type: chr "log"
## $ lower : num [1:186] 0.987 0.966 0.959 0.947 0.941 ...
## $ upper : num [1:186] 1 1 0.997 0.992 0.989 ...
## $ call : language survfit(formula = Surv(time, status) ~ 1, data = lung)
## - attr(*, "class")= chr "survfit"Một số thành phần chính của đối tượng survfit này sẽ được sử dụng để tạo các đường cong sinh tồn bao gồm:
time: các mốc thời gian mà tại đó đường cong có một bước, tức là ít nhất một sự kiện đã xảy ra.
surv: ước tính tỷ lệ sống sót tại thời điểm tương ứng.
Tuy nhiên, để so sánh tỉ lệ sống sót giữa hai giới tính nam và nữ,
tôi sẽ thêm yếu tố giới tính “~sex” vào hàm
survfit() phân nhóm thành các đường cong sống sót theo giới
tính, gán nó cho đối tượng survfit_lung_bysex và xem
cấu trúc bằng hàm str():
survfit_lung_bysex <- survfit(Surv(time, status) ~ sex, data = lung)
str(survfit_lung_bysex)## List of 17
## $ n : int [1:2] 138 90
## $ time : num [1:206] 11 12 13 15 26 30 31 53 54 59 ...
## $ n.risk : num [1:206] 138 135 134 132 131 130 129 128 126 125 ...
## $ n.event : num [1:206] 3 1 2 1 1 1 1 2 1 1 ...
## $ n.censor : num [1:206] 0 0 0 0 0 0 0 0 0 0 ...
## $ surv : num [1:206] 0.978 0.971 0.957 0.949 0.942 ...
## $ std.err : num [1:206] 0.0127 0.0147 0.0181 0.0197 0.0211 ...
## $ cumhaz : num [1:206] 0.0217 0.0291 0.0441 0.0516 0.0593 ...
## $ std.chaz : num [1:206] 0.0126 0.0146 0.018 0.0195 0.021 ...
## $ strata : Named int [1:2] 119 87
## ..- attr(*, "names")= chr [1:2] "sex=1" "sex=2"
## $ type : chr "right"
## $ logse : logi TRUE
## $ conf.int : num 0.95
## $ conf.type: chr "log"
## $ lower : num [1:206] 0.954 0.943 0.923 0.913 0.904 ...
## $ upper : num [1:206] 1 0.999 0.991 0.987 0.982 ...
## $ call : language survfit(formula = Surv(time, status) ~ sex, data = lung)
## - attr(*, "class")= chr "survfit"Ngoài ra, tôi sẽ đưa vào 2 yếu tố khác vào hàm survfit()
để phân nhóm thành các đường cong sống sót theo suy giảm chức năng phổi
và nhóm tuổi:
survfit_lung_byph <- survfit(Surv(time, status) ~ ph.ecog, data = lung)
survfit_lung_byage <- survfit(Surv(time, status) ~ age_group, data = lung)Đường cong sinh tồn Kaplan-Meier
Để vẽ đường cong sinh tồn Kaplan-Meier đơn giản trên R, tôi dùng hàm cơ bản là hàm plot() để vẽ đồ thị như sau:
plot(surv_fit, main = "Đường cong sinh tồn Kaplan Meier", xlab = "Thời gian", ylab = "Xác suất sống sót")Trong đó:
surv_fit: là biến được gán bởi hàmsurvfit()dùng để tính toán đường cong sinh tồn (xem mục 4.4)
Đường cong sinh tồn Kaplan Meier đơn giản
Tiếp tục bộ dữ liệu ung thư phổi - lung. chúng ta có đồ thị về đường cong sinh tồn Kaplan-Meier đơn giản như sau:
plot(survfit_lung, xlab = "Thời gian (ngày)", ylab = "Xác suất sống sót", main = "ĐƯờng cong sinh tồn Kaplan Meier", col = c("blue"), lwd = 2)
Ngoài hàm plot() cơ bản trong R, ta có thể sử dụng
package survminer và sử dụng hàm ggsurvplot()
để vẽ đồ thị:
ggsurvplot(survfit_lung, data = lung, risk.table = TRUE, pval = TRUE, legend.title = "Giới Tính", xlab = "Thời gian (ngày)", ylab = "Xác suất sống sót")
Đường cong sinh tồn Kaplan Meier theo yếu tố giới tính
Nâng cao hơn, ta có thể vẽ đường cong sinh tồn Kaplan-Meier cho từng nhóm được phân loại theo yếu tố. Trong ví dụ này, đường cong sinh tồn Kaplan-Meier được vẽ mỗi nhóm giới tính khác nhau trong dữ liệu ung thư phổi:
plot(survfit_lung_bysex, xlab = "Thời gian (ngày)", ylab = "Xác suất sống sót", main = "Đường cong sinh tồn Kaplan Meier theo giới tính", col = c("blue", "red"), lwd = 2)
legend("bottomleft", legend = c("Nam", "Nữ"), col = c("blue", "red"), lwd = 2, bty = "n")
Ngoài hàm plot() cơ bản trong R, ta có thể sử dụng package survminer và sử dụng hàm ggsurvplot() để vẽ đồ thị:
ggsurvplot(survfit_lung_bysex, data = lung, risk.table = TRUE, pval = TRUE, legend.title = "Giới Tính", xlab = "Thời gian (ngày)", ylab = "Xác suất sóng sót")
Đường cong sinh tồn Kaplan Meier theo yếu tố mức độ suy giảm chức năng phổi
Đường cong sinh tồn Kaplan-Meier được vẽ mỗi nhóm suy giảm chức năng phổi khác nhau trong dữ liệu ung thư phổi:
plot((survfit(Surv(time, status) ~ ph.ecog, data = lung)), xlab = "Thời gian (ngày)", ylab = "Xác suất sống sót", main = "Đường cong sinh tồn Kaplan Meier theo mức độ suy giảm chức năng phổi", col = c("green", "blue", "yellow", "red"), lwd = 2)
legend("bottomleft", legend = c("0", "1", "2", "3"), col = c("green", "blue", "yellow", "red"), lwd = 2, bty = "n")
Đường cong sinh tồn Kaplan Meier theo từng nhóm tuổi
Đường cong sinh tồn Kaplan-Meier được vẽ mỗi nhóm nhóm tuổi khác nhau trong dữ liệu ung thư phổi:
plot((survfit(Surv(time, status) ~ age_group, data = lung)), xlab = "Thời gian (ngày)", ylab = "Xác suất sống sót", main = "Đường cong sinh tồn Kaplan Meier theo nhóm tuổi", col = c("green", "blue", "yellow"), lwd = 2)
legend("bottomleft", legend = c("dưới 50 tuổi", "50 đến 65 tuổi", "trên 65 tuổi"), col = c("green", "blue", "yellow"), lwd = 2, bty = "n")
So sánh các đường cong sinh tồn
Để so sánh các đường cong sinh tồn trên R bằng
log-rank, bạn có thể sử dụng hàm
survdiff() trong gói survival. Hàm này cho
phép bạn tính toán giá trị thống kê log-rank và
giá trị p để so sánh đường cong sinh tồn giữa hai hoặc
nhiều nhóm. Sau đây là một ví dụ cụ thể về cách thực hiện so sánh các
đường cong sinh tồn bằng log-rank trên R:
survdiff(surv_obj ~ x1, data)Trong đó:
surv_obj: là đối tượng sống sót được tạo bởi hàm surv() (xem mục 4.4)x1: là yếu tố được đưa vào để phân nhóm và so sánh
Kết quả của kiểm định log-rank bao gồm giá trị chi-squared và giá trị p. Nếu giá trị p nhỏ hơn mức ý nghĩa được chọn trước (thường là 0,05), bạn có thể kết luận rằng có sự khác biệt đáng kể về tỷ lệ sống sót giữa các nhóm được so sánh.
Trong ví dụ này, đường cong sinh tồn Kaplan-Meier được tính toán cho mỗi nhóm giới tính và mức độ suy giảm chức năng phổi khác nhau trong dữ liệu ung thư phổi như sau:
logrank_test1 <- survdiff(with(lung, Surv(time, status)) ~ sex, data = lung)
logrank_test2 <- survdiff(with(lung, Surv(time, status)) ~ ph.ecog, data = lung)
print(logrank_test1)## Call:
## survdiff(formula = with(lung, Surv(time, status)) ~ sex, data = lung)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138 112 91.6 4.55 10.3
## sex=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001print(logrank_test2)## Call:
## survdiff(formula = with(lung, Surv(time, status)) ~ ph.ecog,
## data = lung)
##
## n=227, 1 observation deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ph.ecog=0 63 37 54.153 5.4331 8.2119
## ph.ecog=1 113 82 83.528 0.0279 0.0573
## ph.ecog=2 50 44 26.147 12.1893 14.6491
## ph.ecog=3 1 1 0.172 3.9733 4.0040
##
## Chisq= 22 on 3 degrees of freedom, p= 7e-05Kết quả kiểm định:
Nhóm giới tính (male-female): giá trị p= 0.001 < 0.05 nên có thể kết luận rằng có sự khác biệt đáng kể giữa nhóm giới tính nam và nữ về tỉ lệ sống sót.
Nhóm Mức độ suy giảm phổi của bệnh nhân (0 = Tốt, 1 = Tàn phế nhẹ, 2 = Tàn phế trung bình, 3 = Tàn phế nặng): giá trị p= 7e-05 < 0.05 nên có thể kết luận rằng có sự khác biệt đáng kể giữa các nhóm mức độ về tỉ lệ sống sót.
Mô hình hồi quy Cox
Mô hình Cox (Cox proportional hazards model) là
một phương pháp hồi quy theo tỉ lệ nguy cơ, được sử dụng trong
phân tích sống còn để đánh giá tác động của các biến độc lập đến tỉ lệ
sống sót. Trong R, chúng ta có thể sử dụng hàm coxph()
trong package survival để tạo một mô hình
Cox. Để tạo một mô hình Cox trên R, tôi dùng hàm sau:
coxph(surv_obj ~ x1 + x2 + ..., data = mydata)Trong đó:
surv_obj: là đối tượng sinh tồn được tạo từ hàm
surv()(phần 4.4)x: là các yếu tố được đưa vào mô hình để phân nhóm các đường cong sinh tồn khác nhau.
Tiếp tục với dữ liệu ung thư phổi - lung theo biến giới tính, chúng ta có mô hình Cox như sau:
coxfit <- coxph(Surv(time, status) ~ sex, data = lung)
coxfit## Call:
## coxph(formula = Surv(time, status) ~ sex, data = lung)
##
## coef exp(coef) se(coef) z p
## sex -0.5310 0.5880 0.1672 -3.176 0.00149
##
## Likelihood ratio test=10.63 on 1 df, p=0.001111
## n= 228, number of events= 165Kiểm định sự phù hợp của mô hình Cox
Trong R, để kiểm tra tính phù hợp của mô hình Cox với giả định độc lập về thời gian (proportional hazards assumption), chúng ta có thể sử dụng hàm cox.zph() trong gói survival.
Cú pháp sử dụng hàm cox.zph() như sau:
cox.zph(model, transform = "km")Trong đó:
model: là đối tượng mô hình Cox đã được ước lượng bằng hàm coxph() trong gói survival.
transform: là phương pháp biến đổi thời gian (time transformation) được sử dụng để tính toán các thống kê kiểm định. Giá trị mặc định là “km” tương ứng với phương pháp đường chéo Kaplan-Meier.
Các biến độc lập được kiểm tra bằng cách tính toán hệ số tương quan giữa các giá trị dự đoán của hàm Cox và các hàm đường chéo Kaplan-Meier tương ứng với từng nhóm giá trị của biến độc lập đó. Nếu giá trị p của kiểm định nhỏ hơn một ngưỡng xác định (thường là 0.05), ta có thể bác bỏ giả định độc lập về thời gian và cần phải tìm cách điều chỉnh mô hình Cox để phù hợp hơn với dữ liệu.
Ví dụ với mô hình Cox đã được tạo từ hàm coxph() (phần 4.8), ta kiểm định như sau:
cox.zph(coxfit, transform = "km")## chisq df p
## sex 2.86 1 0.091
## GLOBAL 2.86 1 0.091Kết quả cho thấy: giá trị p= 0.091 > 0.05 nên mô hình Cox được tạo trên là phù hợp.
So sánh tỉ lệ tử vong giữa hai nhóm bằng mô hình Cox
Từ mô hình Cox, tôi có thể thu được các bảng kết
quả kiểm định giữa hai nhóm giới tính bằng cách sử dụng
hàm tbl_regression() từ package gtsummary, với
tùy chọn exp = TRUE để trả về tỷ lệ rủi ro:
library("gtsummary")
coxph(Surv(time, status) ~ sex, data = lung) %>% tbl_regression(exp = TRUE) | Characteristic | HR1 | 95% CI1 | p-value |
|---|---|---|---|
| sex | 0.59 | 0.42, 0.82 | 0.001 |
| 1 HR = Hazard Ratio, CI = Confidence Interval | |||
HR biểu thị tỷ lệ nguy cơ giữa hai nhóm tại bất kỳ thời điểm cụ thể nào. HR < 1 cho thấy nguy cơ tử vong giảm trong khi HR > 1 cho thấy nguy cơ tử vong tăng lên.
Kết quả từ hàm tbl_regression() cho thấy HR =
0,59 ngụ ý rằng số nam giới sắp chết nhiều gấp
0,59 lần so với nữ giới, tại bất kỳ thời điểm nào. Nói cách
khác, phụ nữ có nguy cơ tử vong thấp hơn đáng kể so với nam giới trong
các dữ liệu này.
Dự đoán tỉ lệ sống sót trong tương lai
Để dự đoán tỉ lệ sống sót của một cá nhân dựa trên mô hình Cox đã
được ước lượng trên R, chúng ta có thể sử dụng hàm
predict() trong package survival.
Cú pháp sử dụng hàm predict() để dự đoán tỉ lệ sống sót
của một cá nhân như sau:
predict(object, newdata, type = "survival")Trong đó:
object: Là đối tượng mô hình Cox đã được ước lượng bằng hàm
coxph()trong package survival.newdata: Là bộ dữ liệu chứa thông tin về các biến độc lập của biến cần dự đoán. Bộ dữ liệu này phải có cùng tên và định dạng với bộ dữ liệu được sử dụng để ước lượng mô hình.
type: Là loại dự đoán cần thực hiện. Trong trường hợp này, ta sử dụng
type = "survival"để dự đoán tỉ lệ sống sót.
Ví dụ: tôi muốn sự báo tỉ lệ sống sót của nhóm giới tính Nam vào ngày thứ 1050, tôi thực hiện các bước dự báo như sau:
mod <- coxph(Surv(time,status) ~+ sex, data = cancer)
pred_dat <- data.frame(time = c(1050), status = c(0), sex = c(0))
preds <- predict(mod, newdata = pred_dat,
type = "survival", se.fit = TRUE)
pred_dat$prob <- preds$fit
pred_dat$lcl <- preds$fit - 1.96*preds$se.fit
pred_dat$ucl <- preds$fit + 1.96*preds$se.fit
pred_dat## time status sex prob lcl ucl
## 1 1050 0 0 0.002573077 -0.005221542 0.0103677Kết quả: vào ngày 1050, với trạng thái sóng sốt (0) và giới tính Nam (0) ta có tỉ lệ sống sót là 0.00257