Phân tích chuỗi trạng thái

1 Giới thiệu

Trong một cuộc thăm dò thường niên và kéo dài từ 1999-2009 của tổ chức Swiss Household Panel, 2612 người dân Thụy Sỹ đã tự đánh giá sức khỏe của bản thân theo thang điểm gồm 5 mức độ: G1=Rất khỏe mạnh (Very well), G2=Khỏe (Well), M=Trung bình (Average), B1=Không tốt và B2=Xấu (Bad).

Nội dung của dữ liệu như sau:

library(tidyverse)
library(PST)

data(SRH)

SRH%>%head()%>%knitr::kable()

idpers	sex	birthy	wp09lp1s	p99c01	p00c01	p01c01	p02c01	p03c01	p04c01	p05c01	p06c01	p07c01	p08c01	p09c01
4101	man	1965	0.6870344	very well	very well	very well	very well	very well	very well	very well	well	very well	very well	well
4102	woman	1968	0.9610977	not very well	well	so, so (average)	very well	very well	very well	so, so (average)	so, so (average)	not very well	not very well	well
5101	man	1961	0.8828794	well	well	well	well	well	NA	well	well	very well	well	well
27101	woman	1968	0.9888913	so, so (average)	so, so (average)	well	well	so, so (average)	well	well	well	well	so, so (average)	so, so (average)
27102	man	1968	0.8139466	NA	very well	well	very well	very well	well	very well	well	so, so (average)	well	very well
34101	woman	1965	0.4232759	very well	very well	very well	well	very well	NA	well	well	well	well	well

Thiết kế nghiên cứu này hết sức thú vị : Ở thời điểm khởi đầu, chưa có một giả thuyết xác định nào: Mục tiêu có thể là khai phá dữ liệu: tìm ra một quy luật diễn tiến (pattern) của sức khỏe của quần thể theo thời gian. Với quy luật này, ta hy vọng có thể tiên lượng được sức khỏe của một cá thể tại một thời điểm bất kì nếu biết thông tin về người này trong quá khứ.

Ý tưởng khảo sát chuỗi trạng thái/sự kiện có thể áp dụng trong nhiều lĩnh vực, thí dụ:

1. Thương mại: Một cửa hàng điện thoại muốn khảo sát khuynh hướng lựa chọn thay đổi smartphone của khách hàng theo thời gian, trong số các thương hiệu, từ đó có thể tiên lượng xác suất một người đang sử dụng máy Samsung sẽ chuyển sang dùng IPhone hoặc tiếp tục trung thành với thương hiệu cũ.

2. Y học lâm sàng: Theo dõi diễn tiến của một bệnh lý theo thời gian, thí dụ lộ trình diễn tiến của bệnh suy tim (theo NYHA), bệnh phổi tắc nghẽn mạn tính, hoặc trạng thái dinh dưỡng (bình thường/thừa cân/béo phì…) ; hoặc sự thay đổi phương pháp điều trị (chuyển từ thuốc dãn phế quản sang corticoid, từ kháng sinh thường sang kháng sinh tân tiến hơn, giữa các loại thuốc giảm đau…)

3. Xã hội học: Theo dõi diễn tiến của tình trạng Gia đình (Độc thân, kết hôn, có con…), học vấn/nghề nghiệp (xác suất vào đại học, và học tiếp lên cao học hay đi làm, nguy cơ thất nghiệp, sự chuyển đổi nghề, công việc), tôn giáo, …

Vấn đề này hoàn toàn mới lạ đối với nhiều bạn, và yêu cầu một phương pháp thống kê đặc biệt có tên là “Phân tích chuỗi trạng thái”, mà mục tiêu là khai thác thông tin từ dữ liệu trường diễn (longitudinal) với outcome là một biến rời rạc nhiều bậc giá trị (trạng thái hay sự kiện).

Bài thực hành hôm nay, Nhi sẽ giới thiệu với các bạn công cụ chuyên dụng cho loại nghiên cứu này, đó là package TRaMineR.

2 Cấu trúc dữ liệu chuỗi

Đầu tiên, chúng ta cần bàn về cấu trúc của dữ liệu cho thiết kế nghiên cứu theo dõi chuỗi trạng thái. Vấn đề là dữ liệu loại này có thể được tạo ra dưới nhiều dạng và chúng khác nhau về cấu trúc. Ngoài các thuộc tính như Giới tính, Năm sinh, … , dữ liệu chuỗi có 3 thành phần chính: Biến định danh cá thể , Thời gian/Thời điểm và trạng thái.

1. STS: Dạng thứ nhất tên là State-sequence, có cấu trúc Wide, mỗi hàng là một cá thể i, cột đầu tiên để định danh, các thời điểm và trạng thái được trình bày bằng nhiều cột nối tiếp nhau, mỗi cột là một thời điểm t và mỗi ô nhận 1 giá trị trạng thái (ij)

data_frame(Id=c(1:5),T1=c("A","A","B","C","D"),
           T2=c("A","B","C","D","D"),
           T3=c("A","B","B","A","C")
           )%>%knitr::kable()

Id	T1	T2	T3
1	A	A	A
2	A	B	B
3	B	C	B
4	C	D	A
5	D	D	C

2. SPS: State-permanence. Dạng dữ liệu này cũng có cấu trúc wide, mỗi hàng cũng trình bày 1 cá thể, nhưng mỗi cột biểu thị cho một bậc trạng thái và nhận giá trị gồm 2 thông tin (Tên trạng thái/thời gian kéo dài cho trạng thái đó).

Ưu điểm của định dạng này là cho phép khảo sát thêm thời gian như một biến liên tục, nhưng khó thao tác và quản lý.

data_frame(Id=c(1:5),State1=c("A/3",NA,"B/4","C/2","D/1"),
           State2=c("C/2","B/3","A/4","D/1",NA),
           State3=c("D/3","C/1",NA,"D/3","A/4")
           )%>%knitr::kable()

Id	State1	State2	State3
1	A/3	C/2	D/3
2	NA	B/3	C/1
3	B/4	A/4	NA
4	C/2	D/1	D/3
5	D/1	NA	A/4

3. TSE: Time-stamped event. Trái với 2 dạng trên, TSE có cấu trúc Long, gồm 3 cột: 1 cột định danh, 1 cột để đo thời gian như con số thực sự, 1 cột chỉ sự kiện. Cấu trúc dữ liệu này rất dễ thao tác, nhập liệu và cho phép khảo sát thời gian như một biến định lượng (thí dụ tuổi, ngày/ tháng), cũng như theo dõi chuỗi sự kiện.

data_frame(Id=c(1,1,1,2,2,2,3,3,3),
           Age=c(24,25,26,22,23,24,31,32,33),
           Event=c("A","B","C","A","C","D","B","C","D")
           )%>%knitr::kable()

Id	Age	Event
1	24	A
1	25	B
1	26	C
2	22	A
2	23	C
2	24	D
3	31	B
3	32	C
3	33	D

4. SPELL: Đây là một dạng mở rộng từ TSE, nó cũng có cấu trúc Long, gồm 4 cột, một cột định danh, một cột chỉ sự kiện, 1 cột chỉ thời điểm phát sinh sự kiện và 1 cột chỉ thời gian chấm dứt sự kiện. Đây là định dạng dữ liệu linh hoạt nhất, nó cho phép khảo sát cả những sự kiện trùng hợp, xảy ra đồng thời, cũng như khảo sát được cả khoảng thời gian và thời điểm.Một biến phụ Duration có thể được tạo ra để đo khoảng thời gian sự kiện diễn ra. Cấu trúc dữ liệu này có thể được tạo ra tự động từ dữ liệu chuỗi thời gian bằng một hàm; khi biết giá trị start point và end point.

data_frame(Id=c(1,1,1,2,2,2,3,3,3),
           From=c(18,23,28,16,26,33,21,35,48),
           To=c(22,27,40,25,32,48,34,47,52),
           Event=c("A","B","C","A","C","D","B","C","D"),
           Duration=c(To-From)
           )%>%knitr::kable()

Id	From	To	Event	Duration
1	18	22	A	4
1	23	27	B	4
1	28	40	C	12
2	16	25	A	9
2	26	32	C	6
2	33	48	D	15
3	21	34	B	13
3	35	47	C	12
3	48	52	D	4

Trở lại với thí dụ trong bài, ta có thể thấy nó có dạng STS.

Package TraMineR cho phép chuyển đổi thoải mái giữa các định dạng dữ liệu, thí dụ ta có thể chuyển từ STS sang SPELL:

library(TraMineR)

srh_seq=seqdef(SRH[,5:15])

seqformat(srh_seq, 
            from = "STS", to = "SPELL",
              process=F,
                with.missing=TRUE)%>%head()%>%knitr::kable()

id	begin	end	states
1	1	7	very well
1	8	8	well
1	9	10	very well
1	11	11	well
2	1	1	not very well
2	2	2	well

Hay chuyển từ STS sang SPS:

seqformat(srh_seq, 
            from = "STS", to = "SPS",
              process=F,
                with.missing=TRUE)%>%head()%>%knitr::kable()

	[1]	[2]	[3]	[4]	[5]	[6]	[7]	[8]	[9]	[10]	[11]
[1]	(very well,7)	(well,1)	(very well,2)	(well,1)	NA	NA	NA	NA	NA	NA	NA
[2]	(not very well,1)	(well,1)	(so, so (average),1)	(very well,3)	(so, so (average),2)	(not very well,2)	(well,1)	NA	NA	NA	NA
[3]	(well,5)	(*,1)	(well,2)	(very well,1)	(well,2)	NA	NA	NA	NA	NA	NA
[4]	(so, so (average),2)	(well,2)	(so, so (average),1)	(well,4)	(so, so (average),2)	NA	NA	NA	NA	NA	NA
[5]	(*,1)	(very well,1)	(well,1)	(very well,2)	(well,1)	(very well,1)	(well,1)	(so, so (average),1)	(well,1)	(very well,1)	NA
[6]	(very well,3)	(well,1)	(very well,1)	(*,1)	(well,5)	NA	NA	NA	NA	NA	NA

Cũng như từ STS sang TSE

ttrans <- seqetm(srh_seq, method='transition')

seqformat(srh_seq, 
            from = "STS", to = "TSE",tevent=ttrans,
              process=F,
                with.missing=TRUE)%>%head()%>%knitr::kable()

id	time	event
1	0	very well
1	7	very well>well
1	8	well>very well
1	10	very well>well
2	0	not very well
2	1	not very well>well

Trên thực tế, do 2 định dạng SPELL và TSE là dễ quản lý nhất, nên chúng ta thường nhập liệu dưới dạng này, nhưng định dạng STS lại là mặc định cho nhiều hàm trong package, nên ta phải chuyển dạng ngược từ SPELL sang STS hay từ TSE sang STS. Để làm việc này, ta chỉ cần khai báo 3 biến: id,begin/end(cho SPELL) hay time (cho TSE) và event/status.

3 Thống kê mô tả

Đầu tiên, từ dữ liệu SRH, ta phải tạo ra chuỗi trạng thái/sự kiện, đây sẽ là nguyên liệu cho tất cả mọi quá trình phân tích tiếp theo nên bước này rất quan trọng:

Bên trong object chuỗi này, ta có thể xác định thêm: Nhãn tên gọi cho mỗi sự kiện/trạng thái, thí dụ B1,B2,M,G1,G2 (theo thứ tự ABC), và phổ màu dùng để vẽ biểu đồ sau này:

srh_seq<-seqdef(SRH,5:15,
                states=c("B1","B2","M","G1","G2"),
                labels=c("NotW","Bad","Average","VeryW","Well"),
                cpal=pals::brewer.rdylgn(5))

summary(srh_seq)

##  [>] sequence object created with TraMineR version 2.0-8 
##  [>] 2612 sequences in the data set, 1901 unique 
##  [>] min/max sequence length: 10/11
##  [>] alphabet (state labels):  
##      1=B1 (NotW)
##      2=B2 (Bad)
##      3=M (Average)
##      4=G1 (VeryW)
##      5=G2 (Well)
##  [>] dimensionality of the sequence space: 44 
##  [>] colors: 1=#D7191C 2=#FDAE61 3=#FFFFBF 4=#A6D96A 5=#1A9641 
##  [>] symbol for missing state: * 
##  [>] symbol for void element: %

Object này có dạng 1 matrix hay dataframe, tức là bạn có thể thao tác với nó thoải mái như bất kì matrix hay dataframe nào khác, thí dụ đổi tên biến bằng chuỗi thời gian từ năm 1999 đến 2009

colnames(srh_seq)=c(1999:2009)

srh_seq%>%head()%>%knitr::kable()

1999	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009
G1	G1	G1	G1	G1	G1	G1	G2	G1	G1	G2
B1	G2	M	G1	G1	G1	M	M	B1	B1	G2
G2	G2	G2	G2	G2	*	G2	G2	G1	G2	G2
M	M	G2	G2	M	G2	G2	G2	G2	M	M
*	G1	G2	G1	G1	G2	G1	G2	M	G2	G1
G1	G1	G1	G2	G1	*	G2	G2	G2	G2	G2

Ta có thể tóm tắt xem có bao nhiêu bậc trạng thái trong dữ liệu này:

seqstatl(srh_seq)

## [1] "%"  "*"  "B1" "B2" "G1" "G2" "M"

alphabet(srh_seq)

## [1] "B1" "B2" "M"  "G1" "G2"

Ta có thể khảo sát một số cá thể ngẫu nhiên trong dữ liệu dưới hình thức SPS:Định dạng SPS như đã giới thiệu bên trên, cho ta biết trình tự các trạng thái trong chuỗi cho mỗi cá thể, và thời gian cho mỗi trạng thái như vậy

print(srh_seq[1:5, ], format = "SPS")

##     Sequence                                                           
## [1] (G1,7)-(G2,1)-(G1,2)-(G2,1)                                        
## [2] (B1,1)-(G2,1)-(M,1)-(G1,3)-(M,2)-(B1,2)-(G2,1)                     
## [3] (G2,5)-(*,1)-(G2,2)-(G1,1)-(G2,2)                                  
## [4] (M,2)-(G2,2)-(M,1)-(G2,4)-(M,2)                                    
## [5] (*,1)-(G1,1)-(G2,1)-(G1,2)-(G2,1)-(G1,1)-(G2,1)-(M,1)-(G2,1)-(G1,1)

Tiếp theo, chúng ta sẽ vẽ một số biểu đồ, để khảo sát trực quan thông tin về chuỗi trạng thái trong dữ liệu

Biểu đồ thứ nhất mô tả tỉ lệ phân bố các trạng thái trong toàn bộ mẫu:nó có bản chất là 1 biểu đồ Stackedbar, xếp cạnh nhau

par(mfrow = c(1,2))

seqdplot(srh_seq, 
         main = "States distribution",
         with.legend = F,border=NA)

seqlegend(srh_seq)

Ta có thể vẽ thủ công biểu đồ này bằng ggplot2:

pat<-seqstatd(SRH.seq)%>%.$Frequencies

colnames(pat)<-c(1:11)

patdf<-pat%>%as.tibble()%>%mutate(Status=rownames(pat))

longdf<-patdf%>%gather(`1`:`11`,key="Timepoint",value="Dist")

longdf$Timepoint=as.integer(longdf$Timepoint)

longdf%>%ggplot(aes(x=Timepoint,
                    y=Dist,fill=Status))+
  geom_bar(stat="identity",position="fill",width=1,col="black")+
  scale_fill_manual(values = pals::brewer.spectral(5))+theme_bw()

longdf%>%ggplot(aes(x=Timepoint,
                    y=Dist,fill=Status))+
  geom_bar(stat="identity",position="fill",width=1,col="black")+
  theme_bw()+coord_polar(theta="x")+
  scale_fill_manual(values = pals::brewer.spectral(5))

longdf%>%ggplot(aes(x=Timepoint,
                    y=Dist,fill=Status))+
  geom_bar(stat="identity",position="identity",width=1,col="black")+
  theme_bw()+coord_polar(theta="y")+
  scale_fill_manual(values = pals::brewer.spectral(5))

Biểu đồ này chỉ mới cung cấp thông tin về tỉ lệ của các trạng thái một cách tổng quát (cho quần thể), theo đó có thể thấy người ta có khuynh hướng tự đánh giá sức khỏe của mình một cách lạc quan vừa phải (mức độ G2 hoặc M).

Một dạng biểu đồ khác có bản chất là 1 heatmap, cho phép khảo sát trình tự diễn biến cho từng cá thể, nhưng nó chỉ có thể trình bày mỗi lần tối đa 10 cá thể:

Kết quả trên 10 người đầu tiên trong danh sách cho thấy đa số có trạng thái sức khỏe hằng năm dao động giữa G1 và G2, chỉ có 1 ngoại lệ là người số 2 có những thay đổi khá đột ngột từ B1 sang G2 và M.

par(mfrow = c(1,2))

seqiplot(srh_seq, 
         main = "Index plot (first 10 persons)",
         with.legend = F)

seqlegend(srh_seq)

Dạng biểu đồ thứ ba, hữu dụng hơn cả, vì nó cho phép khảo sát khuynh hướng diễn tiến chung trong toàn bộ 2612 chuỗi trạng thái. Theo đó, khuynh hướng phổ biến nhất là một sự ổn định ở mức G1-G2

par(mfrow = c(1,2))

seqfplot(srh_seq, 
           main = "Sequence frequency",
            with.legend = F,
         pbarw = T)

seqlegend(srh_seq)

Ta còn có thể so sánh đặc tính phân bố các sự kiện/trạng thái giữa 2 phân nhóm, thí dụ giới tính:

seqmtplot(srh_seq, group = SRH$sex)

Và tương tự, so sánh khuynh hướng diễn tiến chung giữa 2 giới tính

seqiplot(srh_seq, group = SRH$sex, idxs = 0, border = NA,space = 0)

4 Khai thác sâu hơn tính chất chuỗi trạng thái

Trong bước tiếp theo, chúng ta bắt đầu khảo sát sâu hơn về diễn tiến của các chuỗi trạng thái trong dữ liệu: Mục tiêu là tính toán vài trị số thống kê và làm các phân tích xác suất.

Trị số thống kê đầu tiên có tên là Entropy (theo Shannon) cho mỗi chuỗi. Entropy đo lường tính chất “bất ổn định” của việc tiên lượng các trạng thái trong chuỗi. Nếu một chuỗi chỉ chứa toàn 1 trạng thái, entropy=0; entropy của một chuỗi có độ dài 12 thời điểm và có 4 trạng thái sẽ đạt cực đại = 1.386 nếu mỗi trạng thái xảy ra đúng 3 lần.

srh_seq%>%seqient()%>%head()

##     Entropy
## 1 0.2945993
## 2 0.8530953
## 3 0.2019854
## 4 0.4281055
## 5 0.5861353
## 6 0.4181657

srh_seq%>%seqient()%>%apply(.,2,mean)

##   Entropy 
## 0.3714447

Kết quả trên trình bày giá trị Entropy của 6 người đầu tiên, và giá trị Entropy trung bình cho toàn bộ mẫu là 0.3714

Ghi chú là entropy không xét đến thứ tự thay đổi của các sự kiện/trạng thái, nhưng chỉ đếm số lần trạng thái đó xảy ra.

Một trị số khác là Turbulence (theo Elzinga, 2006) cũng mô tả tính bất định của chuỗi, nhưng chịu chi phối của trình tự các trạng thái:

srh_seq%>%seqient()%>%hist(col="gold",xlab = "Entropy")

srh_seq%>%seqST()%>%head()

##   Turbulence
## 1   4.088201
## 2   7.705978
## 3   5.682061
## 4   6.108098
## 5   8.727920
## 6   5.625586

srh_seq%>%seqST()%>%hist(col="red",xlab = "Turbulence")

Ta có thể khảo sát tương quan giữa Entropy và Turbulence:

SRH%>%mutate(Entropy=seqient(srh_seq),
                 Turbulence=seqST(srh_seq))%>%
  ggplot(aes(x=Turbulence,y=Entropy))+
  geom_point(alpha=0.2)+
  geom_smooth(col="red",se=F)+
  theme_bw()

Và so sánh Entropy giữa 10 nhóm tuổi (năm sinh)

SRH%>%mutate(Entropy=seqient(srh_seq),
                 BYGroup=cut(SRH$birthy,b=10,
                        include.lowest = TRUE))%>%
  ggplot(aes(x=BYGroup,y=Entropy))+
  lvplot::geom_lv(col="black",
                  show.legend = F,
                  aes(fill=..LV..),
                  outlier.shape = 3)+
  scale_fill_brewer(palette="Spectral")+
  theme_bw()+coord_flip()

Biểu đồ này cho thấy Những người già và trung niên có sự bất định cao trong diễn tiến sức khỏe, biểu hiện bằng giá trị entropy cao, so với những người trẻ tuổi (sinh sau 1978).

Tiếp theo, chúng ta sẽ làm một phân tích xác suất đơn giản, nhằm khảo sát khả năng biến đổi từ một trạng thái sức khỏe xác định (thí dụ G1) ở thời điểm t sang một trạng thái khác ở thời điểm t+1 (bao gồm sự duy trì chính nó):

Phân tích này gồm 2 bước :

Đầu tiên ta sẽ tạo matrix xác suất cho từng cặp trạng thái:

pmat=seqtrate(srh_seq) %>% as.matrix()

pmat

##             [-> B1]      [-> B2]     [-> M]    [-> G1]   [-> G2]
## [B1 ->] 0.156549521 0.0447284345 0.42172524 0.03514377 0.3418530
## [B2 ->] 0.268292683 0.2682926829 0.26829268 0.04878049 0.1463415
## [M ->]  0.043119848 0.0038046925 0.38237159 0.06055802 0.5101458
## [G1 ->] 0.002486679 0.0000000000 0.03996448 0.46642984 0.4911190
## [G2 ->] 0.006657617 0.0007756448 0.10425958 0.15364230 0.7346649

colnames(pmat)=c("B1","B2","M","G1","G2")
rownames(pmat)=c("B1","B2","M","G1","G2")

Kết quả này rất thú vị, và hợp logic, ta có thể diễn giải nó một cách phong phú, thí dụ: Một người khỏe mạnh và lạc quan (G1) có nhiều cơ may tiếp tục khỏe mạnh vào năm sau, một người bệnh tật (B1) vẫn có đến 42% đạt mức trung bình (M) vào năm tiếp theo, nhưng rất ít cơ may trở về mức khỏe mạnh (3.5%). …

Kết quả này có thể được chuyển thành một biểu đồ mạng (network)

library(ggraph)

cg=data.frame(row=rownames(pmat)[row(pmat)], 
               col=colnames(pmat)[col(pmat)], 
               Prob=pmat%>%as.vector())

names(cg)=c("from","to","Prob")

library(igraph)

graph<-graph_from_data_frame(cg)

ggraph(graph,circular=F,layout = 'lgl')+
  geom_edge_fan(aes(col=Prob,label=round(Prob,2)),
                show.legend = T,
                arrow = arrow(length = unit(5, 'mm')),
                alpha=0.5)+
  geom_edge_loop(aes(col=Prob,label=round(Prob,2)),span=1,alpha=0.5)+
  geom_node_label(aes(label = name))+
  coord_fixed()+
  scale_edge_color_gradient(low="blue",high="red")+
  theme_graph()

Tiếp theo, ta có thể áp dụng phương pháp Clustering để phân cụm 200 cá thể chọn ngẫu nhiên trong mẫu dựa vào đặc tính của chuỗi trạng thái.

Đầu tiên, ta chọn ngẫu nhiên 200 cá thể, rồi dựng matrix giá trị distance theo phương pháp LCS (Longest Common Subsequence). Sau đó ta dùng phương pháp phân cụm của Ward.

Kết quả phân cụm có thể được trình bày bằng dendogram dưới 3 hình thức:

set.seed(123)

library(cluster)
library(factoextra)

srh_seq200<-seqdef(sample_n(SRH,200),5:15,
                  states=c("B1","B2","M","G1","G2"),
                  labels=c("NotW","Bad","Average","VeryW","Well"),
                  cpal=pals::brewer.rdylgn(5))

colnames(srh_seq200)=c(1999:2009)

dmat=seqdist(srh_seq200, method = "LCS",with.missing = T,norm=TRUE)

clusterward <- agnes(dmat, diss = TRUE, method = "ward")

fviz_dend(clusterward , k = 6, # Cut in 6 groups
          cex = 0.5, # label size
          k_colors = c("#24a3f2","#2edb9b","#243ff2","#f2a624","#a759f9", "#f95987"),
          color_labels_by_k = TRUE, # color labels by groups
          rect = TRUE, # Add rectangle around groups
          rect_border = c("#24a3f2","#2edb9b","#243ff2","#f2a624", "#a759f9", "#f95987"),
          rect_fill = TRUE)

fviz_dend(clusterward , k = 6, # Cut in 6 groups
          cex = 0.5, # label size
          k_colors = c("#24a3f2","#2edb9b","#243ff2","#f2a624","#a759f9", "#f95987"),
          color_labels_by_k = TRUE, # color labels by groups
          type = "circular")

library("igraph")

fviz_dend(clusterward, k = 6, # Cut in 6 groups
          k_colors = c("#24a3f2","#2edb9b","#243ff2","#f2a624","#a759f9", "#f95987"),         
          type = "phylogenic", repel = TRUE,
          phylo_layout = "layout.gem")

Lưu ý: Ở đây ta giả định phân dữ liệu thành 6 cụm.

Ta có thể giảm số lượng còn 3 cluster, những biểu đồ sau đây cho phép so sánh đặc tính trình tự chuỗi trạng thái giữa 3 cluster:

cluster3 <- cutree(clusterward, k = 3)
cluster3 <- factor(cluster3, labels = c("Cluster 1", "Cluster 2", "Cluster 3"))

seqfplot(srh_seq200, group = cluster3, pbarw = T)

seqmtplot(srh_seq200, group = cluster3)

Kết quả này cho thấy Cluster 1 gồm những người khỏe mạnh nhất (G1), diễn tiến sức khỏe của họ trong 10 năm dao động giữa G1 và G2 với ít nguy cơ giảm xuống mức B1 hay B2; Cluster 2 gồm những người ít lạc quan hơn nhưng có tình trạng sức khỏe rất tốt và ổn định với sữ duy trì mức G2; Cuối cùng, cluster 3 gồm những người có trạng thái sức khỏe không ổn định và dao động quanh mức trung bình /xấu, với nguy cơ cao xảy ra trạng thái B1, thậm chí B2

5 Kết luận

Bài thực hành ngắn hôm nay đã giới thiệu với các bạn một ý tưởng mới có thể áp dụng cho nghiên cứu lâm sàng/dịch tễ với thiết kế trường diễn/theo dõi kéo dài và kết quả là một biến rời rạc nhiều giá trị. Ý tưởng này độc đáo hơn so với nghiên cứu theo dõi biến định lượng theo thời gian hoặc survival analysis,với nhiều ứng dụng tiềm năng. Nhi cũng giới thiệu sơ lược về một công cụ trong R cho phép khảo sát chuỗi sự kiện/trạng thái. Lưu ý rằng trong bài chúng ta chỉ mới khảo sát khoảng 30% những tính năng mà package TraMineR cung cấp. Ngoài ra, một package khác có tên là PST còn cho phép dựng mô hình cây tiên đoán chuỗi sự kiện mà Nhi sẽ giới thiệu trong một dịp khác.

Chúc các bạn thực hành vui, và áp dụng được ý tưởng mới này vào nghiên cứu của riêng mình.