Tạo thông số lâm sàng từ dữ liệu chuỗi thời gian

1 Đặt vấn đề

Trong một số hoàn cảnh, đại lượng sinh lý/bệnh học cần khảo sát trong nghiên cứu được ghi nhận/theo dõi liên tục trong suốt một khoảng thời gian dài bằng một dụng cụ/máy móc tự động nào đó. Thí dụ:

• Theo dõi sinh hiệu : SaO2, huyết áp, nhịp tim … của bệnh nhân trong phẫu thuật hoặc tại phòng cấp cứu.

• Đếm tần suất cơn ho, số lần sử dụng thuốc dãn phế quản dạng hít của bệnh nhân mắc bệnh lý hô hấp

• Theo dõi cơn gò tử cung ở sản phụ khi chờ chuyển dạ

• Theo dõi biến cố về nhịp tim, độ bão hòa oxy máu trong giấc ngủ tại nhà

Lúc này, dữ liệu thô mà nghiên cứu sinh có trong tay là một chuỗi thời gian. Việc khai thác được thông tin từ dữ liệu dạng này không phải là điều dễ dàng; các bạn cần sự trợ giúp của chuyên viên thống kê để có thể chuyển từ dữ liệu chuỗi thành thông số lâm sàng phù hợp với mục tiêu nghiên cứu.

Trong bài hôm nay, Nhi sẽ trình bày một ý tưởng cho phép các bạn hoán chuyển từ dữ liệu chuỗi thành những thông số lâm sàng có ý nghĩa, từ đó cho phép thực hiện những mô hình và kiểm định thống kê.

2 Thí dụ minh họa

Thí dụ minh họa trong bài như sau:

Một bệnh nhân mắc hội chứng ngưng thở khi ngủ được theo dõi độ bão hòa oxy máu ngoại biên trong khi ngủ (tần số lấy mẫu là 1 Hz, hay 1 giá trị SaO2 mỗi giây) vào đêm trước và trong khi sử dụng một máy thở Oxy áp lực dương. Câu hỏi nghiên cứu là liệu cách điều trị này có hiệu quả cải thiện sự bão hòa oxy máu của bệnh nhân trên hay không ?

Sau khi tải dữ liệu, ta sẽ có 1 dataframe 47595 hàng (tương đương 13.2 giờ), biến SaO2 là một chuỗi giá trị rời rạc, biến Cond là trạng thái Trước (Control) và trong khi điều trị (Treatment)

library(tidyverse)

df = read.csv("Saturation.csv", sep = ";")

head(df)

##      Cond Time SaO2
## 1 Control    1   96
## 2 Control    2   96
## 3 Control    3   96
## 4 Control    4   96
## 5 Control    5   96
## 6 Control    6   96

Việc trước tiên Nhi sẽ làm, đó là nhị phân hóa giá trị SaO2, bằng cách này Nhi tạo ra một Biến cố: Mất bão hòa oxy, với 2 trạng thái: Có (Y) khi SaO2 < 94 và Không (N) khi SaO2 >= 94

Như vậy, chúng ta sẽ tiếp cận bài toán bằng cả 2 loại dữ liệu: Định lượng (SaO2) và Định tính (biến cố Desaturation).

df$Desat = if_else(df$SaO2 < 94, "Y","N")

head(df)

##      Cond Time SaO2 Desat
## 1 Control    1   96     N
## 2 Control    2   96     N
## 3 Control    3   96     N
## 4 Control    4   96     N
## 5 Control    5   96     N
## 6 Control    6   96     N

Bây giờ, mời các bạn cùng Nhi quan sát diễn tiến của SaO2 trong thời giàn hơn 6h, trong 2 điều kiện khác nhau

df %>% ggplot(aes(x = Time/3600, y = SaO2, col = Cond))+
  geom_step(direction = "hv", alpha = 0.4)+
  geom_smooth(se = F)+
  geom_hline(yintercept = 94, linetype = 2)+
  scale_color_manual(values = c("red","blue"))+
  xlab('Duration (h)')+
  theme_minimal()+
  facet_wrap(~Cond,ncol = 1)

Qua hình ảnh này, ta có cảm giác rằng khi điều trị bằng máy thở oxy, SaO2 trung bình có vẻ cao hơn; tuy nhiên làm cách nào để mô tả và kiểm tra điều này bằng con số và thống kê ?

Ta không thể dùng giá trị SaO2 ở thời gian thực trong cả 1 đêm làm biến số để so sánh, như hình vẽ bên dưới cho thấy: rất khó phân biệt 2 histogram phân bố của chuỗi SaO2

df %>% ggplot(aes(x = SaO2, fill = Cond))+
  geom_histogram(alpha = 0.5, binwidth = 1, col = "black")+
  geom_vline(xintercept = 94, linetype = 2)+
  scale_fill_manual(values = c("red","blue"))+
  facet_wrap(~Cond,ncol = 1)+
  theme_minimal()

Ta cũng không thể sử dụng giá trị trung vị hay trung bình làm đại diện cho chuỗi SaO2, như biểu đồ sau đây cho thấy không có sự khác biệt về vị trí trung tâm (cũng như hầu hết bách phân vị) giữa 2 chuỗi SaO2

library(lvplot)

df %>% ggplot(aes(x = Cond, y = SaO2, fill = Cond))+
  geom_lv(aes(alpha = ..LV..), col = "black", show.legend = F)+
  geom_hline(yintercept = 94, linetype = 2)+
  coord_flip()+
  scale_fill_manual(values = c("red","blue"))+
  theme_minimal()

3 Phép đếm tần suất và tái chọn mẫu

Các bạn còn nhớ ở trên ta đã tạo ra 1 chuỗi nhị phân là biến cố Desaturation, bây giờ các bạn thử tưởng tượng: Nếu ta đặt ra một thông số lâm sàng có tên là: Tần suất phát sinh biến cố mất bão hòa Oxy trong 1 giờ, hay số cơn giảm SpO2 trong 1h, đây có thể là một thông số có ý nghĩa đại diện và thực dụng.

Nhưng làm sao xác định được thông số này ? Cách làm đơn giản nhất, đó là ta đếm tổng số cơn mất bão hòa oxy trong cả đêm, rồi chia cho tổng thời gian, thí dụ ta đếm được 100 cơn, và thời gian ngủ là 6.5 giờ, ta sẽ có chỉ số = 15.38 cơn/h

Cách làm này đơn giản, nhưng có nhiều vấn đề như:

• Không xét đến phân bố của biến cố theo thời gian : các biến cố có thể rải rác trong 6h, hoặc chỉ tập trung trong một số giai đoạn nhất định

• Không xét đến thời gian kéo dài của mỗi biến cố: 5-10 biến cố kéo dài vài phút hẳn là nguy hiểm hơn 20-30 biến cố ngắn.

Do đó, Nhi sẽ viết một function với công dụng : Ước tính tần suất phát sinh biến cố Desaturation trong mỗi giờ và độ dài chuẩn hóa của desaturation trong 1h đó dựa vào chọn mẫu ngẫu nhiên.

Như vậy, Nhi xem chuỗi nhị phân Desaturation với hàng ngàn thời điểm là một không gian chọn mẫu, rồi chọn ngẫu nhiên một vị trí bất kì trên chuỗi này, đi đến phía trước 3599 giây; đếm tất cả những biến cố Desaturation = Y trong khoảng thời gian 1h này, tính tổng độ dài của chúng. Lặp lại hàng ngàn lần quy trình này, ta sẽ có thể mô tả được phân bố của 2 đại lượng mà ta cần khảo sát là tần suất phát sinh biến cố mỗi giờ, và độ dài của biến cố.

boot_count = function(org = NULL,replace = T,size = 10, samp_window = 3600){
  
  sample_idx = sample(c(1:(nrow(org)- samp_window - 1)), 
                      replace = replace, 
                      size = size)
  
  out_df = data.frame(Event = NA,
                      Duration = NA,
                      Norm_dt = NA,
                      Freq = NA,
                      Loc = NA)
  
  i = 1
  
  while (i <= length(sample_idx)){
    
    samp_df = org[c(sample_idx[i] : (sample_idx[i]+samp_window - 1)),]
    
    change.idx <- with(samp_df,which(head(Desat, -1) != tail(Desat, -1)))
    
    if(length(change.idx) > 0){
      
      events = samp_df[c(change.idx,max(change.idx)+1),]%>%.$Desat
      
      locs = samp_df[c(change.idx,max(change.idx)+1),]%>%.$Time
      
      start = c(1, change.idx+1)
      
      end = c(change.idx, nrow(samp_df))
      
      duration = 1 + end - start
      
      temp_df = data.frame('Event'= events,
                          'Duration' = duration)
      
      temp_df%>%group_by(Event)%>%
        summarise_at("Duration", sum)%>%
        mutate(Norm_dt = Duration/sum(Duration),
               Freq = as.vector(table(temp_df$Event)),
               Loc = mean(locs))-> sum_df
    }
    
    else{sum_df = data.frame(Event = samp_df[1,]%>%.$Desat,
                             Duration = samp_window,
                             Norm_dt = 1,
                             Freq = 1,
                             Loc = mean(samp_df$Time))}
    
    out_df = bind_rows(out_df, sum_df)
    
    i = i+1
    }
    
  return(out_df %>%na.omit())}

Nhi áp dụng quy trình nói trên 10 ngàn lần cho 2 chuỗi dữ liệu Control và Treatment, kết quả xuất ra như sau:

df %>% filter(Cond == "Control") %>% 
  boot_count(org = .,
             replace = F,
             size = 10000,
             samp_window = 3600)%>%
  mutate(Cond = "Baseline") -> out_base

head(out_base)

##   Event Duration    Norm_dt Freq       Loc     Cond
## 1     N     3425 0.95138889   35  9025.029 Baseline
## 2     Y      175 0.04861111   34  9025.029 Baseline
## 3     N     3426 0.95166667   15 20779.828 Baseline
## 4     Y      174 0.04833333   14 20779.828 Baseline
## 5     N     2494 0.69277778  114 15749.982 Baseline
## 6     Y     1106 0.30722222  114 15749.982 Baseline

df %>% filter(Cond != "Control") %>% 
  boot_count(org = .,
             replace = F,
             size = 10000,
             samp_window = 3600)%>%
  mutate(Cond = "Treatment") -> out_treat

head(out_treat)

##   Event Duration     Norm_dt Freq      Loc      Cond
## 1     N     3580 0.994444444    4 17028.14 Treatment
## 2     Y       20 0.005555556    3 17028.14 Treatment
## 3     N     3600 1.000000000    1 19645.50 Treatment
## 4     N     3470 0.963888889   10 14986.35 Treatment
## 5     Y      130 0.036111111   10 14986.35 Treatment
## 6     N     3148 0.874444444   46  6859.12 Treatment

Lúc này, ta đã có trong tay 2 bộ dữ liệu tái chọn mẫu cho 2 thông số : tần suất phát sinh Desaturation trong 1h, và độ dài của desaturation trong 1h. Với vị trí của cửa sổ chọn mẫu trên trục thời gian, ta có thể mô tả diễn tiến 2 thông số này:

Theo đó, ban đầu bệnh nhân dường như có số cơn Desaturation tương đương hoặc thậm chí cao hơn, nhưng sau khi ngủ 2h thì điều trị bắt đầu cho thấy hiệu quả: số cơn desaturation giảm về 0 trong suốt 2h, có tăng nhẹ lúc 4-5h nhưng sau đó =0 cho đến sáng.

comb_df = bind_rows(out_base, out_treat)

comb_df %>% filter(Event == 'Y')%>% 
  ggplot(aes(x = Loc/3600, y = Freq, col = Cond)) +
  geom_step()+
  xlab("Timeline (h)")+
  ylab('Hourly occuring rate (event/h)')+
  theme_bw()

Trong khi đó, tỉ lệ thời gian Desaturation (Đỏ) so với Bình thường (Xanh) cũng giảm đáng kể khi điều trị

out_base %>% 
  ggplot(aes(x = Loc/3600, 
             y = Norm_dt, 
             fill = Event)) +
  geom_bar(stat = "identity", 
           position = "fill")+
  scale_fill_manual(values = c("green","red"))+
  xlab("Duration (h)")+
  theme_bw()

out_treat %>% 
  ggplot(aes(x = Loc/3600, 
             y = Norm_dt, 
             fill = Event)) +
  geom_bar(stat = "identity", 
           position = "fill")+
  scale_fill_manual(values = c("green","red"))+
  xlab("Duration (h)")+
  theme_bw()

Bản thân thời gian chuẩn hóa của Desaturation cũng theo diễn tiến tương tự như chỉ số tần suất phát sinh biến cố vậy:

comb_df %>% filter(Event == 'Y')%>% 
  ggplot(aes(x = Loc/3600, y = Norm_dt, col = Cond)) +
  geom_step()+ 
  xlab("Timeline (h)")+
  ylab('Normalized event duration')+
  theme_bw()

Ưu điểm của quy trình tái chọn mẫu ngẫu nhiên, đó là nó cho phép mô tả đặc tính phân bố của 2 thông số lâm sàng này, và qua so sánh hình ảnh 2 phân bố này, ta thấy được sự tương phản giữa Trước và Trong khi điều trị:

Theo đó, điều trị bằng máy CPAP làm giảm đáng kể tần suất phát sinh biến cố Desaturation / 1h:

comb_df %>% filter(Event == 'Y')%>% 
  ggplot(aes(x = Freq)) +
  geom_histogram(alpha = 0.5, aes(fill = Cond), col = "black", binwidth = 5)+
  facet_wrap(~Cond, ncol = 1)+
  xlab('Hourly occuring rate (event/h)')+
  scale_fill_manual(values = c("red","blue"))+
  theme_bw()

Nó cũng làm giảm tổng thời gian Desaturation / mỗi giờ

comb_df %>% filter(Event == 'Y')%>% 
  ggplot(aes(x = Norm_dt)) +
  geom_density(alpha = 0.5, aes(fill = Cond), col = "black", binwidth = 5)+
  xlab("Normalized duration")+
  scale_fill_manual(values = c("red","blue"))+
  theme_bw()

Ta có thể thực hiện suy diễn thống kê cho dữ liệu tái chọn mẫu bằng mô hình, thí dụ mô hình sau đây đánh giá hiệu ứng điều trị lên tần suất phát sinh desaturation mỗi giờ - được khảo sát như biến số đếm; mô hình dựa vào phân bố Negative Binomial (Nhị thức âm) có xét thêm Spline function cho thời gian.

Dựa vào suy diễn qua Incidence rate ratio, có thể nói điều trị làm giảm khoảng 50% nguy cơ phát sinh Desaturation mỗi giờ:

library(gamlss)

mod = gamlss(data = comb_df[comb_df$Event == "Y",], 
             formula = Freq ~ Cond + bs(Loc),
             family = "NBII")

## GAMLSS-RS iteration 1: Global Deviance = 178120 
## GAMLSS-RS iteration 2: Global Deviance = 178101.5 
## GAMLSS-RS iteration 3: Global Deviance = 178100.6 
## GAMLSS-RS iteration 4: Global Deviance = 178100.5 
## GAMLSS-RS iteration 5: Global Deviance = 178100.5 
## GAMLSS-RS iteration 6: Global Deviance = 178100.5

summary(mod)

## ******************************************************************
## Family:  c("NBII", "Negative Binomial type II") 
## 
## Call:  
## gamlss(formula = Freq ~ Cond + bs(Loc), family = "NBII", data = comb_df[comb_df$Event ==  
##     "Y", ]) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.172536   0.021760  191.75   <2e-16 ***
## CondTreatment -0.714591   0.009304  -76.81   <2e-16 ***
## bs(Loc)1      -0.537974   0.057974   -9.28   <2e-16 ***
## bs(Loc)2       0.358312   0.027597   12.98   <2e-16 ***
## bs(Loc)3      -0.498868   0.028602  -17.44   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.90988    0.01084   268.6   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  19673 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  19667 
##                       at cycle:  6 
##  
## Global Deviance:     178100.5 
##             AIC:     178112.5 
##             SBC:     178159.8 
## ******************************************************************

Tương tự, ta có thể dùng mô hình với phân bố Beta cho thời gian desaturation /h chuẩn hóa (dao động từ 0 đến 1). Mô hình cho thấy điều trị làm giảm 0.46 lần thời gian Desaturation / h

mod = gamlss(data = comb_df[comb_df$Event == "Y",], 
             formula = Norm_dt ~ Cond + bs(Loc),
             family = "BE")

## GAMLSS-RS iteration 1: Global Deviance = -48227.25 
## GAMLSS-RS iteration 2: Global Deviance = -53240.21 
## GAMLSS-RS iteration 3: Global Deviance = -54749.65 
## GAMLSS-RS iteration 4: Global Deviance = -54993.09 
## GAMLSS-RS iteration 5: Global Deviance = -55018.68 
## GAMLSS-RS iteration 6: Global Deviance = -55020.91 
## GAMLSS-RS iteration 7: Global Deviance = -55021.09 
## GAMLSS-RS iteration 8: Global Deviance = -55021.11 
## GAMLSS-RS iteration 9: Global Deviance = -55021.11 
## GAMLSS-RS iteration 10: Global Deviance = -55021.11

summary(mod)

## ******************************************************************
## Family:  c("BE", "Beta") 
## 
## Call:  
## gamlss(formula = Norm_dt ~ Cond + bs(Loc), family = "BE", data = comb_df[comb_df$Event ==  
##     "Y", ]) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   -1.31864    0.02660 -49.581  < 2e-16 ***
## CondTreatment -0.72049    0.01131 -63.731  < 2e-16 ***
## bs(Loc)1      -1.73454    0.07146 -24.272  < 2e-16 ***
## bs(Loc)2       0.18401    0.03432   5.361 8.37e-08 ***
## bs(Loc)3      -0.53798    0.03406 -15.793  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  logit
## Sigma Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.975601   0.005721  -170.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  19673 
## Degrees of Freedom for the fit:  6
##       Residual Deg. of Freedom:  19667 
##                       at cycle:  10 
##  
## Global Deviance:     -55021.11 
##             AIC:     -55009.11 
##             SBC:     -54961.79 
## ******************************************************************

4 Trung bình và tái chọn mẫu

Ở trên, chúng ta tạo ra thông số lâm sàng dựa vào chuỗi giá trị định tính, bây giờ ta sẽ xét đến chuỗi định lượng, là chính bản thân SaO2.

Áp dụng cơ chế tương tự, Nhi sẽ tạo ra một chỉ số lâm sàng có tên là: SaO2 trung bình / 10 phút; các bạn còn nhớ mình đã ghi nhận sinh hiệu cho bệnh nhân bằng cách đo đếm thủ công trong 1 phút, ở đây ta cũng làm tương tự, khi đặt một cửa sổ chọn mẫu ngẫu nhiên dài 10 phút (600 quan sát) tại những vị trí bất kì trên chuỗi thời gian. Để làm việc này, Nhi viết thêm 1 hàm nữa như sau:

boot_mean = function(org = NULL,replace = T,size = 10, samp_window = 600){
  
  sample_idx = sample(c(1:(nrow(org)- samp_window - 1)), 
                      replace = replace, 
                      size = size)
  
  out_df = data.frame(Mean_SaO2 = rep(NA,size),Loc = rep(NA,size))
  
  i = 1
  
  while (i <= length(sample_idx)){
    samp_df = org[c(sample_idx[i] : (sample_idx[i]+samp_window - 1)),]
    
    out_df$Mean_SaO2[i] = mean(samp_df$SaO2)
    out_df$Loc[i] = mean(samp_df$Time)
    
    i = i + 1
  }
  
  return(na.omit(out_df))
}

Nhi lần lượt áp dụng hàm boot_mean này cho 2 chuỗi Control và Treatment, lưu ý tùy chỉnh replace = T cho thấy đây là 1 quy trình Bootstrap 10 ngàn lượt

df %>% filter(Cond != "Control") %>% 
  boot_mean(org = .,
             replace = T,
             size = 10000,
             samp_window = 600)%>%
  mutate(Cond = "Treatment") -> sa_tm

head(sa_tm)

##   Mean_SaO2     Loc      Cond
## 1  94.37667 15611.5 Treatment
## 2  94.57667  2218.5 Treatment
## 3  94.62500  8858.5 Treatment
## 4  96.07833 18331.5 Treatment
## 5  96.14000 20063.5 Treatment
## 6  95.29833 13773.5 Treatment

df %>% filter(Cond == "Control") %>% 
  boot_mean(org = .,
             replace = T,
             size = 10000,
             samp_window = 600)%>%
  mutate(Cond = "Control") -> sa_base

head(sa_base)

##   Mean_SaO2     Loc    Cond
## 1  93.09500 16690.5 Control
## 2  94.93500  3639.5 Control
## 3  95.32333 17744.5 Control
## 4  95.86833   598.5 Control
## 5  95.01500  2732.5 Control
## 6  93.92833 15188.5 Control

Tương tự, ta có thể khảo sát hình ảnh phân bố của chỉ số lâm sàng vừa tạo ra là SaO2 trung bình trong 10 phút (lưu ý là bạn có thể chọn kích thước cửa sổ này tùy ý, 1 phút, 5 phút, 10 phút, thậm chí 30 phút)

comb_sa = bind_rows(sa_base, sa_tm)

comb_sa %>% ggplot(aes(x = Mean_SaO2, fill = Cond))+
  geom_histogram(alpha = 0.5, col = "black")+
  xlab('SaO2 / 10 min')+
  theme_bw()

comb_sa %>%
  ggplot(aes(x = Cond, y = Mean_SaO2)) +
  geom_lv(alpha = 0.5, aes(fill = Cond), col = "black")+
  coord_flip()+
  theme_bw()

Ta cũng không mất thông tin về yếu tố thời gian, vì vị trí chọn mẫu được lưu lại trong kết quả, ta có thể khảo sát diễn tiến của SaO2/10ph

comb_sa %>% 
  ggplot(aes(x = Loc/3600, y = Mean_SaO2, col = Cond)) +
  geom_step()+ 
  xlab("Timeline (h)")+
  ylab('SaO2 / 10 min')+
  theme_bw()

5 Bàn luận

Thông điệp quan trọng nhất mà Nhi muốn chia sẻ với các bạn qua bài này, đó là kỹ thuật lập trình thống kê có tiềm năng vô cùng lớn. Nó không chỉ cho phép lập kế hoạch phân tích trên những bảng dữ liệu thông số lâm sàng có sẵn, mà còn cho phép tạo ra cả những đại lượng hoàn toàn mới từ dữ liệu thô. Cơ chế được sử dụng trong bài rất đơn giản, đó là tái chọn mẫu ngẫu nhiên hoặc bootstrap, với một ít kỹ thuật về vòng lặp trong R nhưng kết quả tạo ra là khá ấn tượng. Ghi chú thêm là thí dụ minh họa này chỉ mới là một chuỗi nhị phân, tuy nhiên kỹ thuật tái chọn mẫu có thể mở rộng cho chuỗi đa sự kiện (nhiều loại biến cố),

Chúc các bạn vững bước trên hành trình khoa học với R