ĐỀ TÀI: PACKAGE SURVIVAL VÀ PHÂN TÍCH SỐNG CÒN TRÊN R
Nguyễn Khánh An
Lê Phạm Thị Phương Tiền
2023-07-01
1 PHẦN MỞ ĐẦU
1.1 Lời nói đầu
Survival – Phân tích sự tồn tại
_ Phân tích sống còn_ Phân tích sống sót_ Phân tích sự
kiện
Phân tích sống sót (PTSS) hoặc Phân tích sự kiện khi nhà nghiên cứu muốn tìm hiểu ảnh hưởng đến các biến kết cục (biến phụ thuộc) mang tính thời gian.
Khi so sánh 2 phương pháp điều trị cho các bệnh có tần số tử vong cao như bệnh AIDS, các bệnh ung thư. Nếu mô hình phân tích như phân tích hồi qui logistic, chỉ để ý đến biến kết cục (sống/chết hoặc khỏi bệnh/không khỏi bệnh) mà không quan tâm đến yếu tố thời gian thì đôi không tìm thấy sự khác biệt giữa 2 phương pháp điều trị vì tỉ lệ tử vong gần như nhau, nhưng thời gian dẫn đến tử vong ở 2 nhóm có thể khác nhau. Vì vậy chúng ta phải sử dụng mô hình PTSS thì mới thấy sự khác biệt này. Như vậy mô hình nghiên cứu mô tả kết cục là biến nhị phân (sống/chết) tuy quan trọng nhưng không chính xác.
Mô hình chính để thể hiện mối liên hệ giữa thời gian dẫn đến bệnh
(hay không bệnh) và các yếu tố nguy cơ (risk factors)
là mô hình có tên là “survival analysis – Phân tích sự sông
sót”. Cụm từ survival analysis xuất phát từ nghiên
cứu trong bảo hiểm, và giới nghiên cứu y khoa từ đó dùng cụm từ cho bộ
môn của mình. Nhưng như nói trên, sống/chết không phải là biến
duy nhất, vì trong thực tế chúng ta cũng có những biến cố có
bệnh hay không bệnh, xảy ra hay không xảy ra. Ngoài ra, trong các bộ môn
kĩ thuật, người ta dùng một cụm từ khác “reliability analysis –
Phân tích độ tin cậy”, để chỉ cho khái niệm
survival analysis. Tuy nhiên trong đề tài này chúng tôi sẽ
dùng cụm từ Phân tích sống còn.
1.2 Lý do chọn đề tài
Mô hình hóa dữ liệu thời gian đến sự kiện là một
chủ đề quan trọng với nhiều ứng dụng trong các lĩnh vực khác nhau. Tập
hợp các phương pháp để phân tích dữ liệu đó được gọi là phân tích sự tồn
tại, phân tích lịch sử sự kiện hoặc phân tích thời lượng. Phân tích sống
còn được áp dụng rộng rãi vì định nghĩa về một ‘sự kiện’ có thể
rất đa dạng và các ví dụ bao gồm chết và sống, thành công
và phá sản, tiếp tục hoặc hủy bỏ. Do đó, các lĩnh vực ứng
dụng bao gồm từ y học và xã hội học đến tiếp thị và kinh tế. Trong bài
tiểu luận này, chúng ta xem xét các vấn đề cơ bản về lý thuyết của phân
tích sống còn bao gồm các công cụ ước tính cho các chức năng sinh tồn và
rủi ro và cách thực hành trên R dựa trên
packagesurvival.
Trong thời đại ngày nay - Công nghiệp 4.0 là thời đại cho sự phát triển kinh tế xã hội và công nghệ kĩ thuật, việc ứng dụng của phân tích sống sót trong tài chính doanh nghiệp trở thành một vấn đề có ý nghĩa thời sự, đặc biệt trong bối cảnh sự phát triển đa ngành nghề một cách sâu rộng và nhanh chóng như hiện nay.
Thấy được tầm quan trọng cũng như tiềm năng khi phân tích sự tồn
tại của package survival trong thời điểm hiện tại, chúng
tôi mong muốn có thể kết hợp các ngôn ngữ lập trình để giải quyết vấn đề
trong kinh tế xã hội hiện và mô hình này sẽ giúp ích cho việc phân tích
các nguy cơ, sự sống sót có thể xảy ra từ đó có thể dễ dàng giải quyết
được các vấn đề quan trọng.
1.3 Mục đích đề tài
Một ưu điểm của Phân tích sống còn là xử lý được các trường hợp đối tượng nghiên cứu bỏ cuộc giữa chừng (như mất dấu theo dõi, ngưng điều trị do tác dụng phụ của thuốc hoặc tử vong do bệnh lý khác…).
Trong mô hình phân tích này các đối tượng còn sống kể cả đối tượng
bỏ cuộc được gọi là censored hoặc sự kiện chưa xảy
ra. Các đối tượng tử vong hoặc hết sống được gọi là
events hoặc sự kiện đã kết thúc. Mục tiêu suy luận
đối với phân tích sống còn là khoảng thời gian giữa thời điểm bắt đầu và
thời điểm sự kiện xảy ra. Trong nghiên cứu y học hiện nay, phân tích
sống còn được sử dụng rộng rãi trong các nghiên cứu lâm sàng để đánh giá
hiệu quả của một phương phương điều trị hoặc để đánh giá tình trạng sống
còn của một số các biện pháp điều trị ung thư. Trong kinh tế xã hội,
phân tích sống còn được dùng để phân tích rủi ro/ nguy
cơ để đánh giá sự tồn tại, của một doanh nghiệp nào đó.
1.4 Đối tượng, phạm vi nghiên cứu
Đối tượng nghiên cứu
Phân tích sống còn tập trung mô tả cho một cá thể hay một nhóm cá
thể nhất định. Một điểm xác định của một sự kiện được gọi là
failure (như là xuất hiện bệnh, chữa khỏi bệnh, tử
vong, tái phát sau khi đáp ứng với điều trị…) mà xảy ra sau một
khoảng thời gian được gọi là failure time (thời gian
dẫn đến sự kiện) hoặc (follow-up time (thời gian
theo dõi) trong nghiên cứu thuần tập/nghiên cứu dựa vào dân số)
trong suốt thời gian các cá thể được quan sát. Để xác định thời gian dẫn
đến sự kiện, chúng ta cần xác định thời điểm bắt đầu (có thể là ngày
nhận vào, ngày chẩn đoán…).
Phạm vi nghiên cứu
Trong bài báo này, chúng tôi xem xét các vấn đề cơ bản về lý thuyết của phân tích tỷ lệ sống sót bao gồm các công cụ ước tính cho các chức năng sinh tồn và rủi ro. Chúng tôi thảo luận chi tiết về Mô hình mối nguy theo tỷ lệ Cox và cả các phương pháp để kiểm tra giả định về mối nguy theo tỷ lệ (PH). Hơn nữa, chúng tôi thảo luận về các mô hình Cox phân tầng cho các trường hợp khi giả định PH không đúng.
Phần thảo luận của chúng tôi được bổ sung bằng một ví dụ hoạt động sử dụng ngôn ngữ lập trình thống kê R để cho phép ứng dụng thực tế của phương pháp này. Phép ước tính thường được dùng để Phân tích sống còn được gọi là ước tính KaplanMeier. Phép ước tính này giúp ta tính được xác suất sống sót tích lũy tại các mốc thời gian khác nhau . Nếu muốn so sánh sự khác biệt giữa 2 nhóm điều trị, dùng kiểm định log-rank bằng cách so sánh 2 hàm xác suất tích lũy của 2 nhóm.
2 TỔNG QUAN VỀ PACKAGE SURVIVAL
2.1 Thông tin package
Package survival chứa các quy trình phân tích sự sống
sót cốt lõi, bao gồm định nghĩa về các đối tượng Surv, đường
cong Kaplan-Meier và Aalen-Johansen (đa trạng thái),
mô hình Cox và mô hình thời gian thất bại tăng tốc theo tham
số.
Phiên bản mới nhất: 3.5-5
Phiên bản R hỗ trợ: R(>3.5.0)
Nhóm tác giá: Terry M Therneau, Thomas Lumley, Atkinson Elizabeth, Crowson Cynthia
Người phụ trách chính: Terry M Therneau
Tài liệu hướng dẫn: https://cran.r-project.org/web/packages/survival/survival.pdf
2.2 Giới thiệu package survival
Phân tích sống còn được sử dụng phổ biến trong sinh học, y học, kỹ thuật, tiếp thị, khoa học xã hội hoặc khoa học.
Có hai phương pháp đóng góp quan trọng cho sự phát triển của lĩnh vực này. Đầu tiên là phương pháp của Kaplan và Meier, người đã giới thiệu một công cụ “Ước tính xác suất sống sót”. Thứ hai là từ Cox, người đã giới thiệu một mô hình được gọi là “Mô hình rủi ro theo tỷ lệ Cox” (CPHM), là một mô hình hồi quy. Thật đặc biệt khi cả hai mô hình này đều được sử dụng nhiều cho đến nay và luôn là lực chọn của mọi nhà khoa học khi phân tích về dữ liệu sống còn (survival data).
Trong những năm qua, nhiều báo cáo đã được viết về khảo sát phân tích tỷ lệ sống còn. Tuy nhiên, do sự phức tạp của các phương pháp, đặc biệt là đối với dữ liệu đa biến, các cuộc thảo luận dễ dẫn đến nhầm lẫn. Ngoài ra, ở cấp độ giới thiệu có sẵn các đánh giá, đây hoàn toàn là lý thuyết hoặc sử dụng các ngôn ngữ lập trình (như stata hoặc SAS) được sử dụng chủ yếu trong các thống kê sinh học.
Ngược lại, đánh giá của một số người đã kết hợp trình bày lý thuyết với thực tế sử dụng ngôn ngữ lập trình thống kê R. R là ngôn ngữ được sử dụng rộng rãi cho các vấn đề chung trong khoa học dữ liệu vì nó kết hợp các tính năng từ các mô hình lập trình khác nhau. Chúng ta hướng tới một cách trình bày toàn diện và đủ rộng để bao gồm tất cả các chủ đề cần thiết cho phân tích sống còn đa biến nhưng đồng thời cũng dễ hiểu.
Để thực hiện tất cả điều này, cần bổ sung cho việc trình bày các phương pháp với đầy đủ thông tin cơ bản. Giải thích chi tiết việc kiểm chứng và xử lý dữ liệu về các sự kiện thời gian, bởi vì về cơ bản tất cả các phương pháp đều sử dụng kiểm chứng để rút ra các ước tính hiệu quả. Cuối cùng, chúng ta thêm một ví dụ đã hoạt động cho thấy cách thực hiện phân tích khả năng sống sót thực tế với R.
Trong phần tiếp theo, chúng ta cung cấp một dữ liệu để phân tích sự sống còn và ý nghĩa của việc phân tích chúng. Tiếp theo, chúng ta mô tả các đặc điểm chung của các hàm sinh tồn, ước lượng cho các hàm sinh tồn và so sánh hai đường cong sinh tồn. Sau đó, chúng ta giới thiệu Mô hình nguy cơ theo tỷ lệ Cox. Nếu giả định rủi ro theo tỷ lệ không đúng, người ta cần sử dụng mô hình Cox phân tầng. Cuối cùng, chúng ta trình bày một phần về phân tích tỷ lệ sống thực tế bằng R. Bài báo kết thúc với phần tóm tắt ngắn gọn và kết luận.
2.3 Các hàm cơ bản trong R
Dưới đây là một số hàm cơ bản trong package
survival:
Surv(time, event): được sử dụng để tạo ra đối tượngSurvtrong R. Đối tượng này được sử dụng để định nghĩa và lưu trữ dữ liệu sống sót và thời gian theo dõi tương ứng.Survfit(formula, data): được sử dụng để tính toán tỷ lệ sống sót hoặc phân tích thời gian đến sự kiện cho các nhóm dữ liệu khác nhau.Coxph(Surv (time, event) ~ I, data): mô hình cox tiêu chuẩn, mô hình Cox được sử dụng để xác định yếu tố ảnh hưởng đến thời gian sống sót hoặc thời gian đến sự kiện.Survreg(formula, data): Hàm này được sử dụng để phân tích các yếu tố ảnh hưởng đến thời gian sống sót hoặc thời gian đến sự kiện trong mô hình hồi quy.Survexp(): Hàm này được sử dụng để tính toán tỷ lệ sống sót hoặc phân tích thời gian đến sự kiện cho các nhóm dữ liệu khác nhau trong mô hình phân phối mũ.
3 LÝ THUYẾT PHÂN TÍCH SỐNG SÓT - SURVIVAL ANALYSIS
3.1 Phân tích sống còn là gì?
Phân tích sống còn là một nhánh của thống kê để phân tích khoảng thời gian dự kiến cho đến khi một sự kiện xảy ra, chẳng hạn như cái chết ở các sinh vật và sự cố trong các hệ thống máy móc. Chủ đề này được gọi là lý thuyết độ tin cậy hoặc phân tích độ tin cậy trong kỹ thuật, phân tích thời lượng hoặc mô hình thời lượng trong kinh tế học và phân tích lịch sử sự kiện trong xã hội học.
Phân tích sống còn cố gắng trả lời một số câu hỏi, chẳng hạn như tỷ lệ dân số sẽ sống sót qua một thời gian nhất định là bao nhiêu? Trong số những người sống sót, họ sẽ chết hoặc thất bại ở mức nào? Có thể tính đến nhiều nguyên nhân dẫn đến tử vong hoặc hỏng hóc không? Các hoàn cảnh, đặc điểm cụ thể làm tăng hoặc giảm xác suất sống sót như thế nào?
Để trả lời những câu hỏi như vậy, cần phải “xác định trọn đời”. Trong trường hợp tồn tại sinh học, cái chết là rõ ràng, nhưng đối với độ tin cậy cơ học, sự thất bại có thể không được xác định rõ ràng, vì có thể có những hệ thống cơ khí trong đó sự thất bại là một phần, một vấn đề khác nằm ở mức độ thời gian. Lý thuyết được phác thảo dưới đây giả định các sự kiện được xác định rõ ràng vào những thời điểm cụ thể; các trường hợp khác có thể được xử lý tốt hơn bằng các mô hình khác nhau.
Tổng quát hơn, phân tích sự sống còn liên quan đến việc lập mô hình thời gian cho dữ liệu sự kiện; trong bối cảnh này, cái chết hoặc thất bại được coi là một “sự kiện” trong tài liệu phân tích sự sống sót - theo truyền thống, chỉ một sự kiện duy nhất xảy ra đối với mỗi đối tượng, sau đó sinh vật hoặc cơ chế sẽ chết hoặc bị hỏng. Các mô hình sự kiện lặp đi lặp lại làm giảm bớt giả định đó. Nghiên cứu về các sự kiện định kỳ có liên quan đến độ tin cậy của hệ thống trong nhiều lĩnh vực khoa học xã hội và nghiên cứu y học.
Vậy, phân tích sống còn (Survival Analysis) 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ố (time-to-failure/
time-to-event).
Các thuật ngữ cơ bản của phân tích sống còn bao gồm:
• Dữ liệu thời gian và biến cố
• Nhóm kiểm chứng (censoring)
• Hàm sinh tồn và hàm nguy cơ
Dữ liệu thời gian và biến cố
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àn nên SA đôi khi được gọi dưới tên là “Nghiên cứu lâm sàng”, mặc dù số liệu sống còn vẫn có thể có các trong các nghiên cứu khác. 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.
Ví dụ về các dữ liệu thời gian dẫn đến biến cố ( time-to-event)
• Buld: thời gian sống sót của bóng đèn (từ lúc sản xuất
đến khi bị hỏng)
• Deat: tử vong ( thời gian từ lúc bệnh nhân được chuẩn
đoán một loại bệnh nào đó đến khi họ tử vong) Những vấn đề gặp phải khi
sử dụng SA
• Data censoring and truncation: Số liệu bị cắt xén.
• Conditional probability: Xác suất có điều kiện.
Nhóm kiểm chứng
Để thực hiện phân tích sống còn, người ta cần ghi lại thời gian diễn ra sự kiện 𝑡i cho các đối tượng 𝑖∈{1,…,𝑁} của một nhóm. Tuy nhiên, điều này không phải lúc nào cũng khả thi và chúng ta chỉ có một phần thông tin về thời gian diễn ra sự kiện. Trong trường hợp như vậy, người ta nói về kiểm chứng. Cụ thể, một bệnh nhân có thời gian sống sót được kiểm chứng nếu sự kiện chưa xảy ra đối với bệnh nhân này. Điều này có thể xảy ra khi:
• Bệnh nhân bỏ nghiên cứu, ví dụ: ngừng đến phòng khám để tái khám.
• Nghiên cứu các mốc thời gian cố định và sự kiện xảy ra sau thời hạn.
• Bệnh nhân rút khỏi nghiên cứu.
Những ví dụ trên được gọi là right-censoring .Trong
hình bên dưới, chúng ta hình dung ý nghĩa của việc kiểm chứng. Chẳng
hạn, các đối tượng có nhãn ID𝐴 và ID𝐶 trải nghiệm sự kiện trong suốt
thời gian nghiên cứu, được biểu thị bằng một hình tròn màu xanh. Ngược
lại, đối tượng ID𝐵 trải qua sự kiện sau khi kết thúc nghiên cứu, được
biểu thị bằng một vòng tròn màu xanh. Tuy nhiên, điều này không được
quan sát bởi nghiên cứu. Thông tin duy nhất mà chúng ta có thể sử dụng
là vào cuối nghiên cứu, đối tượng ID𝐵 vẫn chưa trải qua sự kiện này.
Do đó, thời gian tồn tại của chủ thể ID𝐵 bị kiểm chứng, được biểu thị bằng dấu X. Điều này có nghĩa là cho đến khi sự kiện kiểm chứng xảy ra (biểu thị bằng dấu X), chủ thể ID𝐵 không trải qua sự kiện này. Ngoài ra, đối với chủ đề ID𝐷 , chúng ta có thời gian tồn tại bị kiểm chứng, tuy nhiên vì một lý do khác trong trường hợp này, nghiên cứu vẫn chưa kết thúc. Một lời giải thích dễ hiểu cho việc kiểm chứng này là đối tượng đã không tham gia các kiểm chứng tiếp theo sau khi sự kiện kiểm chứng xảy ra (biểu thị bằng dấu X). Chính thức, người ta gọi việc kiểm chứng ID chủ thể 𝐵 là cố định, kiểm chứng đúng và kiểm chứng ID chủ thể 𝐷 ngẫu nhiên, kiểm chứng đúng.
Đây là những trường hợp kiểm chứng khác mà người ta có thể phân biệt. Chẳng hạn, người ta nói về kiểm chứng trái nếu sự kiện được quan sát nhưng không phải là điểm bắt đầu của quá trình. Một ví dụ cho điều này là do nhiễm trùng vì thông thường nhiễm trùng được chẩn đoán vào một thời điểm nào đó, nhưng bản thân nhiễm trùng đã xảy ra trước đó tại một điểm khởi đầu không xác định. Sau đây, chúng tôi sẽ giới hạn tập trung vào các chủ đề bị kiểm chứng bên phải. Một bản tóm tắt cho các loại kiểm chứng khác nhau được đưa ra bởi:
• Kiểm chứng loại I: Tất cả các đối tượng bắt đầu và kết thúc nghiên cứu cùng một lúc (thời lượng nghiên cứu cố định). Ví dụ là các thí nghiệm trong phòng thí nghiệm.
• Kiểm chứng loại II: Tất cả các đối tượng bắt đầu nghiên cứu cùng một lúc nhưng nghiên cứu kết thúc thì một số lượng đối tượng cố định được xác định trước đã trải qua sự kiện (thời lượng nghiên cứu linh hoạt). Ví dụ là các thí nghiệm trong phòng thí nghiệm.
• Kiểm chứng loại III: Các đối tượng tham gia nghiên cứu vào các thời điểm khác nhau, nhưng thời lượng của nghiên cứu là cố định. Ví dụ là các thử nghiệm lâm sàng.
3.2 Các hàm được sử dụng trong phân tích sống còn
Hàm sinh tồn
Đường cong tỷ lệ sống còn 𝑆(𝑡) biểu thị tỷ lệ sống sót là một hàm của thời gian (t). Ở đây, 𝑆(𝑡) được gọi là hàm sống còn hay còn gọi là hàm tồn tại hay chỉ đường cong tồn tại. 𝑆(𝑡) là xác suất để biến ngẫu nhiên T lớn hơn một khoảng thời gian (t) xác định, nghĩa là,
𝑆(𝑡)=𝑃𝑟(𝑇>𝑡). (1)
Vì 𝑆(𝑡) được xác định cho một nhóm đối tượng nên 𝑆(𝑡) có thể được hiểu là tỷ lệ đối tượng sống sót cho đến thời gian (t). Từ đó suy ra một công cụ ước lượng cho 𝑆(𝑡) :
(2)
Trong đó N là tổng số đối tượng. Phương trình (1) là ước tính dân số của hàm sinh tồn. Dưới đây, chúng tôi sẽ thảo luận về các ước tính mẫu khác nhau cho điều này có thể được đánh giá bằng số từ dữ liệu. Nói một cách đơn giản, hàm sinh tồn đưa ra xác suất mà một đối tượng (được biểu thị bằng T) sẽ sống sót qua thời gian t.
Hàm sinh tồn có các thuộc tính sau:
Khoảng xác định của thời gian là 𝑡 ∈ [ 0 , ∞ ) .
𝑆(𝑡) không tăng, nghĩa là
𝑆(𝑡1)≥𝑆(𝑡2)với 𝑡1≤𝑡2Tại thời điểm 𝑡=0,
𝑆(𝑡=0)=1, tức là xác suất sống sót qua thời điểm 0 là 1.
Do 𝑆(𝑡) là một xác suất nên tồn tại một mật độ xác suất f với tính chất
(3)
(4)
Hơn nữa, giá trị kỳ vọng của T được xác định bởi
(5)
Đây là trung bình kỳ vọng sống còn của bệnh nhân được biểu thị bằng f.
Sử dụng phương trình (4) và lấy tích phân từng phần, người ta có thể chỉ ra rằng hàm sống sót 𝑆(𝑡) có thể được sử dụng để tính tuổi thọ trung bình,
(6)
Công cụ ước tính cho hàm sinh tồn
- Kaplan–Meier
Công thức ước lượng Kaplan–Meier (KM) của hàm sinh tồn 𝑆𝐾𝑀(𝑡) được xác định bởi
(7)
Ước lượng này đúng với mọi 𝑡>0 và nó chỉ phụ thuộc vào hai biến, 𝑛𝑖 và 𝑑𝑖 là:
𝑛𝑖: số rủi ro tại thời điểm 𝑡𝑖,𝑑𝑖: số sự kiện tại thời điểm 𝑡𝑖
Ở đây, 𝑛𝑖 tương ứng với số đối tượng có mặt tại thời điểm 𝑡𝑖 . Trong hợp đồng, các đối tượng trải qua sự kiện hoặc kiểm duyệt không còn hiện diện. Phần khó khăn của công cụ ước tính này là đối số của tích, chỉ xem xét các sự kiện i xảy ra trước thời điểm t, tức là 𝑡𝑖<𝑡 . Do đó, đường cong sống sót 𝑆𝐾𝑀(𝑡) trong thời gian t xem xét tất cả các sự kiện đã xảy ra trước t.
Điều quan trọng là phải nhận ra rằng, chỉ để đánh giá công cụ ước tính Kaplan–Meier, các sự kiện xảy ra tại {𝑡𝑖} là rất quan trọng. Điều đó có nghĩa là, giữa hai sự kiện, ví dụ: 𝑡𝑖 và 𝑡𝑖+1 , đường cong sống sót không đổi. Điều này cho phép xây dựng lại phương trình (7) đơn giản để viết công cụ ước tính Kaplan–Meier bằng một công thức đệ quy được đưa ra bởi
(8)
Trong Hình 2, chúng tôi hiển thị một ví dụ để đánh giá công cụ ước tính Kaplan–Meier sử dụng Phương trình đệ quy (8). Ví dụ hiển thị bao gồm tổng cộng năm chủ đề và bốn sự kiện.
Hình 2
Hình 2. Ví dụ số cho công cụ ước tính Kaplan–Meier. Đối với mỗi sự kiện, 𝑆𝐾𝑀 được đánh giá đệ quy.
So sánh hai đường cong sinh tồn
Khi một người có nhiều hơn một đường cong sinh tồn, người ta thường quan tâm đến việc so sánh những đường cong này. Trong phần tiếp theo, chúng tôi giả định rằng chúng tôi có hai đường cong sinh tồn tương ứng với hai nhóm đối tượng khác nhau, ví dụ: một nhóm nhận được thuốc, trong khi nhóm kia nhận giả dược.
Về mặt thống kê, việc so sánh có thể được thực hiện bằng một thử nghiệm dựa trên các giả thuyết:
Giả thuyết 1. Không có sự khác biệt về tỷ lệ sống sót giữa (nhóm 1) và (nhóm 2).
Giả thuyết 2. Có sự khác biệt về tỷ lệ sống sót giữa (nhóm 1) và (nhóm 2).
Các thử nghiệm phổ biến nhất để so sánh các đường cong sinh tồn là:
Kiểm tra thứ hạng log,
Kiểm định Wilcoxon (Gehan) (trường hợp đặc biệt của kiểm tra thứ hạng Log có trọng số).
Sự khác biệt giữa cả hai phép thử là phép thử Log-rank có nhiều năng lượng hơn phép thử Wilcoxon trong việc phát hiện sự khác biệt muộn trong các đường cong sinh tồn, trong khi phép thử Wilcoxon có nhiều năng lực hơn phép thử Log-rank để phát hiện sự khác biệt sớm.
Hàm nguy cơ
Hàm nguy cơ (Survival Function) là một khái niệm quan trọng trong phân tích sống còn (Survival Analysis), được sử dụng để mô tả xác suất sống sót của một cá thể trong thời gian tới.
Hàm nguy cơ được định nghĩa là hàm số S(t), trong đó t là thời gian và S(t) là xác suất sống sót của một cá thể sau thời gian t.
Hàm nguy cơ thường được sử dụng để so sánh giữa các nhóm đối tượng khác nhau và định lượng tác động của các yếu tố đến xác suất sống sót của các đối tượng. Nó cũng là một công cụ hữu hiệu để dự đoán xác suất sống sót của một cá thể hoặc nhóm các cá thể trong tương lai.
Hàm nguy cơ không chỉ được sử dụng trong phân tích sống còn mà còn có thể được áp dụng trong nhiều lĩnh vực khác, như sau:
Kinh tế học: được sử dụng để đánh giá rủi ro đối với các quyết định đầu tư hoặc các quyết định kinh doanh khác trong các ngành công nghiệp như bảo hiểm, tài chính, và bất động sản.
Y học: sử dụng để đánh giá tác động của các yếu tố khác nhau đến xác suất mắc bệnh, xác suất phục hồi, hoặc xác suất sống sót của bệnh nhân.
…
Trong lĩnh vực kinh doanh, hàm nguy cơ thường được sử dụng để đánh giá rủi ro và dự đoán xác suất thành công của một dự án hoặc quyết định đầu tư. Việc sử dụng hàm nguy cơ trong kinh doanh có thể giúp các nhà quản lý đánh giá các rủi ro tiềm năng và đưa ra các quyết định thông minh về đầu tư và quản lý rủi ro.
Ví dụ về cách sử dụng hàm nguy cơ trong lĩnh vực kinh doanh:
Đánh giá rủi ro đầu tư
Quản lý rủi ro
Dự đoán doanh thu
Tóm lại, hàm nguy cơ là một công cụ quan trọng trong kinh doanh để đánh giá rủi ro và dự đoán xác suất thành công của các quyết định đầu tư và các dự án. Sử dụng hàm nguy cơ có thể giúp các nhà quản lý đưa ra các quyết định thông minh dựa trên các dữ liệu và thông tin có sẵn.
Có nhiều phương pháp để ước tính hàm nguy cơ, nhưng phương pháp phổ biến nhất là mô hình Cox Proportional Hazards (mô hình tỷ lệ nguy cơ đồng nhất)
- Mô hình Cox Proportional Hazards:
Đây là một phương pháp tham số để tính toán hàm nguy cơ. Mô hình Cox Proportional Hazards được sử dụng để đánh giá tác động của các yếu tố đến xác suất sống sót của một đối tượng hoặc một nhóm đối tượng. Để tính toán hàm nguy cơ bằng mô hình Cox Proportional Hazards, trước tiên ta phải xác định các yếu tố ảnh hưởng đến xác suất sống sót và thu thập các dữ liệu liên quan đến các yếu tố đó. Sau đó, ta sử dụng các phương pháp ước tính hợp lý để tính toán các tham số của mô hình Cox Proportional Hazards. Hàm nguy cơ được tính toán từ các tham số này.
Hàm mật độ sống còn ( Hazard fuction )
Hàm mật độ sống còn (hay còn được gọi là hàm nguy cơ có điều kiện, conditional hazard function) là một khái niệm trong lý thuyết xác suất và thống kê, được sử dụng để đo lường tốc độ xảy ra của các sự kiện (ví dụ như tử vong, tai nạn, hoặc bệnh tật) trong một khoảng thời gian nhất định, dựa trên việc đã xảy ra một sự kiện khác trước đó.
Hàm mật độ sống còn là tỉ lệ giữa xác suất xảy ra sự kiện trong khoảng thời gian rất nhỏ xung quanh thời điểm t đến với độ dài khoảng thời gian đó. Cụ thể, hàm mật độ sống còn tại thời điểm t được xác định bởi tỉ lệ giữa xác suất xảy ra sự kiện trong khoảng thời gian từ t đến t + ∆t đến với độ dài khoảng thời gian đó (trong đó ∆t là một giá trị rất nhỏ).
(9)
Trong đó P(t ≤ T ≤ t + ∆t | T ≥ t) là xác suất xảy ra
sự kiện trong khoảng thời gian từ t đến t + ∆t, biết rằng sự kiện chưa
xảy ra vào thời điểm t, và T là biến cố xảy ra sự kiện. Δt là khoảng
thời gian rất nhỏ và là hàm mật độ xác suất của t.
Hàm mật độ sống còn có thể được sử dụng để đánh giá tốc độ xảy ra sự kiện tại một thời điểm cụ thể và so sánh tốc độ xảy ra sự kiện giữa các nhóm khác nhau. Nó cũng có thể được sử dụng để dự đoán xác suất xảy ra sự kiện trong tương lai.
Mối quan hệ giữa hàm sống còn và hàm nguy cơ
Chúng ta có thể sử dụng hàm rủi ro để tính toán hàm sống sót thông qua công thức:
Hàm rủi ro (h(t)) được định nghĩa là tốc độ xảy ra sự kiện tại thời điểm t. Công thức tổng quát của hàm rủi ro là:
(10)
Hai hàm sống sót và nguy cơ có mối quan hệ chặt chẽ với nhau thông qua các công thức toán học. Sử dụng hai hàm này, chúng ta có thể đánh giá tốc độ xảy ra sự kiện và xác suất sống sót tại một thời điểm cụ thể.
4 THỰC HÀNH PHÂN TÍCH SỐNG CÒN VỚI PACKAGE SURVIVAL
4.1 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.
4.2 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.
4.3 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.
4.4 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, …
4.5 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)4.6 Đườ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")
4.7 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.
4.8 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= 1654.9 Kiể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.
4.10 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.
4.11 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