Phân tích mạng lưới tương quan

1 Khởi hành: Phân tích tương quan cổ điển

Chúng ta bắt đầu câu chuyện bằng một phân tích tương quan theo phương pháp cổ điển, sử dụng hệ số tương quan tuyến tính r (Pearson). Đây là một quy trình phổ biến trong bước thăm dò dữ liệu và những nghiên cứu mô tả.

Thí dụ minh họa trích từ một nghiên cứu cắt ngang, chủ đề lâm sàng nội khoa hô hấp, với 14 thông số chức năng hô hấp được khảo sát ở 2 nhóm bệnh nhân: Có và không có bệnh Xơ phổi mô kẽ.

library(tidyverse)

df = read.csv("RAW_DATA_ORG.csv",dec=",",sep=';')

neg_df = df%>%select(FEV1,Tiffneau,FVC,TLC,AV,RV,TLCO,KCO,TLNO,KNO,DLRatio,DmCO,DmNO,Vcap,ILD)%>%
  filter(.$ILD=="Negative")%>%.[,-15]

pos_df = df%>%select(FEV1,Tiffneau,FVC,TLC,AV,RV,TLCO,KCO,TLNO,KNO,DLRatio,DmCO,DmNO,Vcap,ILD)%>%
  filter(.$ILD=="Positive")%>%.[,-15]

neg_df=sapply(neg_df,as.numeric)
pos_df=sapply(pos_df,as.numeric)

14 thông số gồm: FEV1, tỉ số Tiffneau = FEV1/FVC khảo sát chức năng thông khí, FVC, TLC, RV, AV khảo sát Dung tích phổi, TLCO, KCO, TLNO, KNO, DLratio khảo sát khả năng trao đổi khí, DmCO và DmNO khảo sát dẫn suất của màng phế nang mao mạch, và Vcap khảo sát thể tích mao mạch phổi.

Nhóm Bệnh (Positive) có 147 bệnh nhân,

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

FEV1	Tiffneau	FVC	TLC	AV	RV	TLCO	KCO	TLNO	KNO	DLRatio	DmCO	DmNO	Vcap
154	312	160	106	183	8	264	306	35	336	275	19	239	299
125	155	155	133	187	48	123	83	305	135	325	2	222	119
60	175	68	56	83	46	193	330	304	333	210	315	160	216
9	119	6	112	43	167	31	138	119	78	92	82	302	146
50	296	42	22	41	16	331	49	162	195	369	247	92	39
27	182	25	42	89	66	57	105	152	72	178	128	348	124

Nhóm Negative gồm 229 người không có bệnh.

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

FEV1	Tiffneau	FVC	TLC	AV	RV	TLCO	KCO	TLNO	KNO	DLRatio	DmCO	DmNO	Vcap
105	150	126	177	213	140	274	251	44	346	343	36	255	297
151	269	164	187	173	114	184	182	18	313	342	24	244	181
108	331	97	144	164	127	79	63	222	79	305	204	49	176
45	44	71	123	170	126	25	30	204	47	351	282	127	72
68	276	63	84	123	89	90	106	299	283	358	11	231	97
16	30	30	179	116	181	71	94	233	185	335	280	125	118

Mục tiêu giả định là ta muốn khảo sát mối tương quan giữa 14 thông số chức năng này. Đây là một việc rất đơn giản như ta biết, chỉ cần áp dụng hàm cor trên data matrix có (k) biến số, ta sẽ có kết quả là một correlation matrix với kích thước (k x k)

cormat_pos = cor(pos_df, method="pearson")%>%round(.,3)
cormat_neg = cor(neg_df, method="pearson" )%>%round(.,3)

Ma trận tương quan 14 biến trên nhóm bệnh như sau:

cormat_pos %>% knitr::kable()

	FEV1	Tiffneau	FVC	TLC	AV	RV	TLCO	KCO	TLNO	KNO	DLRatio	DmCO	DmNO	Vcap
FEV1	1.000	0.074	0.939	0.748	0.650	0.150	-0.034	0.188	0.266	0.281	-0.008	-0.129	0.220	-0.015
Tiffneau	0.074	1.000	-0.252	-0.305	-0.214	-0.209	0.106	0.098	-0.052	0.115	0.072	-0.059	-0.022	0.028
FVC	0.939	-0.252	1.000	0.815	0.694	0.194	-0.071	0.144	0.266	0.242	-0.005	-0.106	0.214	-0.048
TLC	0.748	-0.305	0.815	1.000	0.725	0.713	-0.074	0.056	0.305	0.141	-0.012	-0.048	0.153	0.005
AV	0.650	-0.214	0.694	0.725	1.000	0.384	-0.027	-0.086	0.352	0.163	0.177	-0.048	0.149	0.034
RV	0.150	-0.209	0.194	0.713	0.384	1.000	-0.060	-0.063	0.181	-0.020	-0.033	0.044	-0.019	0.075
TLCO	-0.034	0.106	-0.071	-0.074	-0.027	-0.060	1.000	-0.159	-0.257	-0.012	0.207	0.070	-0.198	-0.081
KCO	0.188	0.098	0.144	0.056	-0.086	-0.063	-0.159	1.000	0.285	0.553	-0.555	-0.111	0.175	-0.027
TLNO	0.266	-0.052	0.266	0.305	0.352	0.181	-0.257	0.285	1.000	0.196	0.000	0.022	0.156	0.063
KNO	0.281	0.115	0.242	0.141	0.163	-0.020	-0.012	0.553	0.196	1.000	-0.036	-0.038	0.105	0.009
DLRatio	-0.008	0.072	-0.005	-0.012	0.177	-0.033	0.207	-0.555	0.000	-0.036	1.000	0.089	-0.078	-0.064
DmCO	-0.129	-0.059	-0.106	-0.048	-0.048	0.044	0.070	-0.111	0.022	-0.038	0.089	1.000	-0.444	0.322
DmNO	0.220	-0.022	0.214	0.153	0.149	-0.019	-0.198	0.175	0.156	0.105	-0.078	-0.444	1.000	-0.317
Vcap	-0.015	0.028	-0.048	0.005	0.034	0.075	-0.081	-0.027	0.063	0.009	-0.064	0.322	-0.317	1.000

Còn đây là kết quả trên nhóm Không có bệnh:

cormat_neg %>% knitr::kable()

	FEV1	Tiffneau	FVC	TLC	AV	RV	TLCO	KCO	TLNO	KNO	DLRatio	DmCO	DmNO	Vcap
FEV1	1.000	0.265	0.927	0.555	0.659	-0.153	0.498	0.281	-0.164	0.228	0.023	-0.157	0.136	0.065
Tiffneau	0.265	1.000	-0.059	-0.297	-0.124	-0.384	0.063	0.196	0.045	0.272	0.045	-0.063	-0.048	-0.014
FVC	0.927	-0.059	1.000	0.685	0.716	-0.039	0.498	0.246	-0.162	0.164	-0.005	-0.160	0.167	0.075
TLC	0.555	-0.297	0.685	1.000	0.762	0.659	0.348	0.003	-0.128	-0.054	0.136	-0.090	0.213	0.027
AV	0.659	-0.124	0.716	0.762	1.000	0.345	0.405	-0.093	-0.140	-0.035	0.310	-0.121	0.143	-0.002
RV	-0.153	-0.384	-0.039	0.659	0.345	1.000	-0.029	-0.271	0.003	-0.291	0.160	0.037	0.121	-0.016
TLCO	0.498	0.063	0.498	0.348	0.405	-0.029	1.000	0.467	-0.161	0.538	-0.202	-0.054	0.094	-0.051
KCO	0.281	0.196	0.246	0.003	-0.093	-0.271	0.467	1.000	0.102	0.727	-0.486	-0.022	0.046	-0.007
TLNO	-0.164	0.045	-0.162	-0.128	-0.140	0.003	-0.161	0.102	1.000	0.069	-0.012	-0.029	0.067	-0.082
KNO	0.228	0.272	0.164	-0.054	-0.035	-0.291	0.538	0.727	0.069	1.000	-0.064	-0.054	0.047	-0.077
DLRatio	0.023	0.045	-0.005	0.136	0.310	0.160	-0.202	-0.486	-0.012	-0.064	1.000	-0.076	-0.009	-0.087
DmCO	-0.157	-0.063	-0.160	-0.090	-0.121	0.037	-0.054	-0.022	-0.029	-0.054	-0.076	1.000	-0.462	0.173
DmNO	0.136	-0.048	0.167	0.213	0.143	0.121	0.094	0.046	0.067	0.047	-0.009	-0.462	1.000	-0.210
Vcap	0.065	-0.014	0.075	0.027	-0.002	-0.016	-0.051	-0.007	-0.082	-0.077	-0.087	0.173	-0.210	1.000

Trên thực tế, rất khó đọc kết quả trực tiếp từ correlation matrix, nên ta có thể chuyển những con số này thành hình ảnh (biểu đồ ma trận tương quan, correlograph) với package corrplot như sau:

Trước hết, Nhi viết 1 hàm để tạo ra correlation matrix với giá trị p_values thay vì hệ số tương quan r, kết quả này cho phép lọc bỏ những cặp tương quan không có ý nghĩa thống kê, thí dụ ở ngưỡng 0.05

cor.mtest <- function(mat, ...) {
  mat <- as.matrix(mat)
  n <- ncol(mat)
  p.mat<- matrix(NA, n, n)
  diag(p.mat) <- 0
  for (i in 1:(n - 1)) {
    for (j in (i + 1):n) {
      tmp <- cor.test(mat[, i], mat[, j], ...)
      p.mat[i, j] <- p.mat[j, i] <- tmp$p.value
    }
  }
  colnames(p.mat) <- rownames(p.mat) <- colnames(mat)
  p.mat
}

corsig_pos = cor.mtest(pos_df)%>%round(.,3)
corsig_neg = cor.mtest(neg_df)%>%round(.,3)

Biểu đồ tương quan 14 biến ở nhóm Bệnh nhân:

library(corrplot)

corrplot(cormat_pos, 
         col= pals::coolwarm(100),
         method="pie",
         type="upper",
         order="hclust",
         p.mat = corsig_pos, 
         sig.level = 0.05,
         insig = "blank",
         tl.col="black", tl.srt=45,tl.cex=0.8,
         diag=FALSE
         )

Biểu đồ của nhóm người Bình thường:

corrplot(cormat_neg, 
         col= pals::coolwarm(100),
         method="pie",
         type="upper",
         order="hclust",
         p.mat = corsig_neg, 
         sig.level = 0.05,
         insig = "blank",
         tl.col="black", tl.srt=45,tl.cex=0.8,
         diag=FALSE
)

Dựa trên 2 biểu đồ này, ta có thể cảm nhận trực quan rằng có sự thay đổi về cấu trúc các cặp tương quan giữa 2 trạng thái: Có bệnh và không có bệnh. Nguyên nhân vì bệnh xơ phổi gây tác động khác nhau lên mỗi chức năng sinh lý hô hấp, thí dụ nó làm giảm diện tích khả dụng của màng phế nang mao mạch, tăng bề dày của màng này, giảm khả năng trao đổi khí … trong khi chức năng thông khí có thể vẫn còn bình thường. Như vậy, tương quan giữa 2 biến (đại lượng), ngay cả ở cấp độ quần thể hay mẫu, thực ra là một đặc tính “Động”. Một nghiên cứu cắt ngang chỉ có thể ghi nhận được đặc tính này ở một trạng thái, tại thời điểm nào đó, nhưng kết quả này có thể thay đổi ở một thời điểm / điều kiện khác.

Theo cách làm cổ điển, chúng ta chỉ dừng lại ở đây, hệ số tương quan r và p_value chỉ cho phép chúng ta suy diễn thống kê về độ mạnh, chiều hướng và ý nghĩa thống kê của mối tương quan giữa từng CẶP biến, nhưng không cho phép nhận định về quan hệ phức tạp, chồng chéo giữa TẤT CẢ các biến, dù có correlation matrix trong tay.

Trong bài này, Nhi sẽ giới thiệu với các bạn phương pháp Network analysis cho phép ta đi xa hơn và khai thác thêm nhiều thông tin từ correlation matrix.

2 Từ Correlation matrix đến correlation network

Phân tích mạng lưới (network analysis), hay mô hình mạng (graphical model) là một nhánh của khoa học Thống kê, nhằm giải quyết những vấn đề về mối quan hệ, sự tương tác giữa các cá thể/phần tử trong một tập hợp/quần thể

Khái niệm “mối liên hệ” trong đời sống rất đa dạng, phong phú: việc kết bạn với nhau trên Facebook, đồng tác giả một bài báo, liên lạc điện thoại ,giao dịch ngân hàng, cuộc hôn nhân giữa thành viên 2 gia tộc,một cây cầu nối 2 tỉnh/thành phố, 2 biến cùng xuất hiện trong một mô hình…

Trong trường hợp phân tích tương quan, mối liên kết giữa 2 biến được xác lập nhờ vào giá trị của r và p_value. Ta sẽ mượn các công cụ và lý thuyết Network analysis cho mục tiêu phân tích tương quan. Nhưng trước hết, ta phải tạo ra correlation network. Bằng cách nào ?

Thực ra, một network có thể được dựng từ 2 cấu trúc dữ liệu cơ bản là “adjacency matrix”, hoặc “edges list”. Correlation matrix chính là một adjacency matrix.

Trong R, có 2 nhóm công cụ khác nhau cho Network analysis, đó là packages statnet và igraph. Mỗi nhóm sử dụng dữ liệu network có cấu trúc khác nhau (network và graph). Ta có thể hoán chuyển giữa graph và network bằng package intergraph.

Từ correlation matrix, Nhi sẽ dùng igraph để chuyển adjacency matrix thành graph. Đồng thời, ta cũng sẽ tạo thêm 1 vài thuộc tính cho edges (các mối nối liên kết trong mạng lưới):Direction = hướng quan hệ, nhận giá trị Proportional nếu r>0, Inverse nếu r<0; và Strength: độ mạnh tương quan: Strong nếu abs(r) >0.5 và Weak nếu <0.5.

Làm thử trên nhóm bệnh nhân;

m=cormat_pos
diag(m)<-0
cor_pos_df = data.frame(row=rownames(m)[row(m)[upper.tri(m)]], 
            col=colnames(m)[col(m)[upper.tri(m)]], 
            corr=m[upper.tri(m)])

m=corsig_pos
diag(m)<-0
sig_pos_df = data.frame(row=rownames(m)[row(m)[upper.tri(m)]], 
                        col=colnames(m)[col(m)[upper.tri(m)]], 
                        corr=m[upper.tri(m)])

cor_pos_df$sig = sig_pos_df$corr
rm(sig_pos_df)

names(cor_pos_df)=c("from","to","r","sig")

reduced_pos = cor_pos_df%>%filter(sig<0.05)
reduced_pos<-reduced_pos%>%mutate(Strength = if_else(abs(reduced_pos$r)>0.5,"Strong","Weak"),
                     Direction = if_else(reduced_pos$r>0,"Pro","Inv"))

library(igraph)

graph_pos<-graph_from_data_frame(reduced_pos,directed = F)

Đây là kết quả của package igraph (lưu ý: gần đây có thêm package ggraph cho phép vẽ một dữ liệu igraph thành biểu đồ mạng sử dụng ngữ pháp đồ họa của Wilkinsonnhuè ggplot2, nhưng Nhi không bàn đến nó trong bài này)

plot(graph_pos)

Sau đó, ta dùng intergaph để chuyển graph object thành network object, nguyên nhân Nhi chọn làm việc trên network chứ không phải graph, vì định dạng dữ liệu network tương thích với hệ sinh thái statnet, cho phép các bạn đi xa hơn sau này với rất nhiều phương pháp, mô hình mà hệ sinh thái này cung cấp.

library(intergraph)
library(network)
library(statnet)

pos_net <- asNetwork(graph_pos)

plot(pos_net,displaylabel=T,label.cex=0.8)

Đây là hình ảnh mạng lưới tương quan giữa 14 biến trong nhóm Bệnh nhân, Nhi sẽ đưa thêm vào thuộc tính của vertex (là các nút trên mạng lưới, còn gọi là actor), biến thuộc tính đưa vào là Physio, nhằm phân loại nhóm chức năng sinh lý cho các thông số: Vent = Ventilation, Vol = lung volume, Diff =Diffusion, Memb = membrane conductance, Cap = lung capillary volume:

set.vertex.attribute(pos_net, "Physio",
                     c("vent","vent","Vol","Vol","Vol","Vol",
                       "Diff","Diff","Memb","Memb","Memb",    
                       "Memb","Diff","Cap"))

pos_net %v% "Physio"

##  [1] "vent" "vent" "Vol"  "Vol"  "Vol"  "Vol"  "Diff" "Diff" "Memb" "Memb"
## [11] "Memb" "Memb" "Diff" "Cap"

Ta làm tương tự cho nhóm Negative:

# negative

m=cormat_neg
diag(m)<-0
cor_neg_df = data.frame(row=rownames(m)[row(m)[upper.tri(m)]], 
                        col=colnames(m)[col(m)[upper.tri(m)]], 
                        corr=m[upper.tri(m)])

m=corsig_neg
diag(m)<-0
sig_neg_df = data.frame(row=rownames(m)[row(m)[upper.tri(m)]], 
                        col=colnames(m)[col(m)[upper.tri(m)]], 
                        corr=m[upper.tri(m)])

cor_neg_df$sig = sig_neg_df$corr
rm(sig_neg_df)

names(cor_neg_df)=c("from","to","r","sig")

reduced_neg = cor_neg_df%>%filter(sig<0.05)
reduced_neg<-reduced_neg%>%mutate(Strength = if_else(abs(reduced_neg$r)>0.5,"Strong","Weak"),
                                  Direction = if_else(reduced_neg$r>0,"Pro","Inv"))

graph_neg<-graph_from_data_frame(reduced_neg,directed = F)

plot(graph_neg)

neg_net <- asNetwork(graph_neg)

plot(neg_net,displaylabel=T,label.cex=0.8)

set.vertex.attribute(pos_net, "Physio",
                     c("vent","vent","Vol","Vol","Vol","Vol",
                       "Diff","Diff","Memb","Memb","Memb",    
                       "Memb","Diff","Cap"))

3 Mô tả cấu trúc mạng tương quan

Thực ra, network không chỉ là một hình vẽ trình bày ý tưởng, khái niệm về sự tương quan, nhưng là một mô hình chính xác và toán học về kích thước, cấu trúc phân bố của các nút (node, vertex, actor) và mạng liên kết giữa chúng (edges, arcs).

Sau đây, Nhi sẽ làm một vài thống kê mô tả để khảo sát cấu trúc 2 network mà ta vừa tạo ra:

1. Kích thước tính bằng số nodes hay vertices:

Cả 2 network đều có kích thước như nhau là 14 vertices

network.size(pos_net)

## [1] 14

network.size(neg_net)

## [1] 14

2. Mô tả cấu trúc bằng hàm summary:

Nhóm Positive có 36 edges, 4 edges attributes là r, sig, strength và direction và 2 vertex atribute là vertex.names và Physio

summary(pos_net,print.adj=FALSE)

## Network attributes:
##   vertices = 14
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##  total edges = 36 
##    missing edges = 0 
##    non-missing edges = 36 
##  density = 0.3956044 
## 
## Vertex attributes:
## 
##  Physio:
##    character valued attribute
##    attribute summary:
##  Cap Diff Memb vent  Vol 
##    1    3    4    2    4 
##   vertex.names:
##    character valued attribute
##    14 valid vertex names
## 
## Edge attributes:
## 
##  Direction:
##    character valued attribute
##    attribute summary:
## Inv Pro 
##   9  27 
## 
##  r:
##    numeric valued attribute
##    attribute summary:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.55500  0.07275  0.21700  0.21400  0.36000  0.93900 
## 
##  sig:
##    numeric valued attribute
##    attribute summary:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.001000 0.007667 0.011250 0.049000 
## 
##  Strength:
##    character valued attribute
##    attribute summary:
## Strong   Weak 
##      9     27

Nhóm negative có nhiều edges hơn (45)

summary(neg_net,print.adj=FALSE)

## Network attributes:
##   vertices = 14
##   directed = FALSE
##   hyper = FALSE
##   loops = FALSE
##   multiple = FALSE
##   bipartite = FALSE
##  total edges = 45 
##    missing edges = 0 
##    non-missing edges = 45 
##  density = 0.4945055 
## 
## Vertex attributes:
##   vertex.names:
##    character valued attribute
##    14 valid vertex names
## 
## Edge attributes:
## 
##  Direction:
##    character valued attribute
##    attribute summary:
## Inv Pro 
##  15  30 
## 
##  r:
##    numeric valued attribute
##    attribute summary:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## -0.4860 -0.1600  0.1960  0.1818  0.4670  0.9270 
## 
##  sig:
##    numeric valued attribute
##    attribute summary:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.000000 0.000000 0.006533 0.013000 0.040000 
## 
##  Strength:
##    character valued attribute
##    attribute summary:
## Strong   Weak 
##      9     36

3. Mật độ của mạng lưới:

Mật độ của mạng tương quan cho phép hình dung về khả năng liên kết giữa các phần tử trong tập hợp (ở đây là các xét nghiệm, thông số chức năng sinh lý).

# density

gden(pos_net)

## [1] 0.3956044

gden(neg_net)

## [1] 0.4945055

Mạng tương quan của nhóm negative có mật độ cao hơn so với mạng của nhóm Positive, phù hợp với khác biệt về số cặp tương quan ta thấy ở trên. Nếu giả định nhóm không có bệnh thể hiện trạng thái sinh lý bình thường, thì kết quả này cho thấy bệnh lý đã làm đứt gãy nhiều cặp tương quan sinh lý, thí dụ giữa chức năng thông khí và chức năng trao đổi khí.

4. Đường kính của mạng lưới

Cả 2 mạng lưới đều có diameter=3, đây là con đường nối từ đầu này sang đầu kia (2 vertices) của mạng lưới

# Diameter

diameter= function(network) {
  component.largest(network,result="graph")%>%geodist()->gd
  max(gd$gdist)
  }

diameter(pos_net)

## Node 1, Reach 14, Total 14
## Node 2, Reach 14, Total 28
## Node 3, Reach 14, Total 42
## Node 4, Reach 14, Total 56
## Node 5, Reach 14, Total 70
## Node 6, Reach 14, Total 84
## Node 7, Reach 14, Total 98
## Node 8, Reach 14, Total 112
## Node 9, Reach 14, Total 126
## Node 10, Reach 14, Total 140
## Node 11, Reach 14, Total 154
## Node 12, Reach 14, Total 168
## Node 13, Reach 14, Total 182
## Node 14, Reach 14, Total 196

## [1] 3

diameter(neg_net)

## Node 1, Reach 14, Total 14
## Node 2, Reach 14, Total 28
## Node 3, Reach 14, Total 42
## Node 4, Reach 14, Total 56
## Node 5, Reach 14, Total 70
## Node 6, Reach 14, Total 84
## Node 7, Reach 14, Total 98
## Node 8, Reach 14, Total 112
## Node 9, Reach 14, Total 126
## Node 10, Reach 14, Total 140
## Node 11, Reach 14, Total 154
## Node 12, Reach 14, Total 168
## Node 13, Reach 14, Total 182
## Node 14, Reach 14, Total 196

## [1] 3

5. Chỉ số Transitivity: tỉ lệ giữa số liên kết bộ ba khép kín trên tổng số liên kết bộ ba hở và kín

gtrans(neg_net)

## [1] 0.6

gtrans(pos_net)

## [1] 0.5698324

Mạng của nhóm negative có tỉ số transitivity cao hơn, cho thấy nó có nhiều liên kết bộ ba hơn so với mạng tương quan của nhóm bệnh nhân

Lúc này, cảm nhận mơ hồ của Nhi bên trên về sự khác biệt giữa 2 mạng lưới tương quan của nhóm Bệnh nhân và người bình thường đã được chứng thựcmột cách định lượng: Ở nhóm Negative (không có bệnh), có nhiều mối liên hệ hơn giữa các biến, và một vài liên kết trong số này đã bị đứt gãy (biến mất) do tác động của bệnh lý xơ phổi.

Ta có thể khảo sát một cách trực quan khác biệt này bằng cách so sánh hình ảnh 2 mạng lưới tương quan: Hình ảnh này dễ hiểu và rõ ràng hơn nhiều so với correlogram : Ta có thể thấy những liên kết đã bị mất đi bao gồm: Giữa RV và Dlratio, giữa Tiffneau và KNO, KCO, giữa TLC,FVC,FEV1 và DmCO, … hình ảnh này phù hợp với thực tế, khi đa số các bệnh nhân có bệnh lý hô hấp đều có sự bảo tồn chức năng thông khí trên xét nghiệm hô hấp ký: kết quả xét nghiệm này vẫn bình thường, nhưng thay đổi bệnh lýđược thể hiện ra ở các xét nghiệm khác như đo dung tích phổi và/hoặc TLCO.

# Comparative plots

coords<-plot(pos_net,label=pos_net%v%"vertex.names",label.cex=.5,pad=1)

par(mar=c(0,0,2,0),mfrow=c(1,2))

gplot(pos_net,gmode="graph",
      mode="kamadakawai",
      label=pos_net%v%"vertex.names",
      label.cex=.5,pad=1,
      vertex.col = "red",
      vertex.cex=0.8,
      edge.col="grey50",
      coord = coords)
title("Patients")

gplot(neg_net,gmode="graph",
      mode="kamadakawai",
      label=pos_net%v%"vertex.names",
      label.cex=.5,pad=1,
      vertex.col = "blue",
      edge.col="grey50",
      vertex.cex=0.8,
      coord = coords)
title("Normal")

par(mfrow=c(1,1))

4 Phần tử (biến) nào có vai trò quan trọng nhất ?

Có nhiều ý kiến khác nhau về tầm quan trọng của các xét nghiệm sinh lý hô hấp. Thí dụ, hiện nay hô hấp ký (FEV1, FVC) là xét nghiệm thông dụng nhất, trong khi có quan điểm cho rằng xét nghiệm DLCO nhạy hơn để phát hiện sớm tổn thương chức năng của phổi và cho phép tiên lượng nguy cơ tử vong…Mặt khác, nếu tồn tại mối tương quan giữa 2 thông số xét nghiệm, sẽ dẫn đến suy nghĩ là chỉ cần sử dụng 1 thông số là đủ khai thác thông tin ?

Phần tiếp theo, Nhi sẽ dùng phân tích mạng lưới để nhận diện những biến có vai trò trung tâm, đóng góp quan trọng vào việc kiến tạo mạng liên kết với những biến còn lại, từ đó trả lời câu hỏi: chức năng hô hấp nào là quan trọng nhất trong toàn bộ các xét nghiệm đang được nghiên cứu ?

Một phần tử có vai trò trung tâm khi nó có nhiều liên hệ nhất với các phần tử khác trong mạng lưới. Có nhiều chỉ số để ước lượng tính Trung tâm này, như: Degree, Betweeness, Closeness, Eigen vector, Bonacich power, Information, Harary graph:

# Centrality measures combined

cent_df = data_frame(Vertex = c(neg_net %v% "vertex.names",
                                pos_net %v% "vertex.names"),
                     Physio = c(rep(c("vent","vent","Vol","Vol","Vol","Vol",
                                "Diff","Diff","Memb","Memb","Memb",    
                                "Memb","Diff","Cap"),2)),
                     
                     Status = c(rep("Normal",14),rep("Patients",14)),
                       
                     Degree = c(degree(pos_net, gmode="graph"),
                                degree(neg_net, gmode="graph")),
                     
                     Closeness = c(closeness(pos_net, gmode="graph"),
                                   closeness(neg_net, gmode="graph")),
                     
                     Beetweeness = c(betweenness(pos_net, gmode="graph"),
                                     betweenness(neg_net, gmode="graph")),
                     
                     Eigen = c(evcent(pos_net, gmode="graph"),
                               evcent(neg_net, gmode="graph")),
                     
                     BonPow = c(bonpow(pos_net, gmode="graph"),
                                bonpow(neg_net, gmode="graph")),
                     
                     InfoCent = c(infocent(pos_net, gmode="graph"),
                                  infocent(neg_net, gmode="graph"))
                     )

Có thể thấy, ở từng trạng thái khác nhau: Bệnh lý/Bình thường, mỗi biến có vai trò khác nhau trong việc kiến tạo ra mạng lưới quan hệ với các biến còn lại.

Một số biến có vai trò trung tâm nổi bật, như FEV1, FVC (xét nghiệm hô hấp ký). Ta cũng có thể so sánh vai trò trung tâm giữa các nhóm xét nghiệm chức năng:

cent_df%>%gather(Degree:InfoCent,key="Score",value="Centrality")%>%
  ggplot(aes(x=Vertex,y=Centrality,fill=Status))+
  geom_point(stat="identity",size=2,shape=21,col="black")+
  theme_bw()+coord_flip()+
  facet_wrap(~Score,ncol=3,scales="free")

cent_df%>%gather(Degree:InfoCent,key="Score",value="Centrality")%>%
  ggplot(aes(x=Physio,y=Centrality,fill=Status))+
  geom_boxplot(alpha=0.7)+
  theme_bw()+coord_flip()+
  facet_wrap(~Score,ncol=2,scales="free")

cent_df%>%gather(Degree:InfoCent,key="Score",value="Centrality")%>%
  ggplot(aes(x=Status,y=Centrality,fill=Physio))+
  geom_boxplot(alpha=0.7)+
  theme_bw()+coord_flip()+
  facet_wrap(~Score,ncol=2,scales="free")

Ta còn có thể phát hiện ra phần tử then chốt “cut-point” , có nghĩa là khi loại bỏ phần tử này, cấu trúc mạng lưới sẽ bị phân rã, chia cắt thành nhiều cụm. Trong trường hợp này, có 1 cutpoint chính là biến DmNO. Nếu loại bỏ DmNO, 2 biến DmCO và Vcap sẽ bị tách rời ra khỏi mạng lưới.

# Cutpoint

cp_pos <- cutpoints(pos_net,mode="graph",
                   return.indicator=TRUE)

gplot(pos_net,gmode="graph",
      label=pos_net%v%"vertex.names",
      label.cex=.8,pad=1,
      vertex.col=cp_pos+1,coord=coords,
      jitter=FALSE,displaylabels=TRUE)

5 Yếu tố mỹ thuật của mạng tương quan

Trước khi kết thúc, chúng ta sẽ làm một số thử nghiệm thay đổi tùy chỉnh mỹ thuật cho biểu đồ mạng tương quan:

Có nhiều hình thức trình bày mạng lưới tương quan, sau đây là 6 kiểu thông dụng nhất:

op <- par(mar=c(0,0,2,0),mfrow=c(2,3))
gplot(pos_net,gmode="graph",mode="circle",edge.col="grey",vertex.col = "red",
      vertex.cex=1,main="Circle")
gplot(pos_net,gmode="graph",mode="random",edge.col="grey",vertex.col = "gold",
      vertex.cex=1,main="Random layout")
gplot(pos_net,gmode="graph",mode="fruchtermanreingold",edge.col="grey",vertex.col = 3,
      vertex.cex=1,main="Fruchterman-Reingold")
gplot(pos_net,gmode="graph",mode="spring",edge.col="grey",vertex.col = 4,
      vertex.cex=1,main="Spring")
gplot(pos_net,gmode="graph",mode="eigen",edge.col="grey",vertex.col = 6,
      vertex.cex=1,main="Eigenstructure")
gplot(pos_net,gmode="graph",mode="kamadakawai",edge.col="grey",vertex.col = "black",
      vertex.cex=1,main="Kamadakawai")

par(op)

Ta có thể thấy rằng 3 kiểu trình bày: mạng lưới vòng (Circle), Kamadakawai và Fruchterman-Reingold là rõ ràng và đẹp nhất.

Ta có thể tùy chỉnh kích thước và màu vertices theo attribute của nó: thí dụ 3 chỉ số Dgree, Closeness và Betweeness:

deg <- degree(pos_net,gmode="graph")
cls <- closeness(pos_net,gmode="graph")
bet <- betweenness(pos_net,gmode="graph")

gplot(pos_net, 
      vertex.cex=log(deg), 
      label.cex=0.8,
      label=pos_net%v%"vertex.names",
      vertex.col=rgb(deg/max(deg),0.2,1-deg/max(deg)),
      gmode="graph",
      coord=coords)

gplot(pos_net, 
      vertex.cex=3*cls, 
      label.cex=0.8,
      label=pos_net%v%"vertex.names",
      vertex.col=rgb(cls/max(cls),0.3,1.2-cls/max(cls)),
      gmode="graph",
      coord=coords)

gplot(pos_net, 
      vertex.cex=sqrt(bet+1), 
      label.cex=0.8,
      label=pos_net%v%"vertex.names",
      vertex.col=rgb(bet/max(bet),1-bet/max(bet),0.1),
      gmode="graph",
      coord=coords)

Tương tự, ta có thể tùy chỉnh màu sắc, dán nhãn cho các edges:

edge_dir <- as.factor(pos_net%e%"Direction")
edge_dir_pal <- c("blue","red")

plot(pos_net,label=pos_net%v%"vertex.names",
     displaylabels=T,
     label.cex=0.8,
     vertex.cex=3*cls,
     vertex.col=rgb(cls/max(cls),0.3,1.2-cls/max(cls)),
     edge.col=edge_dir_pal[edge_dir],edge.lwd=1,
     edge.label=edge_dir,edge.label.cex=0.5,
     edge.label.col=edge_dir_pal[edge_dir],
     mode="kamadakawai")

edge_val <- as.factor(pos_net%e%"r")
edge_val_pal <- pals::coolwarm(length(edge_val))

plot(pos_net,label=pos_net%v%"vertex.names",
     displaylabels=T,
     vertex.cex=2*cls,
     vertex.col=rgb(cls/max(cls),0.3,1.2-cls/max(cls)),
     label.cex=0.5,
     edge.col=edge_val_pal[edge_val],edge.lwd=1,
     edge.label=edge_val,
     edge.label.cex=0.5,
     edge.label.col=edge_val_pal[edge_val],
     mode="kamadakawai")

edge_val_pos <- as.factor(pos_net%e%"r")
edge_val_pos_pal <- pals::coolwarm(length(edge_val_pos))

edge_val_neg <- as.factor(neg_net%e%"r")
edge_val_neg_pal <- pals::coolwarm(length(edge_val_neg))

deg_pos <- degree(pos_net,gmode="graph")
deg_neg <- degree(neg_net,gmode="graph")

par(mar=c(0,0,2,0),mfrow=c(1,2))

plot(pos_net,gmode="graph",
      mode="kamadakawai",
      label=pos_net%v%"vertex.names",
      label.cex=.5,pad=1,
      vertex.col=rgb(deg_pos/max(deg_pos),0.2,1-deg_pos/max(deg_pos)),
      vertex.cex=log(deg_pos)*2, 
      edge.col=edge_val_pos_pal[edge_val_pos],edge.lwd=1,
      edge.label=edge_val_pos,
      edge.label.cex=0.5,
      edge.label.col=edge_val_pos_pal[edge_val_pos],
      coord = coords)
title("Patients")

plot(neg_net,gmode="graph",
      mode="kamadakawai",
      label=pos_net%v%"vertex.names",
      label.cex=.5,pad=1,
      vertex.col=rgb(deg_neg/max(deg_neg),0.2,1-deg_neg/max(deg_neg)),
      vertex.cex=log(deg_neg)*2, 
      edge.col=edge_val_neg_pal[edge_val_neg],edge.lwd=1,
      edge.label=edge_val_neg,
      edge.label.cex=0.5,
      edge.label.col=edge_val_neg_pal[edge_val_neg],
      coord = coords)
title("Normal")

par(mfrow=c(1,1))

6 Tổng kết

Bài thực hành đến đây là hết. Trong bài này, Nhi đã kết hợp giữa phân tích tương quan cổ điển và một phương pháp thống kê hiện đại là Network analysis. Kết quả là ta có thể khai thác nhiều thông tin hơn từ ma trận tương quan, khảo sát một cách định lượng cấu trúc của mối quan hệ đa chiều giữa các biến, so sánh 2 mạng lưới tương quan dựa vào các chỉ số về mật độ và vai trò của mỗi biến trong việc kiến tạo nên mạng lưới quan hệ giữa chúng.