Survival Data Manipulation

Giới thiệu

Bài này muốn giới thiệu các bước xử lí data để chuẩn bị cho phân tích survival.

Khi phân tích survival hàm survival chính thường là:

f1 <- survfit(Surv(time, status) ~ 1, data = lung)

Trong đó dataset là “lung”, biến số biến cố sống, chết là “status” và thời gian theo dõi của mỗi bệnh nhân từ khi khởi đầu điều trị cho đến khi có biến cố hoặc lần tiếp xúc cuối cùng của nghiên cứu là “time”.

Tuy nhiên, trong thực tế phân tích dữ liệu nghiên cứu, biến số time không có sẳn như data mẫu. Các bước xử lí data thu thập từ phiếu nghiên cứu để đi đến biến số time không hề đơn giản.

Dưới đây là bảng excel dữ liệu thời gian của các lần khám (visit) được ghi nhận theo định dạng sau. Mỗi bệnh nhân có số lần khám không đồng đều nhau. Tôi sẽ xử lí bảng excel này để cho ra các biến số ID bệnh nhân, thời gian theo dõi (time) và tình trạng sống còn (event). Mỗi row là dữ liệu của một bệnh nhân.

Tải data

library(tidyverse)
# Survival data
df <- read_excel("E:/Statistics/R/R Writings/survdat.xlsx")
df

NA

Trong bảng excel ở trên có những bệnh nhân có biến cố (event = 1) và những bệnh nhân không có biến cố (event = 0). Trong phân tích survival, biến số thời gian (time) được tính từ lúc vào nghiên cứu cho đến khi có biến cố hoặc đến khi hết theo dõi nếu bệnh nhân không có biến cố.

Như vậy, mỗi bệnh nhân, ta phải chọn ra 2 thời điểm. Thời điểm thứ nhất là lần khám 1, vừa vào nghiên cứu và thời điểm thứ 2 là ngày có biến cố (event = 1) hoặc ngày cuối cùng đối với bệnh nhân không có biến cố.

Xử lí data

Đánh dấu các lần khám của mỗi bệnh nhân như sau. Chú ý vào lần khám đầu tiên, lần khám có biến cố và lần khám cuối cùng đối với bệnh nhân không có biến cố.

df2 <- df %>% arrange(ID, visit) %>% 
  group_by(ID) %>%  
  dplyr::mutate(numb = dplyr::row_number()) %>% # Create numbering variable
  filter(numb == 1 | event == 1 | (event == 0 & numb == max(numb)))
df2

Trong bảng kết quả trên, bệnh nhân B và E không có biến cố, chúng ta sẽ chọn ra thời điểm đầu và cuối của họ. Những bệnh nhân A, C và D chúng ta chọn thời điểm vào nghiên cứu và lúc có biến cố (event =1). Vì vậy chúng ta phải loại lần khám cuối cùng của họ (A, C và D) như sau.

df2 <- df2 %>% group_by(ID) %>%
  dplyr::mutate(rank = row_number()) %>%
  filter(rank == 1 | rank == 2)
df2

NA

Đến đây đã có data như mình muốn. Chúng ta làm gọn lại, chỉ chọn những biến số cần thiết mà thôi. Sau đó để tính được khoảng thời gian giữa 2 thời điểm ta biến chúng thành 2 biến số T1 và T2.

df.s <- select(df2, ID,date,event,rank)

# as.factor(c(df.s$event, df.s$rank))
dat <- spread(df.s, rank, date) 

colnames(dat) <- c("ID", "event", "t1", "t2")
dat

Mỗi bệnh nhân vẫn có 2 rows và những bệnh nhân có biến cố thì t1 và t2 không cùng nằm trong một dòng. Ta phải sắp xếp lại lần nữa.

for (i in 1:nrow(dat)) {
  if (dat$event[i] == 1) {
    dat$t1[i] <- dat$t1[i-1]
  }
}

dat <- dat %>% drop_na()
dat

NA

Tới thời điểm này chúng ta đã có bộ data như mình mong muốn. Để tính time cho mỗi bệnh nhân, ta tính khoảng thời gian giữa 2 thời điểm t1 và t2 ở trên.

Package “lubridate” giúp chúng ta xử lí biến số thời gian.

# date process
library(lubridate)

dat <- dat %>% 
  mutate(
    t.first = ymd(t1), 
    t.last = ymd(t2)
  )

dat <- dat %>% 
  mutate(
    time = 
      as.numeric(
        difftime(t.last, 
                 t.first, 
                 units = "days")) / 365.25
  )

dat

NA

Đến đây chúng ta đã có các biến số cần thiết để phân tích survival là event (status) và time.

Cám ơn các bạn đã đọc bài viết này.