ĐỀ TÀI: PHÂN TÍCH SỐNG CÒN - PACKAGE SURVIVAL

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 sót_ Phân tích sống sót_ Phân tích sự kiện

   Phân tích sống sót 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ì tô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 Phân tích sống sót 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 mang 2 giá trị đối lập (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 sự kiện (sống/chết, …) 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ố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 sót.

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 gian. Phân tích sống sót đượ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 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 và cách thực hành trên R dựa trên package survival.

   Đặc biệt, trong thời đại hiện nay vừa trải qua dịch bệnh Covid-19 gây ảnh hưởng rất lớn đến sức khỏe, kinh tế và xã hội của toàn thế giới. Do đó việc áp dụng phân tích sống sót cho dịch bệnh Covid 19 là việc có ý nghĩa rất quan trọng. Cụ thể, phân tích sống sót Covid-19 có thể đưa ra các dự đoán về số ca nhiễm và tử vong trong tương lai dựa trên các mô hình dự đoán, đánh giá khả năng ảnh hưởng của các biện pháp phòng chống dịch bệnh và đưa ra các giải pháp để giảm thiểu tác động của dịch bệnh đến sức khỏe và kinh tế.

   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 sót 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 sót 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 sót đượ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 sót 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 sót đượ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 sót 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 ta sẽ 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.

   Sau đó chúng ta sẽ 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 sót đượ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.

   Nội dung chính của bài tiểu luận này là tìm hiểu về package survival và thực hành trên R với package survival, vì vậy phần lý thuyết về phân tích sống sót sẽ không được giải thích sâu.

2 TỔNG QUAN VỀ PACKAGE SURVIVAL

2.1 Thông tin package survival

   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 sót đượ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 sót (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 sót. 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 sót đ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 sót 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ượng Surv trong 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 sót là gì?

   Phân tích sống sót (survival analysis) là một phương pháp thống kê để đánh giá và dự đoán khả năng sống sót hoặc thời gian đến sự kiện trong một tập dữ liệu.

   Phân tích sống sót được sử dụng rộng rãi trong nhiều lĩnh vực, bao gồm y học, y tế công cộng, kinh tế học, khoa học xã hội, khoa học môi trường, khoa học dữ liệu và nghiên cứu thị trường. Nó cung cấp một cách tiếp cận thống kê để đánh giá khả năng sống sót hoặc thời gian đến sự kiện trong một tập dữ liệu, và đưa ra các dự đoán về khả năng sống sót hoặc thời gian sống của một nhóm bệnh nhân hoặc cá nhân trong tương lai.

   Vậy, phân tích sống sót (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 sự kiện (time-to-failure/ time-to-event).

   Các thuật ngữ cơ bản của phân tích sống sót bao gồm:

• Dữ liệu thời gian và biến cố (time-to-event)

• 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 xảy ra sự kiện 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 sót (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.

   Ví dụ về các dữ liệu thời gian -sự kiện (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 sót, 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.

   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:

   Bản tóm tắt về các loại kiểm chứng:

• 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 sót

   Hàm sinh tồn

   Đường cong tỷ lệ sống sót 𝑆(𝑡) biểu thị tỷ lệ sống sót là một hàm của thời gian (t):

𝑆(𝑡)=𝑃𝑟(𝑇>𝑡)

   Giá trị của hàm S(t) là tỷ lệ sống sót của một nhóm vào thời điểm (t) so với thời điểm bắt đầu quan sá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 , ∞ ) .

• 𝑆(𝑡) luôn giảm, nghĩa là 𝑆(𝑡1)≥𝑆(𝑡2) với 𝑡1≤𝑡2

• Tạ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.

   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

   Ướ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 𝑡𝑖

   Hàm nguy cơ

   Hàm nguy cơ (hazard function) là một khái niệm quan trọng trong phân tích sống sót (survival analysis), được sử dụng để mô tả tốc độ xảy ra sự kiện trong một nhóm dân số theo thời gian.

   Công thức của hàm nguy cơ thường được biểu diễn dưới dạng:

h(t) = lim Δt → 0 [P(t ≤ T < t + Δt | T ≥ t) / Δt]

   Trong đó:

• t là thời gian quan sát, T là thời gian xảy ra sự kiện, Δt là một khoảng thời gian rất nhỏ, tiến tới 0.

• 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 đã không xảy ra trước thời điểm t

   Giá trị của hàm h(t) là tốc độ xảy ra sự kiện (event rate) tại thời điểm t. Ví dụ, nếu h(10) = 0.05 nghĩa là tại thời điểm 10 đơn vị thời gian sau khi bắt đầu quan sát, tốc độ xảy ra sự kiện là 0.05 sự kiện trên một đơn vị thời gian.

   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 để 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.

4 THỰC HÀNH VỚI PACKAGE SURVIVAL - PHÂN TÍCH SỐNG SÓT

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 gọi package survival bằng lệnh sau:

# Gọi package survival từ thư viện
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:

# Hiển thị các thành phần của dữ liệu
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

1. inst: Mã số của bệnh nhân.

2. time: Thời gian sống sót hoặc thời gian theo dõi (ngày).

3. status: Tình trạng sống sau thời gian theo dõi (1 = Sống sót, 2 = Qua đời).

4. age: Tuổi của bệnh nhân.

5. sex: Giới tính của bệnh nhân (1 = Nam, 2 = Nữ).

6. 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).

7. 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)

8. pat.karno: Điểm Karnofsky của bệnh nhân dựa trên người chăm sóc (nếu có).

9. meal.cal: Lượng calo tiêu thụ hàng ngày của bệnh nhân.

10. 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 sót, 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:

# Gọi package dplyr từ thư viện
library(dplyr)
# Chuyển đổi giá trị của biến status 1 => 0 và 2 => 1 
lung <- lung %>% mutate(status = recode(status, '1' = 0, '2' = 1))
# Hiển thị 6 quan sát đầu tiên của 3 biến quan tâm
head(lung[c("time","status","sex")])

##   time status sex
## 1  306      1   1
## 2  455      1   1
## 3 1010      0   1
## 4  210      1   1
## 5  883      1   1
## 6 1022      0   1

   Ngoài ra, chúng ta sẽ chia tuổi của bệnh nhân thành 3 nhóm cho tiện quá trình phân tích:

# Chia tuổi thành 3 nhóm tuổi 
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õ"))
# Hiển thị 6 quan sát đầu tiên của 4 biến quan tâm
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

   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 tuổi, 50 đến 65 tuổi, 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 sót, để 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 sót (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 sót 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:

# Lấy biến time và status từ dataset lung và tạo đối tượng sinh tồn mới
Surv(lung$time, lung$status)[1:10]

##  [1]  306   455  1010+  210   883  1022+  310   361   218   166

   Ta 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 sống sót

   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 sót 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 sót 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 sống sót Surv(lung$time, lung$status) để tạo các đường cong sống sót 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():

# Tính toán đường cong sống sót của tất cả bệnh nhân
survfit_lung <- survfit(Surv(time, status) ~ 1, data = lung)
# Hiển thị các thành phần của dữ liệu 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 sống sót 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: Uớ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():

# Tính toán đường cong sống sót theo giới tính
survfit_lung_bysex <- survfit(Surv(time, status) ~ sex, data = lung)
# Hiển thị các thành phần của dữ liệu 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 (ph.ecog) và nhóm tuổi (age_group):

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 sống sót Kaplan-Meier

   Để vẽ đường cong sống sót 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àm survfit() dùng để tính toán đường cong sinh tồn (xem phần 4.4)

   Đường cong sống sót 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 sống sót 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)

=> Nhận xét đồ thị: đường cong sống sót của dữ liệu ung thư phổi - lung là một đường cong với xu hướng giảm liên tục.

   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 sống sót Kaplan Meier theo yếu tố giới tính

   Nâng cao hơn, ta có thể vẽ đường cong sống sót Kaplan-Meier cho từng nhóm được phân loại theo yếu tố. Trong ví dụ này, đường cong sống sót Kaplan-Meier được vẽ mỗi nhóm giới tính khác nhau trong dữ liệu ung thư phổi:

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")

=> Nhận xét đồ thị: có sự khác biệt về xác suất sống sót giữa 2 nhóm nhớm tính (p=0.0013 < 0.05) và giới tính nam (sex=1) có xác suất sóng sót qua thời gian thấp hơn giới tính nữ (sex=1)

   Đườ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-lung:

ggsurvplot(survfit_lung_byph, data = lung, risk.table = TRUE, pval = TRUE, legend.title = "Mức độ suy giảm chức năng phổi", xlab = "Thời gian (ngày)", ylab = "Xác suất sóng sót")

=> Nhận xét đồ thị: có sự khác biệt đáng kể giữa các mức độ suy giảm chức năng phổi (p~0.0001 < 0.05) và mức độ suy giảm chức năng phổi càng càng cao thì xuất suất sóng sốt qua thời gian càng thấp.

   Đườ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:

ggsurvplot(survfit_lung_byage, data = lung, risk.table = TRUE, pval = TRUE, legend.title = "Nhóm tuổi", xlab = "Thời gian (ngày)", ylab = "Xác suất sóng sót")

=> Nhận xét đồ thị: không có sự khác nhau đáng kể giữa các nhóm tuổi (p= 0.18 > 0.05) khi xác suất sống sót giữa các độ tuổi trong 500 ngày đầu gần như không có sự khác biệt.

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 và mỗi nhóm tuổi khác nhau trong dữ liệu ung thư phổi như sau:

# So sánh 2 nhóm giới tính của dữ liệu lung
logrank_test1 <- survdiff(with(lung, Surv(time, status)) ~ sex, data = lung)
# So sánh 4 nhóm mức độ suy giảm chức năng phổi 
logrank_test2 <- survdiff(with(lung, Surv(time, status)) ~ ph.ecog, data = lung)
# So sánh 3 nhóm độ tuổi 
logrank_test3 <- survdiff(with(lung, Surv(time, status)) ~ age_group, 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.001

print(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-05

print(logrank_test3)

## Call:
## survdiff(formula = with(lung, Surv(time, status)) ~ age_group, 
##     data = lung)
## 
##                            N Observed Expected (O-E)^2/E (O-E)^2/V
## age_group=50 đến 65 tuổi 116       81     87.1     0.428     0.916
## age_group=dưới 50 tuổi    20       11     15.4     1.238     1.373
## age_group=trên 65 tuổi    92       73     62.5     1.753     2.855
## 
##  Chisq= 3.5  on 2 degrees of freedom, p= 0.2

   Kế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.

• Nhóm Độ tuổi ( dưới 50, 50 đến 65, trên 65): giá trị p= 0.2 > 0.05 nên có thể kết luận rằng không có sự khác biệt đáng kể giữa 3 nhóm tuổi được xé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 sót để đá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:

# Tạo mô hình cox theo giới tính
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= 165

4.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 package 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() (xem phần 4.8), ta kiểm định như sau:

# Kiểm định mô hình cox
cox.zph(coxfit, transform = "km")

##        chisq df     p
## sex     2.86  1 0.091
## GLOBAL  2.86  1 0.091

   Kế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:

# Gọi package gtsummary từ thư viện
library("gtsummary")
# Tạo mô hình cox theo giới tính và trả về tỷ lệ rủi ro HR
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: Giá trị p=0.001<0.05 nên kết quả trên là có ý nghĩa thống kê. 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.

   Áp dụng với yếu tố mức độ suy giảm chức năng phổi ph.ecog:

# Tạo mô hình cox theo mức độ suy giảm chức năng phổi và trả về tỷ lệ rủi ro HR
coxph(Surv(time, status) ~ ph.ecog, data = lung) %>% tbl_regression(exp = TRUE) 

Characteristic	HR1	95% CI1	p-value
ph.ecog	1.61	1.29, 2.01	<0.001
1 HR = Hazard Ratio, CI = Confidence Interval

• Kết quả từ hàm tbl_regression() cho thấy: Giá trị p=0.001<0.05 nên kết quả trên là có ý nghĩa thống kê. HR = 1.61 ngụ ý cho thấy nguy cơ tử vong tăng lên theo mức độ suy giảm chức năng phổi và với mỗi mức suy giảm chức năng phổi tăng lên thì xuất sắc bệnh nhân sắp chết có thể tăng theo 1,61 lần.

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:

# Tạo mô hình cox để dự báo
mod <- coxph(Surv(time,status) ~+ sex, data = cancer)
# Tạo 1 dòng quan sát mới để dự báo
pred_dat <- data.frame(time = c(1050), status = c(0), sex = c(0))
# Tiến hành dự báo bằng hàm predict
preds <- predict(mod, newdata = pred_dat,
                 type = "survival", se.fit = TRUE)
# Lọc xác suất sống sót từ biến preds
pred_dat$prob <- preds$fit
pred_dat

##   time status sex        prob
## 1 1050      0   0 0.002573077

Kế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

4.12 Kết quả phân tích dữ liệu ung thư phổi - lung

   Thông qua các bước thực hành phân tích sống sót trên R với dữ liệu ung thư phổi - lung, chúng ta có bản tóm tắt các kết quả:

1. 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,…
2. Biến time, status là tổ hợp biến time-to-event (dữ liệu thời gian - sự kiện) được dùng để phân tích sống sót trong R. Bên cạnh đó, biến age, sex, ph.ecog là 3 biến yếu tố dùng để phân nhóm và đánh giá ảnh hưởng của chúng lên tỉ lệ sống sót.
3. Đồ thị sống sót tổng thể của bệnh nhân ung thư phổi:

• Đồ thị đường cong sóng sót là đồ thị liên tục giảm, không có xu hướng tăng ở bất kỳ thời điểm nào vì bệnh nhân chết đi không thể sống lại.

4. Yếu tố giới tính có ảnh hưởng đến tỉ lệ sóng sót của bệnh nhân, trong đó bệnh nhanh mang giới tính nam có tỉ lệ sóng sót qua thời gian thấp hơn so với bệnh nhân mang giới tính nữ.

5. Yếu tố mức độ suy giảm chức năng phổi có ảnh hưởng đến tỉ lệ sống sót của bệnh nhân, trong đó bệnh nhân có mức độ suy giảm phổi càng lớn thì tỉ lệ sống sót qua thời gian càng thấp.

6. Yếu tố độ tuổi gần như không có ảnh hưởng đến tỉ lệ sống sót của bệnh nhân.

5 THỰC HÀNH NÂNG CAO VỚI PACKAGE SURVIVAL

5.1 Nguy cơ cạnh tranh - competing risks

   Competing risks (nguy cơ cạnh tranh) là một khái niệm trong phân tích sống sót, đề cập đến việc có nhiều hơn một sự kiện có thể xảy ra đối với các cá nhân trong quá trình theo dõi. Trong phân tích sống sót, sự kiện quan tâm thường là sự kiện sớm nhất xảy ra trong quá trình theo dõi, nhưng trong trường hợp có nhiều sự kiện có thể xảy ra, ví dụ như:

• Tái phát.

• Chết do bệnh.

• Chết vì yếu tố khác.

   Như vậy, phân tích nguy cơ cạnh tranh được sử dụng để tính toán tỷ lệ nguy cơ của các sự kiện cạnh tranh và đo lường ảnh hưởng của chúng đến sự kiện quan tâm.

5.1.1 Dataset melanoma - ung thư da

   Để phân tích nguy cơ cạnh tranh, chúng ta sẽ sử dụng dữ liệu về ung thu da - melanoma trong package MASS bằng cách thực hiện như sau:

# Sử dụng dataset melanoma từ package MASS
data(Melanoma, package = "MASS")
# Hiển thị 6 quan sát đầu tiên
head(Melanoma)

##   time status sex age year thickness ulcer
## 1   10      3   1  76 1972      6.76     1
## 2   30      3   1  56 1968      0.65     0
## 3   35      2   1  41 1977      1.34     0
## 4   99      3   0  71 1968      2.90     0
## 5  185      1   1  52 1965     12.08     1
## 6  204      1   1  28 1971      4.84     1

   Bộ dữ liệu gồm 205 quan sát và 7 biến sau:

• time: Thời gian (ngày) từ lúc chẩn đoán ung thư đến sự kiện quan tâm (tái phát ung thư hoặc tử vong) hoặc kết thúc theo dõi.

• status: tình trạng: 1= chết vì ung thư da, 2= còn sống, 3= chết vì nguyên nhân khác.

• sex: Giới tính của bệnh nhân.

• age: Tuổi của bệnh nhân.

• thíckness: độ dày khối u ban đầu (tính bằng mm).

• ulcer: Có hay không vết loét trên da ban đầu.

Mã hóa lại dữ liệu của biến status

   Vì bộ dữ liệu về ung thư da- melanoma có biến status với giá trị chưa được mã hóa chuẩn, chúng ta sẽ chuyển đổi giá trị trong biến thành các giá trị mới để tránh nhầm lẫn:

# Chuyển đổi dữ liệu của biến status 2=>0, 3=>2
Melanoma <- 
  Melanoma %>% 
  mutate(status = as.factor(recode(status, `2` = 0, `1` = 1, `3` = 2)))
# Hiển thị 6 quan sát đầu tiên
head(Melanoma)

##   time status sex age year thickness ulcer
## 1   10      2   1  76 1972      6.76     1
## 2   30      2   1  56 1968      0.65     0
## 3   35      0   1  41 1977      1.34     0
## 4   99      2   0  71 1968      2.90     0
## 5  185      1   1  52 1965     12.08     1
## 6  204      1   1  28 1971      4.84     1

Bây giờ chúng ta có:

• status: 0= còn sống, 1= chết vì ung thư da, 2= chết vì các nguyên nhân khác.

5.1.2 Tỷ lệ tích lũy - Cumulative incidence

   Tỷ lệ tích lũy - Cumulative incidence là một khái niệm trong phân tích nguy cơ cạnh tranh, được sử dụng để tính toán tỷ lệ xảy ra của nhiều sự kiện quan tâm (ví dụ: tái phát ung thư và chết) trong một nhóm cụ thể tại một thời điểm nhất định trong quá trình theo dõi.

Công thức tính Tỷ lệ tích lũy được cho bởi: Cumulative incidence = tổng số sự kiện/số nguy cơ.

Để uớc tính tỷ lệ tích lũy trong rủi ro cạnh tranh, chúng ta thực hiện bằng cách sử dụng hàm cuminc() từ package tidycmprsk:

# Gọi package tydycmprsk từ thư viện
library("tidycmprsk")
# Ước tính tỷ lệ tích lũy
cuminc(Surv(time, status) ~ 1, data = Melanoma)

## 
## time    n.risk   estimate   std.error   95% CI          
## 1,000   171      0.127      0.023       0.086, 0.177    
## 2,000   103      0.230      0.030       0.174, 0.291    
## 3,000   54       0.310      0.037       0.239, 0.383    
## 4,000   13       0.339      0.041       0.260, 0.419    
## 5,000   1        0.339      0.041       0.260, 0.419    
## 
## time    n.risk   estimate   std.error   95% CI          
## 1,000   171      0.034      0.013       0.015, 0.066    
## 2,000   103      0.050      0.016       0.026, 0.087    
## 3,000   54       0.058      0.017       0.030, 0.099    
## 4,000   13       0.106      0.032       0.053, 0.179    
## 5,000   1        0.106      0.032       0.053, 0.179

   Chúng ta có thể sử dụng hàm ggcuminc() từ package ggsurvfit để vẽ biểu đồ tỷ lệ tích lũy. Theo mặc định, nó chỉ vẽ sơ đồ loại sự kiện đầu tiên. Vì vậy, biểu đồ sau đây cho thấy tỷ lệ tử vong tích lũy do ung thư da:

# Gọi package ggsurvfit từ thư viện
library("ggsurvfit")
# Vẽ đồ thị với bảng at rick và khoảng tin cậy 
cuminc(Surv(time, status) ~ 1, data = Melanoma) %>% 
  ggcuminc() + 
  labs(
    x = "Days"
  ) + 
  add_confidence_interval() +
  add_risktable()

   Ngoài ra, chúng ta có thể thể hiện cả hai loại sự kiện là chết vì ung thư da và chết vì nguyên nhân khác, ta thực hiện như sau:

# tỷ lệ tích lũy nguy cơ của chết vì ung thư da và chết vì nguyên nhân khác
cuminc(Surv(time, status) ~ 1, data = Melanoma) %>% 
  ggcuminc(outcome = c("1", "2")) +
  ylim(c(0, 1)) + labs(x = "Days")

   Bây giờ, giả sử chúng ta muốn kiểm tra cái chết do ung thư da hoặc các nguyên nhân khác trong dữ liệu theo biến ulcer: Có hay không vết loét trên da ban đầu (0=không, 1=có). Chúng ta có thể ước tính tỷ lệ tích lũy tại các thời điểm khác nhau theo nhóm và hiển thị tỷ lệ đó trong bảng cách sử dụng hàm tbl_cuminc() từ gói package tidycmprsk để kiểm tra sự khác biệt giữa các nhóm trong toàn bộ thời gian theo dõi bằng cách sử dụng hàm add_p():

# tỷ lệ tích lũy nguy cơ của chết vì ung thư da và chết vì nguyên nhân khác nhưng có yếu tố vết loét
cuminc(Surv(time, status) ~ ulcer, data = Melanoma) %>% 
  tbl_cuminc(
    times = 1826.25, 
    label_header = "**{time/365.25}-year cuminc**") %>% 
  add_p()

Characteristic	5-year cuminc	p-value1
ulcer		<0.001
0	9.1% (4.6%, 15%)
1	39% (29%, 49%)
1 Gray’s Test

   Sau đó, chúng ta có thể thấy đồ thị về cái chết do ung thư da, theo tình trạng vết loét, chúng ta sử dụng ggcuminc() từ package ggsurvfit:

cuminc(Surv(time, status) ~ ulcer, data = Melanoma) %>% 
  ggcuminc() + 
  labs(
    x = "Days"
  ) + 
  add_confidence_interval() +
  add_risktable()

## Plotting outcome "1".

5.1.3 Hồi quy Cox trong nguy cơ cạnh tranh

   Giả sử chúng ta quan tâm đến việc xem xét ảnh hưởng của tuổi tác và giới tính đối với cái chết do ung thư da, với cái chết do các nguyên nhân khác là một sự kiện cạnh tranh.

Hàm crr() từ package tidycmprsk sẽ ước tính các nguy cơ phân phối phụ:

# tạo mô hình cox với 2 yếu tố là giới tính và tuổi
crr(Surv(time, status) ~ sex + age, data = Melanoma)

## 
## Variable   Coef    SE      HR     95% CI       p-value    
## sex        0.588   0.272   1.80   1.06, 3.07   0.030      
## age        0.013   0.009   1.01   0.99, 1.03   0.18

Và chúng ta có thể tạo các bảng kết quả được định dạng bằng cách sử dụng hàm tbl_regression() từ package gtsummary, với tùy chọn exp = TRUE để có được ước tính tỷ lệ rủi ro:

# + trả về tỷ lệ nguy cơ tích lũy
crr(Surv(time, status) ~ sex + age, data = Melanoma) %>% 
  tbl_regression(exp = TRUE)

Characteristic	HR1	95% CI1	p-value
sex	1.80	1.06, 3.07	0.030
age	1.01	0.99, 1.03	0.2
1 HR = Hazard Ratio, CI = Confidence Interval

   Kết quả cho thấy: Giới tính nam (1=nam, 0=nữ) có liên quan đáng kể đến việc tăng nguy cơ tử vong do khối u ác tính (vì p-value=0.03<0.05), trong khi tuổi tác không liên quan đáng kể đến tử vong do khối u ác tính (vì p-value=0.2>0.05).

5.2 Trực quan hóa đường cong tỉ lệ tử vong - Fit complex survival curves

   “Fit complex survival curves” có thể được dịch sang tiếng Việt là “Đường cong tỉ lệ tử vong”. Đây là một thuật ngữ trong phân tích sống còn (survival analysis) để mô hình hóa dữ liệu thời gian sự kiện phức tạp, trong đó sự kiện có thể xảy ra ở nhiều lần khác nhau hoặc có nhiều nguyên nhân khác nhau dẫn đến sự kiện xảy ra.

   Chúng ta vẫn sẽ sử dụng dữ liệu ung thư phổi - lung cho việc trực quan hóa đường tỉ lệ tử vong. Chúng ta sẽ tạo một bảng các đường cong sống sót theo giới tính được phân chia theo 2 yếu tố. Bảng các đồ thị trực quan được thực hiện như sau:

fitx <- survfit( Surv(time, status) ~ sex + age_group + ph.ecog,
                data = lung )

ggsurv <- ggsurvplot(fitx, fun = "event", conf.int = TRUE,
                     ggtheme = theme_bw())
   
ggsurv$plot +theme_bw() + 
  theme (legend.position = "right")+
  facet_grid(age_group ~ ph.ecog)

   Từ bảng các đồ thị, ta thấy đồ thị 1 thuộc nhóm 50 đến 60 tuổi và có mức độ suy giảm chức năng phổi = 0 và có 2 đường cong tỉ lệ tử vong là giới tính nam và nữ.

6 ÁP DỤNG PHÂN TÍCH SỐNG SÓT VỚI PACKAGE SURVIVAL

6.1 Định phí bảo hiểm

Chúng ta sẽ sử dụng dữ liệu về xác suất sống sót của dữ liệu ung thư phổi (xem phần 4.1) để áp dụng vào định giá bảo hiểm điều trị ung thư phổi.

Bằng cách xem xét các đồ thị đường cong sống sót trong R và phân tích các chỉ số liên quan, tôi đã xây dựng một mô hình định giá bảo hiểm hỗ trợ điều trị ung thư phổi một cách chính xác.

Mô hình này cho phép tính toán mức độ rủi ro và xác định các khoản bồi thường phù hợp dựa trên xác suất sống sót của người được bảo hiểm. Việc áp dụng dữ liệu xác suất sống sót vào định giá bảo hiểm giúp tăng tính công bằng và chính xác trong việc định giá các chính sách bảo hiểm điều trị ung thư.

Như vậy, công thức định phí bảo hiểm hỗ trợ điều trị ung thư phổi có thể sử dụng một phương trình đơn giản:

Phí bảo hiểm = Phí điều trị mỗi ngày x Thời gian sống sót trung bình

Chúng ta sẽ đặt giả thiết như sau:

• Không có tình trạng lạm phát xảy ra.

• Phí điều trị như nhau ( ví dụ 500k/ngày)

• Đối tượng mua bảo hiểm: bệnh nhân bị ung thư - hỗ trợ chi phí điều trị cho tới khi chết.

Như vậy, để định phí bảo hiểm hỗ trợ ung thư phổi, ta cần tính thời gian sống sót trung bình:

fit <- coxph(Surv(time, status) ~ 1, data = lung)
print(median(survfit(fit)$time))

## [1] 274

Kết quả: Phí bảo hiểm = 500.000*274= 137.000.000

Như vậy, chúng ta cần bán bảo hiểm hỗ trợ điều trị ung thơ với mức giá trên 137 triệu

Ngoài ra, bằng cách sử dụng package survival và phân tích sống sót, chúng ta có thể đưa ra phí bảo hiểm khác nhau cho từng nhóm khách hàng.

6.2 Ngành sản xuất sản phẩm

Trong sản xuất sản phẩm, phân tích sống sót có thể được áp dụng để xác định chất lượng của sản phẩm dựa trên việc kiểm tra một số lượng sản phẩm mẫu.

Để thực hiện phân tích sống sót cho ngành sản xuất sản phẩm, tôi sẽ dùng dữ liệu lightbulb về 2 nhóm bóng đèn: huỳnh quang và sợi đốt

library(readxl)
lightbulb <- read_excel("lightbulb.xlsx")

Lightbulb là dữ liệu được tôi lấy từ dữ liệu mẫu của phần mềm MATLAB (phần mềm cung cấp môi trường tính toán số và lập trình, do công ty MathWorks thiết kế) và đã được tôi lưu vào excel

lightbulb <- lightbulb %>% mutate(status = recode(status, '1' = 0, '0' = 1))
head(lightbulb)

## # A tibble: 6 × 4
##      id  time status  type
##   <dbl> <dbl>  <dbl> <dbl>
## 1     1  8148      1     0
## 2     2  8503      1     0
## 3     3  5184      1     0
## 4     4 10994      1     0
## 5     5 10463      0     0
## 6     6 12860      1     0

Trong đó:

• Time: Dữ liệu tuổi thọ (tính bằng giờ) của hai loại bóng đèn
• Status: Thông tin kiểm duyệt với 0= chưa xuất hiện lỗi và 1= đã xuất hiện lỗi
• Type: loại bóng đèn với 0= đèn huỳnh quang, 1= đèn sợi đốt

Ta có đồ thị đường cong sống sót của 2 loại bóng đèn như sau:

bongden <- survfit(Surv(time, status) ~ type, data = lightbulb)
plot(bongden, main = "Đường cong không gặp lỗi của bóng đèn ", xlab = "Thời gian (giờ)", ylab = "Xác suất không gặp lỗi")

Chúng ta có thể thấy rằng xác suất gặp lỗi của bóng đèn sợi đốt cao hơn nhiều so với bóng đèn huỳnh quang.

7 PHẦN KẾT LUẬN

7.1 Ưu điểm của package survival

• Phù hợp với nhiều loại dữ liệu: Package “survival” có thể được sử dụng cho nhiều loại dữ liệu phân tích sống còn, bao gồm dữ liệu censored, dữ liệu sự kiện rời rạc và dữ liệu không censored.

• Độ tin cậy cao: Các mô hình survival phân tích dữ liệu thời gian đến sự kiện thường được sử dụng trong các nghiên cứu y tế và kinh tế học, và package “survival” cung cấp các công cụ và hàm để phân tích các mô hình này với độ tin cậy cao.

• Cung cấp các công cụ để đánh giá mối quan hệ giữa các biến và thời gian cho đến sự kiện: Package “survival” cung cấp các công cụ để đánh giá mối quan hệ giữa các biến và thời gian cho đến sự kiện, bao gồm phân tích Cox proportional hazards và tính toán survival function.

• Tích hợp tốt với các package khác trong R: Package “survival” tích hợp tốt với các package khác trong R, cho phép phân tích dữ liệu phân tích sống còn và đánh giá mối quan hệ giữa các biến và thời gian cho đến sự kiện trong nhiều ngữ cảnh khác nhau.

7.2 Nhược điểm của package survival

• Khó sử dụng cho người mới bắt đầu: Package “survival” có nhiều hàm và tham số phức tạp, do đó có thể khó sử dụng cho người mới bắt đầu sử dụng R hoặc phân tích dữ liệu thời gian đến sự kiện.

• Đòi hỏi kiến thức về phân tích sống sót: Package “survival” đòi hỏi kiến thức chuyên sâu về phân tích dữ liệu phân tích sống sót và các mô hình survival, do đó không phù hợp cho những người không có nền tảng kiến thức đầy đủ.

• Có thể gặp vấn đề với dữ liệu thiếu: Package “survival” có thể gặp vấn đề khi phân tích dữ liệu thời gian đến sự kiện có dữ liệu thiếu hoặc dữ liệu censored nhiều, do đó cần phải xử lý kỹ trước khi phân tích.

• Không phù hợp cho các nghiên cứu quan sát ngắn hạn: Package “survival” thường được sử dụng để phân tích dữ liệu thời gian đến sự kiện trong các nghiên cứu quan sát dài hạn, do đó không phù hợp cho các nghiên cứu quan sát ngắn hạn.

7.3 Ý nghĩa của phân tích sống còn

• Trong y học: Phân tích sống sót được sử dụng trong nghiên cứu y học để đánh giá tác động của các yếu tố khác nhau đến sự kiện quan trọng như tử vong, bệnh lý, tái phát bệnh, và thời gian sống sót. Nó có thể giúp các nhà nghiên cứu phân tích và đưa ra dự báo về nguy cơ mắc bệnh, tốc độ tiến triển bệnh và dự báo thời gian sống sót.

• Trong kinh tế học: Phân tích sống sót cũng được sử dụng trong kinh tế học để đánh giá tác động của các chính sách và các yếu tố khác đến sự kiện quan trọng như thời gian thất nghiệp, thời gian nghỉ hưu, và thời gian kinh doanh của một doanh nghiệp.

• Trong nghiên cứu xã hội: Phân tích sống sót cũng có giá trị trong nghiên cứu xã hội, giúp các nhà nghiên cứu đánh giá tác động của các yếu tố khác đến sự kiện quan trọng như thời gian chuyển đổi việc làm, thời gian kết hôn hoặc ly hôn, và thời gian sống sót của một nhóm người.

  Tuy nhiên, bài tiểu luận này chỉ thực hành phân tích survival analysis với dữ liệu thuộc về y học.

8 TÀI LIỆU THAM KHẢO

1. Phân tích sống còn - Cẩm nang dịch tễ học với R (Applied Epi)

2. A package for survival analysis in R - Terry Therneau

3. Survival Analysis in R - Emily C.Zabor