Saturday, October 3, 2015

Biểu đồ sống sót (survival graph)

Trong phân tích sống còn hay phân tích sự kiện (survival analysis), việc thể hiện dữ liệu bằng biểu đồ có khi khá nan giải. Trong R có vài hàm có thể giúp chúng ta thể hiện dữ liệu sống còn với chất lượng cao. Trong bài này tôi sẽ giới thiệu các hàm đó.


Chúng ta hãy xem một dữ liệu đơn giản của 18 bệnh nhân.

ID
Thời gian (tuần)
Bệnh (1=có, 0=không)
1
18
0
2
10
1
3
13
0
4
30
1
5
19
1
6
23
0
7
38
0
8
54
0
9
36
1
10
107
1
11
104
0
12
97
1
13
107
0
14
56
0
15
59
1
16
107
0
17
75
1
18
93
1

Chúng ta muốn thể hiện dữ liệu trên bằng biểu đồ. Trước hết, nhập dữ liệu:

Weeks = c(10, 13, 18, 19, 23, 30, 36, 38, 54,
             56, 59, 75, 93, 97, 104, 107, 107, 107)
Status = c(1, 0, 0, 1, 0, 1, 1,0, 0, 0, 1, 1, 1,
              1, 0, 1, 0, 0)

dat = data.frame(Weeks, Status)

survtime = Surv(Weeks, status==1)

Hiện nay, hàm tốt nhất để vẽ biểu đồ sống còn là "ggfortify". Có thể cài đặt package này từ github (không phải từ R) bằng devtools như sau:

library(devtools)
install_github('sinhrks/ggfortify')
library(ggfortify)

Bước kế đến là dùng package survival để tính xác suất sống còn:

library(survival)
fit = survfit(Surv(Weeks, Status) ~ 1, data = dat)

và dùng hàm autoplot trong ggfortify như sau:

autoplot(fit, xlab="Time (weeks)", ylab="Survival probability")


Thay vì vẽ xác suất không mắc bệnh, chúng ta có thể vẽ xác suất mắc bệnh, dùng option fun="event":

autoplot(fit, xlab="Time (weeks)", ylab="Survival probability", fun="event") + theme_bw()



Biểu đồ sống sót theo nhóm

Trong ví dụ trên, chúng ta chỉ có 1 nhóm bệnh nhân. Trong thực tế, có nhiều trường hợp có nhiều nhóm bệnh nhân. Chẳng hạn như chúng ta muốn tìm hiểu nguy cơ bệnh theo giới tính trong nghiên cứu dưới đây:

ID: Mã số của cá nhân
Gender: Giới tính (giá trị Male, Female)
Age: Độ tuổi  (biến liên tục)
Disease: Bệnh (yes, no)
TimeDis: Thời gian từ lúc vào nghiên cứu đến lúc mắc bệnh

Trong đó, Disease là tình trạng bệnh (yes/no) và TimeDis là thời gian từ lúc vào nghiên cứu đến lúc mắc bệnh (tính theo tháng). Chúng ta muốn thể hiện nguy cơ mắc bệnh cho nam và nữ trong một biểu đồ.

Đọc dữ liệu vào R:

dat = read.csv("~/Google Drive/_QA Book (Vietnamese)/Data/ex-data.csv", header=T)

# Code biến Disease thành biến số với giá trị 0,1.
dat$Dis = as.numeric(dat$Disease)-1

attach(dat)

Dùng survfit để tính xác suất mắc bệnh:

fit = survfit(Surv(TimeDis, Dis) ~ Gender, data = dat)

Sau đó là vẽ biểu đồ:

autoplot(fit, xlab="Time", ylab="Survival probability")



Biểu đồ xác suất tích luỹ mắc bệnh:
  
autoplot(fit, xlab="Time", ylab="Survival probability", fun="event")


Biểu đồ xác Aalen's additive regression (aareg): 


autoplot(aareg(Surv(TimeDis, Dis) ~ Age + Sex, data = dat))




1 comment:

  1. Thưa thầy Tuấn, thầy cho em hỏi mình nên diễn giải đường cong Kaplan-Meier như thế nào cho đúng ạ? Đó là xác suất (hay phần trăm) tổng số các đối tượng sống sót đến thời gian T hay là xác suất mà 1 cá nhân sống sót đến thời gian T ạ? Và cột Y-axis mình nên chú giải là " Cumulative survival probability" hay là " Survival probability" ạ? Em cảm ơn thầy ạ.

    ReplyDelete