Phân Tích Survival Trong Trường Hợp Biến Cố Xảy Ra Nhiều Lần - Time-Dependent Covariates in Survival Analysis
1 Survival Analysis: Time - event analysis
Phân tích survival rất thường gặp trong nhiều nghiên cứu y học. Đó là loại phân tích dữ liệu về một biến cố xảy ra tại những thời điểm được ghi nhận, sau một can thiệp điều trị nào đó hoặc sau thời điểm bệnh lí được chẩn đoán.
Ví dụ, tính toán nguy cơ tử vong của những bệnh nhân sau khi được phẩu thuật cắt bỏ khối ung thư tuyến giáp. Người ta theo dõi những bệnh nhân được thực hiện phẩu thuật trên sau một khoảng thời gian 3-10 năm, thống kê những trường hợp tử vong để xác định tỉ lệ tử vong, khả năng sống của bệnh nhân tại các thời điểm sau phẩu thuật.
Bên cạnh đó, người ta còn so sánh xác suất sống của bệnh nhân sau khi được can thiệp của hai hay nhiều phương thức điều trị nào đó, ví dụ phẩu thuật so với phẩu thuật kết hợp hóa trị.
Ngoài ra, người ta còn tìm ra mối liên hệ giữa xác suất còn sống của bệnh nhân với các yếu tố khác như tuổi lúc can thiệp, giai đoạn bệnh, liều lượng hóa trị thông qua việc xây dựng mô hình tỉ lệ nguy cơ Cox.
Thông thường trong survival analysis người ta theo dõi một biến cố (event) có xảy ra hay không và thời gian tính từ khi can thiệp cho đến khi xảy ra biến cố đó (time). Ví dụ biến cố tử vong, chỉ xảy ra một lần đối với một bệnh nhân.
Với những biến cố có thể xảy ra nhiều lần, chẳng hạn cơn kịch phát hen nặng có thể xảy ra vài lần trong thời gian theo dõi nghiên cứu. Trong trường hợp này người ta theo dõi biến cố đầu tiên và thời gian đến khi xảy ra biến cố đầu tiên này.
Nếu chỉ sử dụng dữ liệu của lần xảy ra biến cố đầu tiên mà thôi thì rất có thể không phản ánh được đúng mức độ quan trọng của nguy cơ, xác suất sống còn khác nhau giữa các bệnh nhân. Có những bệnh nhân không xảy ra cơn kịch phát hen nặng, nhưng có những bệnh nhân khác đã xảy ra một lần, hai lần, ba lần hoặc nhiều hơn nữa. Các cơn kịch phát hen nặng đã xảy ra thưa hay gần nhau sẽ phản ánh những mức độ nguy cơ khác nhau và mối liên hệ khác nhau với các biến số khác.
Một bệnh nhân chỉ xảy ra biến cố 1 lần thì nguy cơ sẽ khác nhau rất nghiều so với các bệnh nhân có 3-4 biến cố xảy ra gần nhau trong thời gian nghiên cứu.
Bài viết này muốn trình bày cách xử lí dữ liệu survival trong trường hợp biến cố có thể xảy ra nhiều lần và số lượng biến cố khác nhau giữa các đối tượng nghiên cứu.
Tuy nhiên, trước khi trình bày đến vấn đề trên, tôi vẫn phải phân tích theo tình huống biến cố chỉ xảy ra một lần như thông lệ để có sự so sánh.
2 Chuẩn bị dữ liệu
Dữ liệu sử dụng để phân tích trong bài này được lấy từ nguồn OpenML (phpOAYun7.csv), có chỉnh sửa để phù hợp với chủ đề của bài viết.
Bộ dữ liệu này có 120 rows và 14 variables, trong đó có 1 biến số event là Class, thời gian đến khi xảy ra biến cố đầu tiên (time), tổng thời gian theo dõi (futime), và thời gian xảy ra các biến cố từ etime1 đến etime4. Giá trị NA về thời gian chỉ thời điểm không xảy ra biến cố.
library(readr)
df <- read_csv("D:/Data/survival_df.csv")
head(df)## # A tibble: 6 x 14
## ID V1 V2 V3 V4 V5 V6 Class time futime etim~ etim~
## <int> <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <int> <int> <int> <int>
## 1 1 35.5 no no no no no no 365 365 NA NA
## 2 2 35.9 no no yes yes yes no 365 365 NA NA
## 3 3 35.9 no no no no no no 365 365 NA NA
## 4 4 36.0 no no yes yes yes no 365 365 NA NA
## 5 5 36.0 no no no no no no 365 365 NA NA
## 6 6 36.0 no no no no no no 365 365 NA NA
## # ... with 2 more variables: etime3 <int>, etime4 <int>3 Phân tích survival dựa vào thời gian đến khi xảy ra biến cố đầu tiên
Trước hết chúng ta load các packages cần thiết cho phân tích survival và khai báo biến số.
library(tidyverse)
library("survival")
library("survminer")
df$Class <- ifelse(df$Class == "yes", 1 ,0)
fit <- survfit(Surv(time, Class) ~ 1, data = df)
ggsurvplot(fit,
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800"))
V2 đến V6 là những biến số factor, chúng ta phân tích xem xác suất xảy ra biến cố (Class) có khác biệt giữa các nhóm của những biến số này hay không.
Trước tiên với V2
fit2 <- survfit(Surv(time, Class) ~ V2, data = df)
ggsurvplot(fit2,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"))
Khác biệt xác suất không xảy ra biến cố giữa hai nhóm “no” và “yes” là có ý nghĩa thống kê. Tương tự như vậy với các biến số từ V3 đến V6.
Ngoại trừ với V4 và V6 là không có sự khác biệt có ý nghĩa.
fit6 <- survfit(Surv(time, Class) ~ V6, data = df)
ggsurvplot(fit6,
pval = TRUE, conf.int = TRUE,
risk.table = TRUE, # Add risk table
risk.table.col = "strata", # Change risk table color by groups
linetype = "strata", # Change line type by groups
surv.median.line = "hv", # Specify median survival
ggtheme = theme_bw(), # Change ggplot2 theme
palette = c("#E7B800", "#2E9FDF"))
3.1 Cox proportional hazards models
Ngoài các biến số factor, thì V1 là biến numeric. Chúng ta xem xét model Cox biểu hiện sự tương quan giữa xác suất không có biến cố của các đối tượng với các biến số này.
Coxfit <- coxph(Surv(time, Class) ~ V1 + V2 + V3 + V4+ V5 + V6, data = df)
summary(Coxfit)## Call:
## coxph(formula = Surv(time, Class) ~ V1 + V2 + V3 + V4 + V5 +
## V6, data = df)
##
## n= 120, number of events= 51
##
## coef exp(coef) se(coef) z Pr(>|z|)
## V1 1.1840 3.2674 0.2161 5.478 4.31e-08 ***
## V2yes 1.8930 6.6394 0.5801 3.263 0.001102 **
## V3yes 3.1759 23.9487 0.9014 3.523 0.000426 ***
## V4yes -0.1036 0.9015 0.4237 -0.245 0.806774
## V5yes -1.6625 0.1897 0.5296 -3.139 0.001694 **
## V6yes -0.2057 0.8141 0.4366 -0.471 0.637584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## V1 3.2674 0.30605 2.13905 4.9909
## V2yes 6.6394 0.15062 2.12962 20.6991
## V3yes 23.9487 0.04176 4.09279 140.1341
## V4yes 0.9015 1.10920 0.39291 2.0686
## V5yes 0.1897 5.27258 0.06717 0.5355
## V6yes 0.8141 1.22837 0.34595 1.9157
##
## Concordance= 0.925 (se = 0.042 )
## Rsquare= 0.71 (max possible= 0.979 )
## Likelihood ratio test= 148.6 on 6 df, p=0
## Wald test = 44.09 on 6 df, p=7.105e-08
## Score (logrank) test = 131.2 on 6 df, p=0Kết quả này phù hợp với phép phân tích ở trên đối với v4 và V6 không có ý nghĩa thống kê. Khi loại hai biến số này ra chúng ta có mô hình Cox như sau:
Coxfit_2 <- coxph(Surv(time, Class) ~ V1 + V2 + V3 + V5, data = df)
summary(Coxfit_2)## Call:
## coxph(formula = Surv(time, Class) ~ V1 + V2 + V3 + V5, data = df)
##
## n= 120, number of events= 51
##
## coef exp(coef) se(coef) z Pr(>|z|)
## V1 1.1546 3.1729 0.2162 5.340 9.29e-08 ***
## V2yes 1.8663 6.4642 0.5820 3.207 0.001342 **
## V3yes 2.8852 17.9074 0.8086 3.568 0.000360 ***
## V5yes -1.4853 0.2264 0.4478 -3.317 0.000911 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## V1 3.1729 0.31517 2.07686 4.8474
## V2yes 6.4642 0.15470 2.06607 20.2247
## V3yes 17.9074 0.05584 3.67040 87.3682
## V5yes 0.2264 4.41614 0.09414 0.5447
##
## Concordance= 0.925 (se = 0.042 )
## Rsquare= 0.709 (max possible= 0.979 )
## Likelihood ratio test= 148 on 4 df, p=0
## Wald test = 43.03 on 4 df, p=1.018e-08
## Score (logrank) test = 128.5 on 4 df, p=0Khi loại được hai biến số không đóng góp nhiều cho mô hình, Rsquare giảm không đáng kể là điều chấp nhận được.
Tuy nhiên, mô hình Cox này đang sử dụng dữ liệu từ biến cố đầu tiên mà thôi. Có thể việc loại trừ dữ liệu của những biến cố khác đã xảy ra trong thời gian nghiên cứu làm cho các mối tương quan chưa được phản ánh đúng thực chất.
Bây giờ ta xây dựng mô hình Cox có sử dụng dữ liệu của tất cả những biến cố đã xảy ra cho các đối tượng nghiên cứu.
4 Mô hình Cox cho Time-dependent covariates
Để phân tích được dữ liệu đồng thời nhiều biến cố, chúng ta phải:
- Sắp xếp lại dữ liệu để mỗi biến cố thành một row
- Thêm những argument cần thiết cho hàm Cox
Ở đây chúng ta sử dụng hàm tmerg để merge dữ liệu theo từng biến cố đã xảy ra. Một đối tượng có k biến cố thì sẽ có k dòng.
Số biến cố lớn nhất là k = 4, tức là một bệnh nhân nghiên cứu đã có số biến cố nhiều nhất là k = 4. Vì thế sẽ có 4 biến số mới tương ứng với 4 biến cố là etime1 đến etime4.
newdf <- tmerge(data1 = df[,1:10], data2 = df, id=ID, tstop=futime,
elapse = event(etime1), elapse = event(etime2),
elapse = event(etime3), elapse = event(etime4))
newdf <- tmerge(newdf, newdf, ID, enum=cumtdc(tstart))
dim(newdf)## [1] 202 15tail(newdf)## ID V1 V2 V3 V4 V5 V6 Class time futime id tstart tstop elapse
## 197 119 41.5 yes yes yes no yes 1 24 350 119 221 350 0
## 198 120 41.5 yes yes yes no yes 1 17 352 120 0 17 1
## 199 120 41.5 yes yes yes no yes 1 17 352 120 17 80 1
## 200 120 41.5 yes yes yes no yes 1 17 352 120 80 228 1
## 201 120 41.5 yes yes yes no yes 1 17 352 120 228 300 1
## 202 120 41.5 yes yes yes no yes 1 17 352 120 300 352 0
## enum
## 197 4
## 198 1
## 199 2
## 200 3
## 201 4
## 202 5Khoảng thời gian theo dõi của mỗi bệnh nhân được chia ra thành các giai đoạn tương ứng với các biến cố đã xảy ra. Mỗi giai đoạn này có tstart là điểm khởi đầu và tstop là thời điểm kết thúc. Tương ứng với các giai đoạn, elapse bằng 1 khi có biến cố xảy ra và elapse = 0 khi không có biến cố xảy ra trong giai đoạn tương ứng.
Bộ dữ liệu mới này, newdf, có 202 dòng và 15 biến số.
Bây giờ chúng ta viết code cho model Cox với dữ liệu của tất cả các biến cố, không chỉ biến cố đầu tiên, được đưa vào mô hình để phân tích.
Coxfit2 <- coxph(Surv(tstart, tstop,elapse) ~ V1 + V2 + V3 + V4+ V5 + V6 + cluster(ID), newdf)
summary(Coxfit2)## Call:
## coxph(formula = Surv(tstart, tstop, elapse) ~ V1 + V2 + V3 +
## V4 + V5 + V6 + cluster(ID), data = newdf)
##
## n= 202, number of events= 82
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## V1 0.53523 1.70785 0.11163 0.06485 8.253 < 2e-16 ***
## V2yes 0.51367 1.67142 0.33745 0.21533 2.385 0.01706 *
## V3yes 2.23139 9.31284 0.66445 0.73542 3.034 0.00241 **
## V4yes -0.06073 0.94107 0.32534 0.17465 -0.348 0.72803
## V5yes -0.49625 0.60881 0.32079 0.15342 -3.234 0.00122 **
## V6yes -0.09448 0.90985 0.34895 0.17951 -0.526 0.59868
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## V1 1.7078 0.5855 1.5040 1.9393
## V2yes 1.6714 0.5983 1.0960 2.5490
## V3yes 9.3128 0.1074 2.2034 39.3615
## V4yes 0.9411 1.0626 0.6683 1.3252
## V5yes 0.6088 1.6425 0.4507 0.8224
## V6yes 0.9098 1.0991 0.6400 1.2935
##
## Concordance= 0.817 (se = 0.032 )
## Rsquare= 0.457 (max possible= 0.979 )
## Likelihood ratio test= 123.5 on 6 df, p=0
## Wald test = 186.7 on 6 df, p=0
## Score (logrank) test = 100.5 on 6 df, p=0, Robust = 60.91 p=2.943e-11
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).Kết quả cũng cho thấy các biến số đều có ý nghĩa ngoại trừ V4 và V6. Tuy nhiên, so với mô hình ban đầu chỉ với biến cố đầu tiên thì Rsquare đã giảm từ 0.709 xuống còn 0.457. Rất có thể mô hình này phản ánh đúng thực chất các mối quan hệ trong bộ dữ liệu trên.
Để có thể rút gọn mô hình chúng ta loại bỏ hai biến số V4 và v6 và phân tích tiếp tục.
Coxfit3 <- coxph(Surv(tstart, tstop,elapse) ~ V1 + V2 + V3 + V5 + cluster(ID), newdf)
summary(Coxfit3)## Call:
## coxph(formula = Surv(tstart, tstop, elapse) ~ V1 + V2 + V3 +
## V5 + cluster(ID), data = newdf)
##
## n= 202, number of events= 82
##
## coef exp(coef) se(coef) robust se z Pr(>|z|)
## V1 0.52261 1.68642 0.10882 0.06506 8.033 9.99e-16 ***
## V2yes 0.51137 1.66757 0.33675 0.21634 2.364 0.01809 *
## V3yes 2.15440 8.62271 0.64900 0.76283 2.824 0.00474 **
## V5yes -0.42181 0.65586 0.26116 0.13078 -3.225 0.00126 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## V1 1.6864 0.5930 1.4845 1.9158
## V2yes 1.6676 0.5997 1.0913 2.5482
## V3yes 8.6227 0.1160 1.9334 38.4561
## V5yes 0.6559 1.5247 0.5076 0.8475
##
## Concordance= 0.816 (se = 0.032 )
## Rsquare= 0.457 (max possible= 0.979 )
## Likelihood ratio test= 123.2 on 4 df, p=0
## Wald test = 185.3 on 4 df, p=0
## Score (logrank) test = 98.9 on 4 df, p=0, Robust = 60.81 p=1.965e-12
##
## (Note: the likelihood ratio and score tests assume independence of
## observations within a cluster, the Wald and robust score tests do not).Tinh gọn mô hình, Rsquare không thay đổi (0.457), cũng có thể là một chọn lựa chấp nhận được.
Cần có những datasets phù hợp hơn với chủ đề này để phân tích và đánh giá kết quả của các phân tích này.
Cám ơn các bạn đã đọc bài viết. Mong nhận được những góp ý, phân tích của các bạn.