Mô hình Markov ẩn

1 Mở đầu

Câu chuyện này bắt đầu với bệnh COPD (Bệnh phổi tắc nghẽn mạn tính), vào năm 2017 tổ chức GOLD đã cập nhật tài liệu khuyến cáo về chẩn đoán, tiên lượng và điều trị bệnh COPD, trong đó đặt ra hai hệ thống phân loại riêng biệt: Độ nặng dựa vào kết quả hô hấp ký, có 4 mức độ là 1,2,3,4; và phân loại định hướng điều trị gồm 4 loại là A,B,C,D dựa vào triệu chứng lâm sàng và các biến chứng.

Trong một nghiên cứu (giả định), chúng ta có 100 bệnh nhân COPD được theo dõi trong khoảng thời gian dài qua 10 lần tái khám, mỗi lần cách nhau 4 tháng. Ở mỗi lần khám, trạng thái phân loại độ nặng theo GOLD được đánh giá lại. Như vậy, ở mỗi bệnh nhân chúng ta có dữ liệu gồm 2 chuỗi trạng thái có độ dài = 10 (Severity và Grade), đều là biến định tính rời rạc (categorical).

Câu hỏi đặt ra là: Liệu chúng ta có thể rút ra được một quy luật diễn tiến của độ nặng và trạng thái lâm sàng theo thời gian hay không ?

Nhi sẽ giới thiệu với các bạn một phương pháp thống kê cho phép trả lời câu hỏi trên, đó là mô hình Markov ẩn (Hidden Markov model, HMM), thông qua package seqHMM của Satu Helske và Jouni Helske, mới được họ công bố vào đầu năm 2019.

Trước hết chúng ta tải dữ liệu COPD_HMM.csv vào R

library(tidyverse)
library(seqHMM)
library(igraph)

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

head(df)

##   Visit Severity Grade
## 1     1       G2     B
## 2     2       G3     A
## 3     3       G2     C
## 4     4       G2     C
## 5     5       G4     B
## 6     6       G4     B

Dữ liệu có dạng bảng dọc (long format), mỗi hàng là 1 thời điểm (Visit), mỗi 10 hàng là 1 bệnh nhân, 2 cột là 2 chuỗi trạng thái Severity và Grade.

Qua thống kê mô tả với hàm group_by, ta thấy có 12 cặp tổ hợp giữa Severity và Grade:

df %>% mutate(ID = rep(c(1:100),each = 10))%>%
  group_by(Severity,Grade) %>%
  summarise(frequency = n())

## # A tibble: 12 x 3
## # Groups:   Severity [4]
##    Severity Grade frequency
##    <fct>    <fct>     <int>
##  1 G1       A            12
##  2 G1       B            51
##  3 G2       A            12
##  4 G2       B           130
##  5 G2       C           108
##  6 G3       A            11
##  7 G3       B           257
##  8 G3       C           100
##  9 G4       A             3
## 10 G4       B            26
## 11 G4       C           169
## 12 G4       D           121

2 Mô hình HMM

Sơ đồ trên trình bày về các thành phần trong một mô hình Markov ẩn, đầu ra của mô hình là 2 chuỗi trạng thái được quan sát (Y), ta giả định rằng chúng là kết quả của một quá trình Markov mang tính ngẫu nhiên (stochastic process) gồm các tham số trạng thái chưa xác định (trạng thái ẩn, X) với xác suất ban đầu (initial probability). Có sự chuyển tiếp giữa các tham số trạng thái ẩn, và xác suất đầu ra cho phép liên kết mỗi trạng thái ẩn với kết quả quan sát được.

Trong bài toán này, chúng ta có 2 chuỗi kết quả đầu ra là Severity và Grade, nên đây là một mô hình đa kênh (multichannels), ta ấn định một số lượng tham số trạng thái ẩn X (hidden states) = 12. Ta sẽ dùng mô hình HMM để xác định những xác suất chuyển tiếp.

Trước hết, ta chuyển dạng dữ liệu từ long thành wide format, là định dạng yêu cầu bởi package seqHMM, với 1 vòng lặp while :

# Setting
period = 10
start = 1

Severity = matrix(ncol = period)
Class = matrix(ncol = period )

while ((start + period -1)<= nrow(df)){
  end = (start + period) - 1
  temp_df = df[c(start:end),]
  Severity = rbind(Severity,as.vector(temp_df$Severity))%>%na.omit()
  Class = rbind(Class,as.vector(temp_df$Grade))%>%na.omit()
  start = end + 1
}

colnames(Severity) = c(paste('t',rep(1:period), sep=""))
colnames(Class) = c(paste('t',rep(1:period), sep=""))

Ta cũng tách rời 2 phần Severity và Grade trong 2 matrices riêng :

head(Severity)

##      t1   t2   t3   t4   t5   t6   t7   t8   t9   t10 
## [1,] "G2" "G3" "G2" "G2" "G4" "G4" "G4" "G4" "G3" "G3"
## [2,] "G3" "G2" "G4" "G3" "G1" "G3" "G3" "G4" "G4" "G3"
## [3,] "G2" "G3" "G4" "G4" "G4" "G2" "G3" "G4" "G3" "G3"
## [4,] "G3" "G1" "G3" "G4" "G3" "G4" "G1" "G1" "G4" "G3"
## [5,] "G2" "G4" "G3" "G3" "G3" "G1" "G4" "G4" "G3" "G2"
## [6,] "G2" "G4" "G4" "G2" "G3" "G2" "G4" "G4" "G4" "G4"

head(Class)

##      t1  t2  t3  t4  t5  t6  t7  t8  t9  t10
## [1,] "B" "A" "C" "C" "B" "B" "B" "B" "A" "A"
## [2,] "A" "B" "B" "A" "A" "A" "A" "B" "B" "A"
## [3,] "B" "A" "A" "A" "A" "B" "A" "B" "A" "B"
## [4,] "B" "B" "B" "B" "B" "B" "B" "B" "B" "B"
## [5,] "C" "B" "B" "B" "B" "B" "B" "B" "B" "B"
## [6,] "C" "B" "B" "B" "B" "C" "B" "B" "B" "B"

Bước tiếp theo, ta dùng hàm seqdef (bạn còn nhớ package TraMineR) để thiết lập 2 chuỗi dữ liệu đầu ra cho mô hình là GOLD và ABCD_Class, từ 2 matrices Severity và Class ở trên

GOLD = seqdef(Severity,start=1)
ABCD_Class = seqdef(Class, start=1)

Chuỗi dữ liệu là 1 object chứa các thông tin về cấu trúc như độ dài =10, số quan sát = 100, tên các nhãn giá trị, và phổ màu để vẽ đồ thị.

summary(GOLD)

##  [>] sequence object created with TraMineR version 2.0-11.1 
##  [>] 100 sequences in the data set, 100 unique 
##  [>] min/max sequence length: 10/10
##  [>] alphabet (state labels):  
##      1=G1 (G1)
##      2=G2 (G2)
##      3=G3 (G3)
##      4=G4 (G4)
##  [>] dimensionality of the sequence space: 30 
##  [>] colors: 1=#7FC97F 2=#BEAED4 3=#FDC086 4=#FFFF99

summary(ABCD_Class)

##  [>] sequence object created with TraMineR version 2.0-11.1 
##  [>] 100 sequences in the data set, 99 unique 
##  [>] min/max sequence length: 10/10
##  [>] alphabet (state labels):  
##      1=A (A)
##      2=B (B)
##      3=C (C)
##      4=D (D)
##  [>] dimensionality of the sequence space: 30 
##  [>] colors: 1=#7FC97F 2=#BEAED4 3=#FDC086 4=#FFFF99

Ta có thể thay đổi phổ màu của mỗi object để đạt hiệu ứng mỹ thuật tùy ý:

attr(GOLD, "cpal") <- c('#ffc2c5',
                         '#fa757c',
                         '#fc303a',
                         '#8a0007')

attr(ABCD_Class, "cpal") <- c('#98c429',
                              '#ffce08',
                              '#ff8903',
                              '#ff0385')

Hàm ssplot cho phép vẽ các biểu đồ sau, nhằm tóm tắt về phân bố và quy luật chuyển tiếp giữa các trạng thái

ssplot(list("GOLD Severity" = GOLD, 
            "ABCD_Grade" = ABCD_Class),
       plots = 'obs',
       type = 'd')

ssplot(list("GOLD Severity" = GOLD, 
            "ABCD_Grade" = ABCD_Class),
       plots = 'obs',
       sort.channel = "from.start",
       type = 'I')

Ta có thể thấy rằng ở mỗi bệnh nhân (hàng), độ năng và trạng thái lâm sàng của bệnh COPD không diễn tiến một cách tuyến tính, nhưng có thể chuyển tiếp ngược từ một trạng thái nặng hơn về trạng thái nhẹ hơn:

Để dựng mô hình HMM multichannel (đa kênh), các chuỗi Y được đóng gói vào trong 1 list và đưa vào hàm build_hmm, sau đó ta xác định số hidden states mong muốn (12), rồi dùng tiếp hàm fit_model để dựng mô hình.

Tùy chỉnh threads cho phép chạy quy trình trên nhiều CPU để tăng tốc độ tính toán, trên máy của Nhi có 8 cores nên Nhi set threads = 8.

list(GOLD, ABCD_Class) %>%
  build_hmm(observations = ., 
            n_states = 12,
            channel_names = c("GOLD","Class")) %>% 
  fit_model(., 
            em_step = TRUE, 
            global_step = TRUE,
            local_step = TRUE,
            threads = 8) -> mc_fit

Sau khi mô hình dựng xong, bạn có thể dùng hàm plot để vẽ ra mạng lưới chuyển tiếp giữa các hidden states, là kết quả cần tìm.

plot(mc_fit$model,
     edge.arrow.size = 0.5, 
     interactive = TRUE, 
     loops  = TRUE, 
     vertex.size = 15,
     edge.curved = FALSE,
     layout =  layout.star, 
     legend.prop = 0.1,
     trim = 0.01,
     with.legend = "right")

Packages HMM vẽ mạng lưới này bằng package igraph, nên bạn có thể gọi bất cứ định dạng mỹ thuật nào cho mô hình mạng từ package igraph trong tùy chỉnh layout.

plot(mc_fit$model,
     edge.arrow.size = 0.5, 
     interactive = TRUE, 
     loops  = TRUE, 
     vertex.size = 15,
     edge.curved = TRUE,
     layout =  layout.fruchterman.reingold, 
     legend.prop = 0.1,
     trim = 0.01,
     with.legend = "right")

Hàm BIC cho phép đánh giá phẩm chất mô hình:

BIC(mc_fit$model)

## [1] 4748.835

Nội dung bên trong mô hình HMM gồm những ma trận xac suất, như: matrix Initial probabilities là xác suất ban đầu của mỗi hidden state (tại thời điểm t=1), matrix Transition probabilities chứa xác suất chuyển tiếp giữa các hidden states (từ đó biểu đồ mạng được vẽ ra trong igraph, và bạn cũng có thể vẽ thủ công bằng ggraph), Sau cùng là Emission probabilities cho mỗi chuỗi đầu ra

mc_fit$model

## Initial probabilities :
##  State 1  State 2  State 3  State 4  State 5  State 6  State 7  State 8 
##   0.0430   0.1200   0.1772   0.2644   0.0174   0.0000   0.1738   0.1700 
##  State 9 State 10 State 11 State 12 
##   0.0000   0.0343   0.0000   0.0000 
## 
## 
## Transition probabilities :
##           to
## from       State 1 State 2 State 3 State 4 State 5 State 6 State 7 State 8
##   State 1   0.0637   0.000 0.00000   0.456 0.00000  0.0000  0.0000  0.4800
##   State 2   0.0000   0.000 0.33709   0.000 0.00000  0.0000  0.0000  0.0000
##   State 3   0.0000   0.114 0.49461   0.000 0.00000  0.0000  0.0000  0.0000
##   State 4   0.0343   0.000 0.00438   0.375 0.00465  0.0999  0.1123  0.3692
##   State 5   0.0889   0.000 0.30491   0.567 0.00000  0.0387  0.0000  0.0000
##   State 6   0.0000   0.000 0.00000   0.803 0.04448  0.0297  0.1226  0.0000
##   State 7   0.0440   0.000 0.00000   0.305 0.03072  0.1506  0.0930  0.1750
##   State 8   0.0333   0.000 0.00000   0.473 0.00000  0.1078  0.0981  0.2881
##   State 9   0.0171   0.214 0.51619   0.000 0.03399  0.0000  0.0283  0.0000
##   State 10  0.0000   0.000 0.00000   0.000 0.00000  0.0292  0.0000  0.0204
##   State 11  0.0000   0.000 0.54627   0.000 0.00000  0.0000  0.0000  0.0000
##   State 12  0.0000   0.575 0.00000   0.000 0.00000  0.0000  0.0000  0.0000
##           to
## from        State 9 State 10 State 11 State 12
##   State 1  0.00e+00   0.0000    0.000    0.000
##   State 2  1.43e-01   0.0000    0.232    0.288
##   State 3  2.26e-01   0.0000    0.165    0.000
##   State 4  0.00e+00   0.0000    0.000    0.000
##   State 5  0.00e+00   0.0000    0.000    0.000
##   State 6  0.00e+00   0.0000    0.000    0.000
##   State 7  0.00e+00   0.0952    0.000    0.107
##   State 8  0.00e+00   0.0000    0.000    0.000
##   State 9  1.90e-01   0.0000    0.000    0.000
##   State 10 0.00e+00   0.9504    0.000    0.000
##   State 11 3.50e-01   0.0000    0.103    0.000
##   State 12 9.53e-05   0.0000    0.000    0.425
## 
## 
## Emission probabilities :
## GOLD :
##            symbol_names
## state_names     G1     G2    G3    G4
##    State 1  1.0000 0.0000 0.000 0.000
##    State 2  0.0000 0.0000 0.000 1.000
##    State 3  0.0000 0.3342 0.666 0.000
##    State 4  0.0445 0.1215 0.834 0.000
##    State 5  1.0000 0.0000 0.000 0.000
##    State 6  0.0000 0.8715 0.128 0.000
##    State 7  0.0000 1.0000 0.000 0.000
##    State 8  0.0000 0.0000 0.000 1.000
##    State 9  0.0000 0.0000 0.000 1.000
##    State 10 0.0836 0.0678 0.390 0.459
##    State 11 0.1546 0.8454 0.000 0.000
##    State 12 0.2169 0.3229 0.460 0.000
## 
## Class :
##            symbol_names
## state_names      A     B     C      D
##    State 1  0.0758 0.924 0.000 0.0000
##    State 2  0.0000 0.000 0.000 1.0000
##    State 3  0.0000 0.000 1.000 0.0000
##    State 4  0.0000 1.000 0.000 0.0000
##    State 5  1.0000 0.000 0.000 0.0000
##    State 6  0.1661 0.834 0.000 0.0000
##    State 7  0.0000 0.208 0.792 0.0000
##    State 8  0.0000 0.000 0.975 0.0247
##    State 9  0.0000 0.000 0.000 1.0000
##    State 10 0.2231 0.777 0.000 0.0000
##    State 11 0.1147 0.885 0.000 0.0000
##    State 12 0.0000 1.000 0.000 0.0000

3 Kết luận:

Qua thí dụ mô phỏng này, các bạn đã làm quen với một số tính năng cơ bản của package seqHMM để dựng mô hình Markov ẩn đa kênh. Đây là một phương pháp có tiềm năng ứng dụng rộng lớn, chứ không chỉ dành cho bài toán phân tích chuỗi giá trị categorical. Package seqHMM còn cho phép bạn đưa vào mô hình các hiệp biến (covariates) thí dụ tuổi, giới tính, … của bệnh nhân để thực hiện một phân tích Clustering kết hợp mô hình Markov.Như vậy bạn đã có trong tay một công cụ độc đáo trong nghiên cứu trường diễn hoặc các bài toán liên quan đến dữ liệu chuỗi.