Phân Tích Sinh Tồn và Mô Hình Cox: Hướng Dẫn Toàn Diện cho Người Mới Bắt Đầu

Bài viết này là cẩm nang khởi đầu cho bất kỳ ai muốn nắm vững phân tích sinh tồn, từ đường cong Kaplan-Meier trực quan đến mô hình hồi quy Cox mạnh mẽ. Chúng ta sẽ khám phá cách xử lý dữ liệu kiểm duyệt, diễn giải tỷ lệ rủi ro và kiểm tra các giả định quan trọng của mô hình, kèm theo mã nguồn Python chạy được ngay. Đây là kỹ năng không thể thiếu cho các nhà phân tích dữ liệu trong lĩnh vực y tế, tài chính, công nghệ và nhiều ngành khác.

Câu Hỏi Mà Hồi Quy Thông Thường Không Thể Trả Lời

Hãy xem xét một bộ dữ liệu về 432 người được thả khỏi tù và được theo dõi trong một năm. Họ được theo dõi để phát hiện một sự kiện duy nhất: bị bắt lại. Một số bị bắt lại ở tuần thứ 8, số khác ở tuần thứ 30. Phần lớn (318 người) kết thúc nghiên cứu mà không bao giờ bị bắt lại.

Câu hỏi đơn giản đặt ra là: Mất bao lâu để một người tái phạm?

Chúng ta không thể chỉ lấy trung bình thời gian bị bắt. Ba phần tư số người không bao giờ bị bắt, vì vậy họ không có thời gian bị bắt để tính trung bình. Chúng ta cũng không thể loại bỏ những người này. Bởi vì nếu chỉ phân tích những người thất bại, chúng ta sẽ kết luận rằng tất cả mọi người cuối cùng đều tái phạm, điều này vừa sai vừa bi quan. Chúng ta cũng không thể mã hóa những người sống sót như bị bắt sau tuần 52, vì họ đã không bị bắt trong thời gian đó. Họ có thể bị bắt ở một tuần nào đó sau tuần thứ 52 mà chúng ta không bao giờ thấy, hoặc có thể không bao giờ bị bắt.

Tình huống cuối cùng—khi chúng ta biết ai đó sống sót ít nhất đến một thời điểm, nhưng không biết tổng thời gian—được gọi là kiểm duyệt (censoring). Đây là lý do toàn bộ phân tích sinh tồn tồn tại như một lĩnh vực riêng. Hồi quy tuyến tính thông thường không thể đại diện cho câu trả lời "ít nhất là 52." Nó muốn một con số cụ thể. Phân tích sinh tồn là tập hợp các công cụ được xây dựng cho trường hợp mà đối với một phần lớn dữ liệu, đồng hồ vẫn đang chạy khi chúng ta ngừng theo dõi.

Bài viết này tập trung vào vấn đề cốt lõi đó. Chúng ta sẽ bắt đầu với ba khái niệm chính cần hiểu trong Phân tích Sinh tồn. Sau đó, chúng ta sẽ tiến hành ba hoạt động sau:

1. Ước tính đường cong sống sót trực tiếp từ dữ liệu với Kaplan-Meier.
2. Dành phần lớn thời gian để tìm hiểu mô hình hồi quy Cox về tỷ lệ rủi ro.
3. Cuối cùng, chạy một mô hình trong Python và đọc kết quả dưới dạng tỷ lệ rủi ro.

Ba Khái Niệm Cốt Lõi trong Phân Tích Sinh Tồn

Hầu hết mọi thứ trong phân tích sinh tồn được xây dựng từ ba đối tượng sau:

1. Sự kiện và thời gian: Hãy chọn một sự kiện được xác định rõ ràng. Đó có thể là cái chết, bị bắt lại, hỏng hóc máy móc, hủy đăng ký hoặc vỡ nợ khoản vay. Đối với mỗi đối tượng, chúng ta ghi lại hai điều: thời gian quan sát (duration) và chỉ báo sự kiện (event indicator), cho biết liệu quan sát kết thúc vì sự kiện xảy ra hay vì chúng ta ngừng theo dõi (mã nhị phân 1 hoặc 0). Hai cột này là thông tin quan trọng nhất cho mọi mô hình dưới đây.

2. Hàm sống sót, S(t): Đây là xác suất một đối tượng vượt qua thời điểm t mà không có sự kiện xảy ra. Nó bắt đầu ở mức 1 (tại thời điểm 0) và giảm dần về 0 theo thời gian. Ví dụ, câu "Xác suất sống sót sau 12 tháng là 0.7" có nghĩa là 70% đối tượng dự kiến vẫn không có sự kiện sau một năm.

3. Hàm rủi ro (hazard), h(t): Rủi ro là tốc độ tức thời của sự kiện tại thời điểm t, với điều kiện là đối tượng đã sống sót đến t. Nói một cách dễ hiểu: "Trong số những người đã đến được thời điểm này, tỷ lệ thất bại ngay lúc này là bao nhiêu?"

Sự khác biệt giữa S(t) và h(t) rất quan trọng. Sống sót là tích lũy; rủi ro là nhất thời. Xác suất tổng thể chúng ta còn sống ở tuổi 80 là thấp (hàm sống sót nhỏ), nhưng rủi ro tại chính xác thời điểm chúng ta bước sang tuổi 80 lại là một đại lượng hoàn toàn khác. Sống sót giống như lượng nước còn lại trong bể. Rủi ro, mặt khác, là tốc độ nước đang chảy ra ngay lúc này. Chúng được liên kết với nhau bằng một chút giải tích. Rủi ro là mật độ thất bại chia cho tỷ lệ sống sót. Do đó, tích phân của hàm rủi ro sẽ cho chúng ta hàm sống sót.

Bây giờ, bạn có thể tự hỏi, tại sao chúng ta cần tập trung vào rủi ro thay vì mô hình hóa trực tiếp sự sống sót? Bởi vì rủi ro là nơi tự nhiên để gắn các biến đồng hành (covariates). Thật dễ dàng để nói "Hỗ trợ tài chính nhân tỷ lệ tái phạm của chúng tôi với 0.68 tại mọi thời điểm." Câu nói này là một phát biểu về rủi ro, và đó chính xác là những gì mô hình Cox sẽ cung cấp cho chúng ta.

Minh họa khái niệm cơ bản về phân tích sinh tồn

Kaplan-Meier: Đường Cong Sống Sót Không Cần Giả Định

Trước khi mô hình hóa bất cứ điều gì, chúng ta có thể ước tính S(t) trực tiếp từ dữ liệu. Ước lượng Kaplan-Meier (Kaplan & Meier, 1958, một trong những bài báo được trích dẫn nhiều nhất trong Thống kê) làm điều này mà không cần giả định bất kỳ hình dạng cụ thể nào của đường cong.

Ý tưởng rất đơn giản. Chúng ta tiến về phía trước theo thời gian. Tại mỗi thời điểm một sự kiện thực sự xảy ra, chúng ta xem có bao nhiêu người vẫn còn trong nhóm rủi ro ngay trước đó và bao nhiêu người thất bại. Nhân tất cả các "tỷ lệ sống sót ngay lập tức" này lại với nhau khi chúng ta tiến tới. Những người bị kiểm duyệt đóng góp vào số lượng rủi ro cho đến khi họ rời đi, sau đó âm thầm rời khỏi nhóm mà không bao giờ gây ra một bước giảm trên đường cong.

Bây giờ, hãy áp dụng điều này vào bộ dữ liệu tái phạm. Bộ dữ liệu cụ thể được sử dụng trong bài viết này là bộ dữ liệu tái phạm Rossi, được tích hợp trong thư viện lifelines của Python. Nó đến từ một thí nghiệm ngẫu nhiên năm 1980 của Rossi, Berk và Lenihan. Nghiên cứu đã theo dõi 432 tù nhân được thả trong một năm và ghi lại việc họ có bị bắt lại hay không, cùng với các biến đồng hành như chủng tộc, trình độ học vấn, v.v. Chúng ta sẽ chia mọi người thành hai nhóm: những người nhận hỗ trợ tài chính sau khi ra tù và những người không nhận. Trong nghiên cứu ban đầu, khoản hỗ trợ này được chỉ định ngẫu nhiên, vì vậy việc so sánh là công bằng.

from lifelines import KaplanMeierFitter
from lifelines.datasets import load_rossi
import matplotlib.pyplot as plt

df = load_rossi()          # 432 hàng: 'week', 'arrest' (1=sự kiện), + các biến đồng hành

kmf = KaplanMeierFitter()
for value, label in [(0, "Không hỗ trợ tài chính"), (1, "Có hỗ trợ tài chính")]:
    g = df[df.fin == value]
    kmf.fit(g["week"], g["arrest"], label=label)
    kmf.plot_survival_function()

plt.ylabel("S(t): xác suất không bị bắt lại")
plt.show()

Hình 1. Đường cong Kaplan-Meier theo nhóm hỗ trợ tài chính. Mỗi bước giảm là một lần bị bắt lại; các dải bóng mờ là khoảng tin cậy 95%. Nhóm được hỗ trợ (màu xanh) luôn cao hơn trong suốt thời gian. Đến tuần 52, khoảng 22% nhóm được hỗ trợ đã bị bắt lại so với 31% nhóm không được hỗ trợ.

Để kiểm tra xem khoảng cách này có phải chỉ là nhiễu hay không, công cụ tiêu chuẩn là kiểm định log-rank. Nó so sánh số lượng sự kiện quan sát được trong mỗi nhóm với những gì chúng ta mong đợi nếu hai đường cong thực sự giống nhau.

from lifelines.statistics import logrank_test

results = logrank_test(
    df[df["fin"] == 1]["week"],
    df[df["fin"] == 0]["week"],
    event_observed_A=df[df["fin"] == 1]["arrest"],
    event_observed_B=df[df["fin"] == 0]["arrest"],
)
results.p_value

Trên dữ liệu này, nó trả về p ≈ 0.05, ở mức ranh giới. Đây là một lời nhắc nhở tốt rằng Kaplan-Meier cộng với log-rank là mô tả một biến tại một thời điểm. Nó không thể điều chỉnh cho tuổi tác, tiền án, hoặc bất cứ thứ gì khác. Ngay khi bạn muốn kiểm soát các biến đồng hành, chúng ta cần một mô hình Cox.

Mô Hình Cox: Hồi Quy Trên Hàm Rủi Ro

Khi chúng ta muốn một thứ gì đó giống như hồi quy—ví dụ, đưa các biến đồng hành vào, nhận ra tác động của chúng—nhưng cho hàm rủi ro, cách tiếp cận ngây thơ là viết ra một công thức đầy đủ cho h(t) và ước tính mọi thứ. Nhưng h(t) có một hình dạng theo thời gian, và chúng ta thường không biết hình dạng đó là gì và cũng không thực sự muốn cam kết với một hình dạng cụ thể. Sáng kiến của Cox là chúng ta có thể ước tính tác động của các biến đồng hành mà không cần xác định hình dạng cơ bản. Mô hình được xác định như sau:

h(t|x) = h₀(t) * exp(βᵀx)

Hãy đọc nó như hai phần nhân với nhau:

• h₀(t), rủi ro cơ bản: rủi ro theo thời gian cho một đối tượng giả định với tất cả các biến đồng hành bằng 0. Hình dạng hoàn toàn không được xác định. Nó có thể có bất kỳ hình dạng nào. Đây là phần "phi tham số".
• exp(βᵀx), hiệu ứng biến đồng hành: một số duy nhất nhân toàn bộ rủi ro cơ bản lên hoặc xuống tùy thuộc vào đặc điểm cá nhân. Đây là phần "tham số".

Bởi vì mô hình là một nửa phi tham số và một nửa tham số, nó được gọi là mô hình bán tham số (semiparametric).

Đây là phần làm cho mô hình Cox rất hữu ích. Chúng ta có thể lấy hai đối tượng và lập tỷ lệ rủi ro của chúng. Rủi ro cơ bản h₀(t) là giống nhau cho cả hai, vì vậy khi chúng ta lập tỷ lệ, rủi ro cơ bản xuất hiện ở tử số và mẫu số và triệt tiêu lẫn nhau:

HR = h(t|x₁) / h(t|x₂) = exp(βᵀx₁) / exp(βᵀx₂)

Vế phải không có t. Rủi ro cơ bản không xác định và thay đổi theo thời gian đã biến mất. Do đó, tỷ lệ rủi ro là như nhau ở tuần 1, tuần 20 và tuần 52. Chúng ta không bao giờ phải xác định hình dạng của rủi ro cơ bản. Sự triệt tiêu đó là toàn bộ phép màu.

Cox sau đó đã biến điều này thành một phương pháp gọi là khả năng một phần (partial likelihood). Tại mỗi thời điểm sự kiện, mô hình đặt câu hỏi:

"Trong số những người vẫn còn rủi ro ngay bây giờ, khả năng người cụ thể này thất bại cao hơn bao nhiêu so với những người khác?"

Những người bị kiểm duyệt phù hợp một cách tự nhiên trong thiết lập này. Họ ở trong nhóm rủi ro cho đến khi họ rời khỏi nghiên cứu, sau đó âm thầm rời đi. Họ không bao giờ cần có sự kiện. Họ vẫn đóng góp thông tin bằng cách cho chúng ta biết họ không có sự kiện cho đến thời điểm đó.

Có hai điều cần nhắc đến trước khi chạy mô hình:

• Liên kết (Ties): Khả năng một phần giả định chúng ta có thể sắp xếp thứ tự các sự kiện. Khi hai sự kiện xảy ra trong cùng một tuần, phần mềm áp dụng một sự hiệu chỉnh (thường là phương pháp Efron, là mặc định hiện đại và chính xác hơn phương pháp Breslow).
• exp(β) là tỷ lệ rủi ro (hazard ratio). Hệ số β = 0 có nghĩa là exp(β) = 1, tức là không có tác động. β > 0 làm tăng rủi ro. Chúng ta gần như luôn báo cáo phiên bản được mũ hóa này.

Chạy Mô Hình Cox trong Python và Diễn Giải Tỷ Lệ Rủi Ro

Thư viện lifelines làm cho việc chạy mô hình Cox trở nên đơn giản.

from lifelines import CoxPHFitter

cph = CoxPHFitter()
cph.fit(df, duration_col="week", event_col="arrest")

cph.print_summary()

Các tỷ lệ rủi ro được ước tính được trình bày dưới đây:

Biến đồng hành	Tỷ lệ rủi ro exp(β)	95% CI	p-value
Hỗ trợ tài chính	0.68	0.47 – 1.00	0.047
Tuổi tác (mỗi năm)	0.94	0.90 – 0.99	0.009
Tiền án (mỗi lần)	1.10	1.04 – 1.16	0.001
Chủng tộc	1.37	0.75 – 2.50	0.31
Kinh nghiệm làm việc	0.86	0.57 – 1.31	0.48
Đã kết hôn	0.65	0.31 – 1.37	0.26
Đang được ân xá	0.92	0.63 – 1.35	0.67

Tóm lại:

• Hỗ trợ tài chính, HR 0.68. Nhận hỗ trợ có liên quan đến việc giảm 32% rủi ro bị bắt lại tại bất kỳ thời điểm nào. Đây là hiệu quả điều trị mà thí nghiệm ban đầu được xây dựng để tìm ra.
• Tuổi tác, HR 0.94. Mỗi năm tuổi tác tăng thêm làm giảm rủi ro khoảng 6%. Người được thả lớn tuổi hơn tái phạm chậm hơn.
• Tiền án, HR 1.10. Mỗi tiền án làm tăng rủi ro khoảng 10%, và nó nhân lên: năm tiền án là khoảng 1.10⁵ ≈ 1.6 lần rủi ro.

Tỷ lệ rủi ro không phải là sự khác biệt về xác suất sống sót. Nó là một số nhân cho tỷ lệ tức thời, được giả định là không đổi trong suốt thời gian theo dõi. Điều này đưa chúng ta đến câu hỏi tiếp theo: điều đó có thực sự không đổi không? Phần còn lại của bài viết sẽ giải thích điều này.

Kết quả mô hình Cox trong Python

Giả Định Trong Chính Cái Tên: Tỷ Lệ Rủi Ro Tương Xứng

Mô hình được gọi là tỷ lệ rủi ro tương xứng (proportional hazards) là có lý do. Chúng ta hãy xem lại tỷ lệ đã làm cho mọi thứ hoạt động:

HR = exp(βᵀx₁) / exp(βᵀx₂)

Vế phải không có 't'. Đó là giả định tỷ lệ rủi ro tương xứng (PH) được phát biểu chính xác: tỷ lệ rủi ro giữa hai đối tượng bất kỳ là không đổi theo thời gian. Hỗ trợ tài chính làm giảm 32% rủi ro ở tuần 2 và chính xác 32% ở tuần 50. Rủi ro của hai nhóm tăng và giảm cùng nhau. Một nhóm luôn là một bội số cố định của nhóm kia. Đường cong sống sót của họ không bao giờ có thể cắt nhau.

Đôi khi điều đó là đúng. Thường thì không. Một phương pháp điều trị có thể giúp ích rất nhiều lúc đầu và sau đó giảm tác dụng. Một yếu tố rủi ro chỉ có thể ảnh hưởng về lâu dài. Khi tỷ lệ rủi ro thực sự thay đổi theo thời gian, một mô hình Cox đơn giản sẽ tính trung bình nó thành một con số duy nhất sai ở cả hai đầu. Vì vậy, chúng ta phải kiểm tra.

Sự kiểm tra mà hầu hết mọi người nên sử dụng dựa trên phần dư Schoenfeld (Schoenfeld, 1982; và phiên bản được chia tỷ lệ với kiểm định chính thức của Grambsch & Therneau, 1994). Tại mỗi thời điểm sự kiện, phần dư Schoenfeld đo lường khoảng cách giữa giá trị biến đồng hành của người thực sự thất bại và giá trị trung bình của biến đồng hành trong số những người có rủi ro. Nếu giả định PH đúng, các phần dư này không nên cho thấy xu hướng theo thời gian. Nếu chúng ta có thể thấy một độ dốc, hiệu ứng đang thay đổi theo thời gian và tính tương xứng bị phá vỡ.

cph.check_assumptions(df, p_value_threshold=0.05, show_plots=True)

Trên dữ liệu của chúng ta, kiểm định chỉ ra hai vi phạm rõ ràng. Các biến 'tuổi tác' (age) và 'kinh nghiệm làm việc' (wexp) đã thất bại trong kiểm định phi tương xứng với p-values lần lượt là 0.0007 và 0.0063.

Ví dụ này cho thấy kết quả thực sự trên một bộ dữ liệu nổi tiếng, có thật. Nó chính xác là loại điều sẽ không được nhìn thấy nếu chúng ta chạy mô hình, đọc p-values và bỏ đi. Hệ số cho tuổi tác (HR 0.94) không sai nhưng chưa đầy đủ. Về cơ bản, nó thu gọn một hiệu ứng thực sự thay đổi theo thời gian thành một con số.

Phải Làm Gì Khi Tính Tương Xứng Bị Phá Vỡ

Việc tìm ra một vi phạm không phải là ngõ cụt. Nó thường là phát hiện thú vị nhất trong phân tích, và có ba chiến lược khắc phục tiêu chuẩn.

1. Phân tầng (Stratify): Nếu một biến vi phạm PH nhưng chúng ta chỉ cần điều chỉnh cho nó (không ước tính hiệu ứng của nó), thì chúng ta có thể đặt nó vào các tầng (strata). Phân tầng phù hợp với một rủi ro cơ bản riêng cho từng mức của biến đó và không bao giờ ép buộc hiệu ứng của nó phải tương xứng.

cph_strat = CoxPHFitter()
cph_strat.fit(
    df,
    duration_col="week",
    event_col="arrest",
    strata=["wexp"]
)
cph_strat.print_summary()

Sau khi phân tầng theo kinh nghiệm làm việc để khắc phục vi phạm tỷ lệ rủi ro tương xứng, các kết quả chính vẫn được giữ nguyên, tức là hỗ trợ tài chính (HR 0.68), tuổi tác lớn hơn (HR 0.94 mỗi năm) và tiền án (HR 1.09 mỗi lần) vẫn có ý nghĩa thống kê, trong khi chủng tộc, tình trạng hôn nhân và tình trạng ân xá thì không. Sự nhất quán (concordance) của mô hình ở mức 0.61 cho thấy khả năng khiêm tốn trong việc xếp hạng ai sẽ bị bắt lại sớm hơn.

2. Để hiệu ứng thay đổi theo thời gian. Nếu chúng ta thực sự quan tâm đến cách hiệu ứng thay đổi, chúng ta cần thêm một tương tác giữa biến đồng hành và một hàm của thời gian. Điều này cho phép tỷ lệ rủi ro của nó tăng hoặc giảm trong suốt quá trình theo dõi. Trong lifelines, đây là CoxTimeVaryingFitter với dữ liệu được chia tại các thời điểm sự kiện. Về cơ bản, chúng ta đang phù hợp với β(t) thay vì một β duy nhất.

3. Kiểm tra dạng hàm số trước tiên. Kiểm định Schoenfeld nhạy cảm với một biến đồng hành được xác định sai. Nếu hiệu ứng thực sự của tuổi tác là phi tuyến tính (rủi ro giảm nhanh rồi ổn định), việc đưa nó vào như một đường thẳng có thể làm cho kiểm định PH thất bại ngay cả khi các rủi ro thực sự tương xứng. Trước khi sử dụng các mô hình thời gian biến thiên, hãy thử thêm một số hạng bình phương cho biến gây ra vi phạm.

Những Điểm Chính Cần Ghi Nhớ

• Kiểm duyệt là một đặc điểm, không phải dữ liệu bị thiếu. Toàn bộ điểm của các phương pháp này là sử dụng chính xác thông tin một phần trong các đối tượng bị kiểm duyệt. Nếu bạn thấy mình đang loại bỏ những người không có sự kiện, bạn về cơ bản đang đưa ra chính xác sự sai lệch mà lĩnh vực này được phát minh để tránh.
• Kaplan-Meier để quan sát, Cox để điều chỉnh. Vẽ biểu đồ các đường cong KM trước. Chúng ít giả định và xây dựng trực giác về hình dạng của sự sống sót. Chuyển sang Cox khi bạn cần kiểm soát các biến đồng hành. Báo cáo các tỷ lệ rủi ro với khoảng tin cậy, không chỉ p-values.
• Một tỷ lệ rủi ro nhân rủi ro của bạn, nó không trừ đi. HR 0.68 có nghĩa là hỗ trợ cắt giảm tỷ lệ bị bắt lại xuống còn khoảng hai phần ba tại mọi thời điểm. Nó không cho bạn biết có bao nhiêu người tránh bị bắt, hoặc họ ở miễn phí lâu hơn bao lâu. Ngoài ra, nó âm thầm giả định rằng hai phần ba đó giữ nguyên trong suốt thời gian, đó là điều bạn phải kiểm tra.
• Luôn kiểm tra tỷ lệ rủi ro tương xứng. Nó chỉ là hai dòng mã và là điều duy nhất tách biệt một phân tích sinh tồn có thể bảo vệ được khỏi một phân tích dễ vỡ. Khi nó thất bại, bạn có thể phân tầng, mô hình hóa sự biến thiên theo thời gian hoặc sửa dạng hàm số.
• Chỉ số nhất quán (Concordance), không phải R², để đánh giá độ phù hợp. Chỉ số nhất quán đo lường mức độ tốt của mô hình trong việc xếp hạng ai thất bại sớm hơn.

Phân tích sinh tồn thực chất chỉ là hồi quy trung thực về những gì nó chưa biết. Nắm vững tỷ lệ rủi ro, tôn trọng giả định trong chính tên gọi của nó, và bạn sẽ có được những câu trả lời đúng đắn ngay cả sau khi nghiên cứu kết thúc.

Đường cong Kaplan-Meier cho các mô hình khác nhau

Tài Liệu Tham Khảo

Nền tảng

• Cox, D. R. (1972). Regression Models and Life-Tables. Journal of the Royal Statistical Society, Series B, 34(2), 187–220.
• Kaplan, E. L., & Meier, P. (1958). Nonparametric Estimation from Incomplete Observations. Journal of the American Statistical Association, 53(282), 457–481.

Chẩn đoán và mô hình Cox chuyên sâu

• Therneau, T. M., & Grambsch, P. M. (2000). Modeling Survival Data: Extending the Cox Model. Springer.
• Grambsch, P. M., & Therneau, T. M. (1994). Proportional Hazards Tests and Diagnostics Based on Weighted Residuals. Biometrika, 81(3), 515–526.
• Schoenfeld, D. (1982). Partial Residuals for the Proportional Hazards Regression Model. Biometrika, 69(1), 239–241.