Kaplan meier là gì

PHÂN TÍCH SỐNG CÒN (Survival Analysis) Lê Tấn Phùng* NHỮNG KHÁI NIỆM CƠ BẢN Phân tích sống còn là gì? Trong nghiên cứu, đôi khi người ta gặp những số liệu liên quan đến đánh giá một sự kiện nào đó theo thời gian. Ví dụ, sự kiện “tử vong” mô tả thời gian từ khi chẩn đoán đến khi chết. Sự kiện “tái phát” mô tả thời gian từ khi áp dụng một liệu pháp điều trị nào đó cho đến khi xảy ra tái phát, chẳng hạn sử dụng phương đặt sten mạch vành và theo dõi cho đến khi tái hẹp (tái phát). Tuy nhiên, cho đến khi kết thúc thời gian nghiên cứu, không phải tất cả các trường hợp đều xảy ra “sự kiện”. Ví dụ nghiên cứu về tử vong do ung thư gan từ lúc được chẩn đoán, thời gian nghiên cứu trong 24 tháng, không phải tất cả bệnh nhân nghiên cứu đều tử vong sau khi kết thúc nghiên cứu và do đó thời gian tử vong của họ là chưa biết. Do đó số liệu nghiên cứu loại này có đặc điểm hiếm khi có phân phối bình thường (normal distribution) vì một số số liệu hoàn chỉnh, một số số liệu bị cắt (censored). Số liệu của dạng nghiên cứu này được phân tích theo phương pháp riêng gọi là phân tích sống còn (survival analysis). a) Phân tích sống còn (viết tắt là SA) liên quan đến biến thời gian. Biến này ghi nhận thời gian từ lúc bắt đầu theo dõi cho đến khi xảy ra biến cố (Trong tiếng Anh thường sử dụng 2 thuật ngữ: time-to-failure hoặc time-to-event). Các thuật ngữ khác của SA có thể dùng trong tiếng Anh gồm: duration analysis, transition analysis, failure time analysis. - Không giống như hồi quy tuyến tính, SA có biến phụ thuộc là biến nhị biến. Có nghĩa là biến này chỉ có 2 giá trị. Ví dụ: chết hoặc sống, thành công hay thất bại, tái phát hay không tái phát v.v... - Không giống như hồi quy logistic, SA liên quan đến việc xây dựng mô hình biến thời gian – biến cố (modelling the time to an event). - SA liên quan đến áp dụng các biểu đồ Kaplan Meier (KM), Log rank test, và hồi quy Cox (sẽ đề cập đến ở phần sau). b) Những ưu điểm, thuận lợi của SA - Có thể giải thích cho những số liệu bị cắt (censored data). Khái niệm này sẽ được giải thích ở phần sau. - Có thể so sánh giữa 2 nhóm hoặc nhiều nhóm với nhau - Có thể đánh giá mối liên quan giữa các biến độc lập và thời gian sống (survival time). *

Bác sĩ nội khoa, Tiến sĩ YTCC, Sở Y tế tỉnh Khánh Hòa. [email protected]

1 Số liệu về thời gian – biến cố (time-to-event) Trong nghiên cứu y học, rất hay gặp những trường hợp biến số liên quan đến thời gian. Ví dụ: thời gian chết, thời gian khỏi bệnh sau khi điều trị, thời gian tái phát v.v... Tất cả các số liệu liên quan đến những biến số như vậy gọi là số liệu sống còn (survival data), mặc dù thuật ngữ này có vẻ mô tả không chính xác bản chất của số liệu. Các số liệu dạng “sống còn” rất hay gặp trong nghiên cứu lâm sàng nên SA đôi khi được gọi dưới tên là “Nghiên cứu lâm sàng” (Clinical Research), mặc dù số liệu sống còn vẫn có thể có trong các nghiên cứu khác. Khi nào thì sử dụng SA? - SA được sử dụng trong những phân tích liên quan đến các số liệu về thời gian Khi người ta cho rằng một hoặc nhiều biến độc lập là lý do làm cho có sự khác biệt về thời gian đối với một biến cố. Đặc biệt, SA được sử dụng khi thời gian theo dõi không được hoàn thành hoặc có sự khác nhau về thời gian theo dõi này. Những vấn đề gặp phải khi sử dụng SA - Số liệu bị cắt hay bị xén (Data censoring and truncation): Đây là các số liệu không hoàn chỉnh. Phần sau sẽ giải thích nội dung này. Xác suất có điều kiện (conditional probability): Khi xem xét xác suất để cho một trường hợp “bị chết” tại một thời diểm nào đó, ta cần phải đặt điều kiện rằng trường hợp đó “còn sống” cho đến ngay trước thời điểm đó. Ví dụ, để tính xác suất cho A sẽ “chết” vào ngày 30 thì phải đặt điều kiện rằng đối tượng đó sẽ “còn sống” cho đến ngày 29. Số liệu bị cắt (Censoring) Số liệu bị cắt là những số liệu không hoàn chỉnh (incompleted data). Người ta phân biệt 2 loại số liệu bị cắt: Cắt bên phải và cắt bên trái a) Cắt bên phải (Right censored data): - Đây là dạng bị cắt phổ biến nhất trong SA. Hầu hết các số liệu không hoàn chỉnh trong SA đều thuộc dạng này. - Dạng số liệu này xảy ra khi một đối tượng nào đó vẫn “còn sống” khi nghiên cứu đã kết thúc, hoặc bị mất theo dõi (lost-to-follow-up), hoặc bỏ cuộc (withdrawn). Ví dụ, trong nghiên cứu tái phát nhồi máu cơ tim khi điều trị với loại thuốc A nào đó trong thời gian 1 năm, những số liệu bị cắt bên phải là những số liệu sau: Bệnh nhân không bị tái phái sau khi theo dõi đủ 1 năm; bệnh nhân theo dõi không đủ 1 năm vì bệnh nhân di chuyển chổ ở hay vì lý do gì đó mà không thể tiếp tục theo dõi được; hoặc bệnh nhân không tham gia nghiên cứu này nữa khi chưa đủ 1 năm theo dõi. b) Cắt bên trái (Left Censored data): Rất ít gặp. Dạng số liệu này khi “biến cố” đã xảy ra, trước khi nghiên cứu bắt đầu. Ví dụ, nghiên cứu về theo dõi sự xuất hiện của ung thư 2 gan, số liệu bị cắt bên trái nếu bệnh nhân tham gia nghiên cứu này đã bị ung thư gan từ trước (mà vì lý do nào đó nhà nghiên cứu không phát hiện được khi chọn nhóm nghiên cứu). Ví dụ về số liệu bị cắt Sơ đồ của 5 trường hợp dưới đây diễn tả một nghiên cứu kéo dài 12 tháng với 5 đối tượng theo dõi, từ A đến E. Các trường hợp này có ý nghĩa như sau: Tháng 4 2 T=5 A 6 8 10 x T = 12 B Hết NC T=3 C Rút khỏi NC T=8 D T=6 E 12 Hết NC Mất theo dõi x F Các đối tượng không xảy ra “biến cố”: B, C, D, E Theo dõi không hoàn thành:   Mất theo dõi: E Không tham gia nghiên cứu nữa: C Bị cắt bên phải:   B và D (vì nghiên cứu đã kết thúc); C và E (vì mất theo dõi hay rút khỏi nghiên cứu) Phân phối thời gian sống còn (distribution of survival time) Một số khái niệm và ký hiệu sau đây sử dụng trong SA: f (t) ký hiệu cho xác suất sống còn tại một thời điểm t nào đó F (t) ký hiệu cho xác suất sống còn cộng dồn (cumulative probability) 3 Ví dụ: f (23) là xác suất sống còn tại ngày thứ 23; còn F(23) là xác suất sống còn đến ngày thứ 23 và cả ngày 23. Tuy nhiên, trong SA, người ta quan tâm nhiều hơn đến xác suất mà đối tượng nào đó còn “sống” khi qua khỏi một mốc thời gian nhất định. Do đó, hàm số quan tâm ở đây là hàm sống còn (Survival function) S(t), và hàm Hazard h(t). Hàm sống còn S(t) xác định xác suất “còn sống” (tức là không xảy ra sự kiện) lâu hơn thời điểm t. Ví dụ, xác suất sống lâu hơn 1 năm.  S (t )  P(T  t )   f (u) du  1  F (t ) t Hàm Hazard, h(t) là nguy cơ tức thì (instantaneous) của sự kiện tại thời điểm t h(t )  f (t ) S (t ) Ví dụ liên quan đến sự kiện tử vong do ung thư phổi. Hàm sống còn S (t) sau 2 năm tức là xác suất còn sống sau 2 năm của một người bị ung thư phổi, trong khi hàm hazard mô tả xác suất để người đó chết tại thời điểm 2 năm sau. Nói cách khác, hàm sống còn S (t) phản ánh tỉ lệ dồn của việc không xảy ra sự kiện (cumulative non-occurence) sau thời điểm t trong khi hàm hazard h (t) phản ánh tỉ lệ xảy ra sự kiện tại thời điểm t xác định nào đó (incidence). Xác suất sống còn S (t) được tính bằng tích số của các xác suất sống còn tại thời điểm 1, 2, 3, ... đến thời điểm t, tức là S (t) = S(1)* S(2)* S(3)*...* S(t). Hàm sống còn S (t) được biểu diễn bởi biểu đồ Kaplan Meier. Biểu đồ Kaplan Meier có thể áp dụng để thể hiện xác suất sống còn giữa 2 nhóm khác nhau, ví dụ như giữa nam và nữ, giữa phác đồ A và phác đồ B v.v... Ngoài ra, để so sánh 2 nhóm với nhau trong phân tích sống còn, Log-rank test được sử dụng với giả thuyết H0 rằng xác suất sống còn của 2 nhóm là không khác nhau thông qua test χ2 đi kèm với giá trị p. PHÂN TÍCH SỐNG CÒN Đặc điểm: SA sẽ cho những kết quả phân tích sau đây: - Các thống kê mô tả thông thường (trung bình, mode, độ lệch chuẩn, phương sai...) Phân tích nhị biến (univariate) Phân tích đa biến (ví dụ như Hồi quy Cox) Đường biểu diễn sống còn Kaplan-Meier - Kaplan-Meier là ước lượng (estimator) của đường biểu diễn sống còn 4 - Nó còn được gọi là công thức product-limit (product limit formula) Là phương pháp phi tham số (non-parametric), do đó không cần giả định phân phối (distribution assumptions) cho biến thời gian. Có thể giải thích cho các số liệu bị cắt Tạo nên một biểu đồ dạng “bậc thang” đặc trưng Không thể lý giải cho nhiễu hay tương tác (confounding or effect modification) do các biến độc lập khác Đường biểu diễn được giải thích là xác suất “còn sống” (tức là không xảy ra biến cố) tại một thời điểm t nào đó. Giải thích này sẽ được nói rõ hơn trong ví dụ tiếp theo. Ví dụ minh họa sử dụng R Bộ số liệu ovarian sẽ được sử dụng để minh họa cho SA sử dụng R. Bộ số liệu mô tả kết quả nghiên cứu của 26 bệnh nhân bị ung thư buồng trứng với các biến số sau đây: - futime: Làn thời gian còn sống, tính bằng ngày fustat: Là tình trạng có hay không dữ liệu bị cắt (censored), với 1 là bị cắt và zero là không bị cắt age: Tuổi, tính bằng năm resid.ds: Khối u có thoái triển hay không, với 1 là không thoái triển và 2 là có thoái triển rx: Phác đồ điều trị 1 và 2 ecog.ps: Đánh giá bệnh nhân theo thang đo ECOG, với 1 là tốt và 2 là không tốt Bộ số liệu này gắn liền với package “survival”. Ngoài ra, để có thể vẽ biểu đồ Kaplan Meier cần có thêm package “survminer”. Để có thể mô tả bộ số liệu, package “dplyr” nên được sử dụng. Như vậy, để có thể bắt đầu phân tích, 3 packages này phải được đưa vào. > library(survival) > library(survminer) Loading required package: ggplot2 Loading required package: ggpubr Loading required package: magrittr > library(dplyr) Attaching package: ‘dplyr’ The following objects are masked from ‘package:stats’: filter, lag The following objects are masked from ‘package:base’: intersect, setdiff, setequal, union Gọi bộ số liệu ovarian: > data(ovarian) 5 Mô tả các biến trong bộ số liệu ovarian thông qua lệnh glipmse: > glimpse(ovarian) Observations: 26 Variables: 6 $ futime 769,... $ fustat 0, 0, ... $ age 56.9370, ... $ resid.ds 2, 1, ... $ rx 2, 2, ... $ ecog.ps 1, 2, ... 59, 115, 156, 421, 431, 448, 464, 475, 477, 563, 638, 744, 1, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 72.3315, 74.4932, 66.4658, 53.3644, 50.3397, 56.4301, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 2, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 2, 2, 2, 1, 2, 2, 1, 2, 1, 1, 2, 2, 1, 1, Kết quả cho thấy bộ số liệu có 26 quan sát, 6 biến với các tên ở bên trái và một số giá trị đại diện ở bên phải. Để thuận tiên cho việc phân tích theo 2 nhóm sử dụng Log-rank test, một số biến sẽ được mã hóa thành biến phân loại 2 giá trị (nhị biến), còn gọi là dummy variable. Trước tiên, xem xét biến age, biến này sẽ chia tuổi ra làm 2 nhóm. Để có căn cứ chia biến này thành 2 nhóm, trước hết vẽ biểu đồ histogram cho biến này, kết quả như sau: > hist(ovarian$age) Hình trên cho thấy tuổi được phân bố 2 nhóm rõ rệt với điểm cắt là 50 tuổi. Do đó, các biến được mã hóa trở lại với hàm cơ bản factor như sau: > ovarian$age_group = 50,"old", "young")) 6 > ovarian$rx ovarian$resid.ds ovarian$ecog.ps # Tao object sống còn, trong đó biến thời gian là futime và biến sự kiện là fustat > surv_object # Tạo object fit0 để tạo hàm Kaplan Meier không phân nhóm (biến nhóm =1) > fit0 # Tạo object fit1 để tạo hàm Kaplan Meier theo nhóm điều trị rx > fit1 #Thống kê cơ bản cho 2 phân tích sống còn fit0 và fit1 > fit0 Call: survfit(formula = surv_object ~ 1, data = ovarian) n events 26 12 median 0.95LCL 0.95UCL 638 464 NA > fit1 Call: survfit(formula = surv_object ~ rx, data = ovarian) n events median 0.95LCL 0.95UCL rx=A 13 7 638 268 NA rx=B 13 5 NA 475 NA > # Vẽ biêu đồ Kaplan Meier chung > ggsurvplot(fit0,data = ovarian) 7 > # Vẽ biêu do Kaplan Meier theo nhom dieu tri rx, có hiển thị giá trị p > ggsurvplot(fit1,data = ovarian, pval= TRUE) Các vạch dọc trên biểu đồ đại diện cho dữ liệu bị cắt. Nếu tồn tại sự khác nhau có ý nghĩa thống kê thì 2 biểu đồ sẽ không cắt nhau (nhưng không phải trường hợp ngược lại là luôn luôn đúng) > # Test log rank > survdiff(surv_object ~ rx,data=ovarian) 8 Call: survdiff(formula = surv_object ~ rx, data = ovarian) N Observed Expected (O-E)^2/E (O-E)^2/V rx=A 13 7 5.23 0.596 1.06 rx=B 13 5 6.77 0.461 1.06 Chisq= 1.1 on 1 degrees of freedom, p= 0.303 Giải thích kết quả như sau: - - - Kết quả thống kê fit0 với trung vị bằng 638: Có 26 quan sát được phân tích với 12 sự kiện được xảy ra. Có 50% số phụ nữ ung thư buồng trứng sống được đến ít nhất 638 ngày. Kết quả thống kê fit1: Có 2 nhóm điều trị theo phác đồ A và phác đồ B. Ở phác đồ A có 13 phụ nữ và 7 trường hợp xảy ra sự kiện. Trung vị của nhóm phác đồ A là 268 nói lên rằng 50% số phụ nữ còn sống ít nhất đến 638 ngày. Phác đồ B có 13 phụ nữ và 5 trường hợp xảy ra sự kiện. Không thu được giá trị trung vị trong phác đồ này (đường biểu diễn Kaplan Meier thể hiện rõ điều này). Test Log-rank: Test cho kết quả χ2 = 1,1 và giá trị p = 0,303. Có thể kết luận hiệu quả điều trị giữa phác đồ A và B không khác nhau. Đường biểu diễn Kaplan Meier cũng thể hiện giá trị p này. PHÂN TÍCH SỐNG CÒN ĐA BIẾN Phần trước đã đề cập đến SA với Kaplan Meier, chỉ dành cho phân tích nhị biến, giữa 1 biến thời gian và một biến độc lập khác (ví dụ như biến giới tính ở ví dụ trên). Trong trường hợp SA với 2 hay nhiều biến độc lập, ta không thể dùng Kaplan Meier với Log rank test mà phải dùng một mô hình khác. Mô hình hay gặp nhất là hồi quy proportional hazard Cox (Cox proportional hazards Regression), hay gọi tắt là hồi quy Cox. Hồi quy Cox sử dụng ước lượng Maximum likelihood (MSE). Do đó, có thề sử dụng MSE để giải thích kết quả, tương tự như áp dụng trong hồi quy tuyến tính tổng quát hoá (generalized linear regression). Như trên đã nói, không giống như ước lượng Kaplan Meier phân tích xác suất “còn sống”, hồi quy Cox phân tích xác suất “chết” tại thời điểm t nào đó và xác suất này tương đương với xác suất điểm (incidence). Trở lại với ví dụ trên, người ta muốn phân tích ảnh hưởng đồng thời của phác đồ điều trị, sự thoái triển của khối u, nhóm tuổi và điểm ECOG đối với xác suất “chết” của phụ nữ ung thư buống trứng. Để thực hiện phân tích này, hồi quy proportional hazards Cox được áp dụng với hàm coxph và vẽ biểu đồ bằng hàm ggforest. Biểu đồ này sẽ thể hiện tỉ số harzard (Hazard 9 ratio: HR). Nói một cách tổng quát trong một mô hình đa biến, khi HR > 1 thì biểu hiện nguy cơ “chết” tăng. Ngược lại, nếu HR < 1 thì thể hiện nguy cơ “chết” giảm. Sau đây là các bước thực hiện trên R với bộ số liệu ovarian: > # Tạo object fit.coxph với hàm coxph để thực hiện mô hình hồi quy Cox > fit.coxph ggforest(fit.coxph, data = ovarian) Giải thích kết quả của biểu đồ: Mô hình gồm 4 biến độc lập là các biến nhị biến: phác đồ điều trị rx (A hoặc B), sự thoái triển của khối u resid.ds (có hay không), nhóm tuổi age_group (già hay trẻ), điểm số ECOG ecog.ps (tốt hay xấu). Mỗi biến độc lập có 1 giá trị tham chiếu, và giá trị HR là so sánh với giá trị tham chiếu này. Biến rx có giá trị tham chiếu là phác đồ A; biến resid.ds có giá trị tham chiếu là “không có thoái triển”; biến nhóm tuổi có giá trị tham chiếu là già; và biến ecog.ps có giá trị tham chiếu là “tốt”. Kết quả phân tích HR cho thấy: - Nguy cơ chết của bệnh nhân dùng phác đồ B thấp hơn phác đồ A vì HR = 0,25 < 1 và có ý nghĩa thống kê (p = 0,032) Nguy cơ chết của bệnh nhân có thoái triển khối u cao hơn so với không thoái triển vì HR = 4,25 > 1 và có ý nghĩa thống kê (p = 0,045). Không rõ vì sao lại có kết luận hơi ngược này. 10 - Nguy cơ chết của bệnh nhân trẻ thấp hơn so với bệnh nhân già vì HR = 0,11 < 1 và có ý nghĩa thống kê (p = 0,047). Nguy cơ chết của bệnh nhân có điểm ECOG xấu cao hơn so với bệnh nhân có điểm ECOG tốt vì HR = 1,80 > 1 nhưng không có ý nghĩa thống kê (p = 0,355). Các đường ngang biều hiện khoảng tin cậy 95%. Đường ngang nào có chứa giá trị 1 thì sự khác biệt không có ý nghĩa thống kê. Tham khảo chính: https://www.datacamp.com/community/tutorials/survival-analysis-R

11