Mạng mô hình BAYES

1 Giới thiệu

Thuyên tắc phổi là tình trạng một hay nhiều động mạch trong phổi bị nghẽn bởi một cục máu đông di chuyển từ một phần khác của cơ thể, như chi dưới (huyết khối tĩnh mạch sâu (DVT)). Thuyên tắc phổi có thể gây tử vong nếu bệnh nhân không được chẩn đoán và điều trị kịp thời. Chụp động mạch phổi với thuốc cản quang là xét nghiệm chính xác nhất để chẩn đoán tắc mạch phổi, nhưng có tính xâm lấn và nguy cơ cao nên chỉ được thực hiện khi các xét nghiệm khác không cung cấp đủ bằng chứng.Lượng D-dimer trong máu tăng có thể gợi ý về sự hiện diện của cục máu đông, tuy nhiên D-dimer cũng có thể tăng bởi các yếu tố khác, như phẫu thuật, mang thai.

Đầu tiên, chúng ta trích xuất bộ dữ liệu PE từ package HydeNet, cần ghi chú đây chỉ là một dữ liệu mô phỏng với mục tiêu minh họa chứ không dựa trên nghiên cứu có thực.

library(tidyverse)
library(HydeNet)

data(PE)

PE%>%head()

str(PE)

## 'data.frame':    10000 obs. of  7 variables:
##  $ wells   : num  3 2 3 3 2 6 4 6 4 2 ...
##  $ pregnant: Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ pe      : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 2 2 1 2 1 ...
##  $ angio   : Factor w/ 2 levels "Negative","Positive": 1 2 1 1 1 2 2 1 1 1 ...
##  $ d.dimer : num  206 180 197 192 281 318 303 234 229 235 ...
##  $ treat   : Factor w/ 2 levels "No","Yes": 1 2 1 2 2 2 2 1 1 1 ...
##  $ death   : Factor w/ 2 levels "No","Yes": 1 1 1 1 1 1 2 1 2 1 ...

Như ta thấy, PE là 1 dataframe với 10000 trường hợp (bệnh nhân), với 7 biến lần lượtco& nội dung:

1. wells là 1 biến số rời rạc, trình bày điểm lâm sàng theo tiêu chí Wells.

2. pregnant là 1 biến nhị phân, chỉ tình trạng Có/không mang thai.

3. pe là 1 biến nhị phân, cho biết chẩn đoán của bệnh thuyến tắc phổi.

4. angio là 1 biến nhị phân, trình bày kết quả xét nghiệm xạ hình động mạch phổi với thuốc cản quang, với 2 giá trị Positive/Negative.

5. d.dimer là biến số liên tục, trình bày kết quả xét nghiệm định lượng D.dimer trong máu.

6. treat là 1 biến nhị phân, là quyết định điều trị : Có/không.

7. death là 1 biến nhị phân, cho biết kết cục của bệnh nhân: No=sống sót, Yes=Tử vong

2 Bước 1: Thiết kế cấu trúc mạng mô hình

Với dữ liệu hiện thời, rất nhiều hướng nghiên cứu có thể đặt ra, với các giả thuyết/câu hỏi khác nhau. Thí dụ: Ta có thể dựng một mô hình nhằm tiên lượng nguy cơ thuyên tắc phổi (kết quả = biến pe) dựa vào các thông tin lâm sàng (điểm wells, kết quả xét nghiệm d.dimer, angiography); ta cũng có thể nhắm đến kết cục “tử vong” của bệnh nhân và khảo sát mối liên hệ giữa kết cục này và chẩn đoán (pe) và chỉ định điều trị (treat);Cũng có thể ta muốn khảo sát sự ảnh hưởng của kết quả các xét nghiệm và điểm wells với quyết định điều trị của bác sĩ…

Tuy nhiên, trong bài thực hành hôm nay Nhi sẽ cùng các bạn đặt ra một cách tiếp cận vấn đề toàn diện và phức tạp hơn, với mục tiêu/ý tưởng:

Dựng mô hình cho phép khảo sát được toàn bộ (mạng lưới) những mối liên hệ / ảnh hưởng khả dĩ giữa các yếu tố: Bệnh lý (pe), kết quả xét nghiệm cận lâm sàng (d.dimer, angio), cảm nhận chủ quan (wells), quyết định điều trị (treat), và kết cục lâm sàng (nguy cơ tử vong) của bệnh nhân. Kết quả mô hình bao gồm chiều hướng/kích thước của mỗi hiệu ứng /liên hệ.

Mục tiêu này có thể hiểu như: tôi muốn dùng phân tích thống kê để suy diễn nhân/quả: giữa các yếu tố, và tôi muốn kể một câu chuyện có ý nghĩa dựa trên dữ liệu này.

Khác với những phương pháp phân tích mô hình đa biến quy ước - vốn thường quá chú trọng về thống kê, thí dụ chọn lọc mô hình / biến số dựa vào các chỉ số BIC, AIC (Stepwise) hoặc BMA, các phương pháp rút gọn mô hình như Lasso…, trong bài toán này chúng ta sẽ ưu tiên xây dựng mạng lưới tương tác giữa các yếu tố một cách chủ quan, dựa vào kiến thức y học, kinh nghiệm và quy luật hợp lý khác.

Ta có thể bắt đầu việc này như sau:

Trước hết, mạng lưới mô hình của chúng ta có 7 nút , trong đó 2 nút quan trọng nhất là Chẩn đoán (pe) và Điều trị (treat), nút tận cùng là kết cục: sống sót hay tử vong (death).

Ở thời điểm nhập viện, chưa có xét nghiệm nào được làm, bệnh nhân được khai thác bệnh sử và khám bởi bác sĩ, do đó điểm khởi đầu của mạng lưới chính là nút wells. Như ta biết, Wells là một thang điểm lâm sàng tiên lượng nguy cơ thuyên tắc phổi, dựa vào các tiêu chí như bệnh sử (bất động kéo dài, phẫu thuật, bệnh lý tĩnh mạch, điều trị ung thư), triệu chứng (dấu hiệu huyết khối tĩnh mạch, ho ra máu, nhịp tim nhanh…). Ta có thể nói chẩn đoán lâm sàng (pe) phụ thuộc thang điểm Wells, từ đó vạch ra liên kết đầu tiên trong mạng lưới: pe | wells. Ta có thể trình bày quan hệ điều kiện này bằng một mô hình logistic đơn biến.

Tiếp theo, ta xét nút D.dimer như biến kết quả. Có thể nói giá trị d.dimer phụ thuộc vào tình trạng có thuyên tắc phổi hay không, do đó có quan hệ giữa d.dimer và pe; tuy nhiên kết quả xét nghiệm d.dimer có thể bị sai lệch khi có thai (pregnant = Yes, với nguy cơ dương tính giả, do đó cả pe và pregnant cùng tham gia vào mô hình tiên lượng d.dimer.

kết quả xét nghiệm Angiography chỉ phụ thuộc vào sự hiện diện thuyên tắc phổi (pe) hay không, nên chỉ có 1 liên kết từ pe đến angio.

Tiếp theo, ta xét đến nút: chỉ định điều trị (treat): Bệnh nhân chỉ được điều trị khi có xét nghiệm d.dimer và/hoặc angio dương tính, xác định chẩn đoán; do đó ta có thể dùng một mô hình logistic khác với 2 biến là angio và d.dimer.

Cuối cùng, nút tận cùng của mạng lưới là kết cục tử vong của bệnh nhân, đây là kết quả cửa sự tương tác giữa 2 điều kiện: bệnh lý thuyên tắc phổi (pe) và chỉ định điều trị (treat). Thí dụ: ở bệnh nhân có thuyên tắc phổi , nguy cơ tử vong có thể thấp hơn nếu được điều trị kịp thời, ngược lại nguy cơ tử vong cao nếu không được điều trị.

Sau khi phác thảo mạng mô hình trên giấy, Nhi dùng package HydeNet để thiết kế cấu trúc của mạng lưới mô hình và vẽ sơ đồ :

net <- HydeNetwork(~ wells
                   + pe | wells
                   + d.dimer | pregnant*pe
                   + angio | pe
                   + treat | d.dimer*angio
                   + death | pe*treat)

HydePlotOptions(variable=list(fillcolor = "gold"))

plot(net)

Bên trong object net này có chứa sẵn một ma trận DAG (Directed Acyclic graph), Nhi muốn bật mí cho các bạn, vì từ matrix DAG này ta có thể vẽ ra sơ đồ mạng tùy thích bằng package igraph và ggraph.

net$dag

##          wells pe d.dimer pregnant angio treat death
## wells        0  1       0        0     0     0     0
## pe           0  0       1        0     1     0     1
## d.dimer      0  0       0        0     0     1     0
## pregnant     0  0       1        0     0     0     0
## angio        0  0       0        0     0     1     0
## treat        0  0       0        0     0     0     1
## death        0  0       0        0     0     0     0

3 Bước 2: Chuyển sơ đồ mạng thành mạng mô hình Bayes

Dựa vào sơ đồ mạng mô hình như trên, ta có thể hình dung về nó như một tập hợp 5 mô hình tương ứng với 5 nút kết quả: pe, angio, d.dimer, treat và death. Như vậy nếu ta xác định thành công 5 mô hình con này, xem như mục tiêu đã đạt được.

3.1 Phân tích hồi quy với hàm GLM

3.1.1 PE

Mô hình cho liên kết Wells - PE có dạng một mô hình logistic đơn biến, với giả định kết quả xác suất PE được ước tính bằng phân bố binomial, link function là hàm logit

binomial_smooth <- function(...) {
  geom_smooth(method = "glm", method.args = list(family = "binomial"), 
              col="black",fill="gold",linetype=2)}

PE%>%mutate(diag=as.numeric(pe)-1)%>%
  ggplot(aes(x=wells,y=diag))+
  geom_point(aes(col=pe))+
  binomial_smooth(formula = y ~ x,alpha=0.3)+
  theme_bw()

pe_lm<-glm(pe ~ wells, family="binomial",data=PE)

summary(pe_lm)

## 
## Call:
## glm(formula = pe ~ wells, family = "binomial", data = PE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9995  -0.6063  -0.4636  -0.3517   2.5935  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.90355    0.08308  -46.98   <2e-16 ***
## wells        0.57570    0.01754   32.81   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9524.2  on 9999  degrees of freedom
## Residual deviance: 8238.9  on 9998  degrees of freedom
## AIC: 8242.9
## 
## Number of Fisher Scoring iterations: 5

pe_lm%>%broom::tidy()%>%
  mutate(OR=exp(estimate),
         UL=exp(estimate+1.96*std.error),
         LL=exp(estimate-1.96*std.error))%>%
  select(term,OR,LL,UL,p.value)

3.1.2 Angio

Mô hình cho liên kết pe-angio cũng là một mô hình logistic đơn biến:

angio_lm<-glm(angio ~ pe, family="binomial",data=PE)

summary(angio_lm)

## 
## Call:
## glm(formula = angio ~ pe, family = "binomial", data = PE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6469  -0.4528  -0.4528  -0.4528   2.1580  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -2.22585    0.03730  -59.67   <2e-16 ***
## peYes        3.28411    0.06516   50.40   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 10425.6  on 9999  degrees of freedom
## Residual deviance:  7309.1  on 9998  degrees of freedom
## AIC: 7313.1
## 
## Number of Fisher Scoring iterations: 4

angio_lm%>%broom::tidy()%>%
  mutate(OR=exp(estimate),
         UL=exp(estimate+1.96*std.error),
         LL=exp(estimate-1.96*std.error))%>%
  select(term,OR,LL,UL,p.value)

3.1.3 D.Dimer

Mô hình cho liên kết (d.dimer | pregnant*pe) là một mô hình hồi quy tuyến tính, với giả định biến kết quả d.dimer được ước lượng bằng phân bố Gaussian, phụ thuộc vào 2 biến pregnant và pe:

ddimer_lm<-glm(d.dimer ~ pregnant + pe,data=PE)

summary(ddimer_lm)

## 
## Call:
## glm(formula = d.dimer ~ pregnant + pe, data = PE)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -103.243   -19.622    -0.243    19.757   107.757  
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 210.2425     0.3450  609.46   <2e-16 ***
## pregnantYes  29.2950     1.0021   29.23   <2e-16 ***
## peYes        68.3794     0.7723   88.54   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 892.5316)
## 
##     Null deviance: 16687526  on 9999  degrees of freedom
## Residual deviance:  8922638  on 9997  degrees of freedom
## AIC: 96324
## 
## Number of Fisher Scoring iterations: 2

3.1.4 Treat

Mô hình cho liên kết (treat | d.dimer*angio) có dạng logistic, 2 biến :

treat_lm<-glm(treat ~ angio + d.dimer,family="binomial",data=PE)

summary(treat_lm)

## 
## Call:
## glm(formula = treat ~ angio + d.dimer, family = "binomial", data = PE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3985  -0.6798  -0.5073   0.6114   2.6499  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -5.8931574  0.1652900  -35.65   <2e-16 ***
## angioPositive  1.7335361  0.0596181   29.08   <2e-16 ***
## d.dimer        0.0199375  0.0007072   28.19   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 12247.6  on 9999  degrees of freedom
## Residual deviance:  9466.8  on 9997  degrees of freedom
## AIC: 9472.8
## 
## Number of Fisher Scoring iterations: 4

treat_lm%>%broom::tidy()%>%
  mutate(OR=exp(estimate),
         UL=exp(estimate+1.96*std.error),
         LL=exp(estimate-1.96*std.error))%>%
  select(term,OR,LL,UL,p.value)

3.1.5 Death

death_lm<-glm(death ~ pe + treat,family="binomial",data=PE)
             
summary(death_lm)

## 
## Call:
## glm(formula = death ~ pe + treat, family = "binomial", data = PE)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7525  -0.1736  -0.1736  -0.0662   3.5002  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -4.1876     0.1008  -41.53   <2e-16 ***
## peYes         5.4808     0.1406   38.99   <2e-16 ***
## treatYes     -1.9358     0.1163  -16.64   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 6382.2  on 9999  degrees of freedom
## Residual deviance: 3278.3  on 9997  degrees of freedom
## AIC: 3284.3
## 
## Number of Fisher Scoring iterations: 7

death_lm%>%broom::tidy()%>%
  mutate(OR=exp(estimate),
         UL=exp(estimate+1.96*std.error),
         LL=exp(estimate-1.96*std.error))%>%
  select(term,OR,LL,UL,p.value)

3.2 Mạng mô hình Bayes sử dụng brms và STAN

Bây giờ, Nhi sẽ dựng một mạng mô hình Bayes thực sự dùng ngôn ngữ STAN qua giao thức là package brms.

Trong các bài trước, Nhi đã từng giới thiệu cho các bạn package brms là một công cụ rất mạnh mẽ, linh hoạt cho phân tích Bayes. Một trong những tính năng mà brms hỗ trợ là dựng mô hình mạng lưới cấu trúc (bạn có thể gọi phương pháp này bằng nhiều cái tên: SEM, DAG, Latent variable analysis hay Bayesian network…, nhưng cách thực hiện cũng như nhau). Tính năng này có bản chất là ước lượng (MCMC sampling) đồng thời cho hàng loạt mô hình nằm trong một mạng lưới cấu trúc phức tạp. brms cho phép áp dụng priors, family và công thức riêng biệt cho mỗi mô hình bên trong mạng lưới, do đó hết sức linh hoạt khi đối diện với các bài toán phức tạp.

Trong bài thực hành này, Nhi đơn giản vấn đề bằng cách chỉ áp dụng 2 kiểu mô hình: logistic và Gaussian, và không tùy chỉnh sâu về priors.

Trước hết, ta viết công thức cho 5 mô hình theo cú pháp brms:

library(brms)

pe_mod<-bf(pe ~ wells, family=bernoulli(link = "logit"))
angio_mod<-bf(angio ~ pe, family=bernoulli(link = "logit"))
ddimer_mod<-bf(d.dimer ~ pregnant + pe)
treat_mod<-bf(treat ~ angio + d.dimer,family=bernoulli(link = "logit"))
death_mod<-bf(death ~ pe + treat,family=bernoulli(link = "logit"))

Sau đó, ta dùng hàm brm để đóng gói tất cả mô hình trong cùng 1 chế độ lấy mẫu MCMC: với 1 chuỗi, 1000 lượt khởi động, 2000 lượt lấy mẫu, sau đó rút gọn 50%.

brms sẽ compile mô hình thành STAN code, sau đó thi hành. Tiến trình tương đối chậm, tùy vào máy tính mạnh hay yếu (Nhi mất khoảng 15 phút), nhưng mô hình rồi cũng sẽ converge xong.

k_fit_brms <- brm(pe_mod +
                    angio_mod +
                    ddimer_mod +
                    treat_mod +
                    death_mod+
                    set_rescor(FALSE),
                  data=PE,
                  iter = 3000, 
                  warmup = 1000,
                  chains = 1, 
                  thin = 2,
                  refresh = 0,
                  control=list(adapt_delta=0.95,
                               max_treedepth=20))

## 
## Gradient evaluation took 0.01 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 100 seconds.
## Adjust your expectations accordingly!
## 
## 
## 
##  Elapsed Time: 1021.12 seconds (Warm-up)
##                528.894 seconds (Sampling)
##                1550.01 seconds (Total)

4 Bước 3: Khảo sát nội dung mạng mô hình Bayes

Trước hết, Nhi trích nội dung STAN code của mô hình vừa dựng, code này chỉ nhằm thỏa mãn sự tò mò của các bạn thích học viết code thủ công bằng STAN

stancode(k_fit_brms )

## // generated with brms 2.3.1
## functions { 
## } 
## data { 
##   int<lower=1> N;  // total number of observations 
##   int Y_pe[N];  // response variable 
##   int<lower=1> K_pe;  // number of population-level effects 
##   matrix[N, K_pe] X_pe;  // population-level design matrix 
##   int Y_angio[N];  // response variable 
##   int<lower=1> K_angio;  // number of population-level effects 
##   matrix[N, K_angio] X_angio;  // population-level design matrix 
##   vector[N] Y_ddimer;  // response variable 
##   int<lower=1> K_ddimer;  // number of population-level effects 
##   matrix[N, K_ddimer] X_ddimer;  // population-level design matrix 
##   int Y_treat[N];  // response variable 
##   int<lower=1> K_treat;  // number of population-level effects 
##   matrix[N, K_treat] X_treat;  // population-level design matrix 
##   int Y_death[N];  // response variable 
##   int<lower=1> K_death;  // number of population-level effects 
##   matrix[N, K_death] X_death;  // population-level design matrix 
##   int prior_only;  // should the likelihood be ignored? 
## } 
## transformed data { 
##   int Kc_pe = K_pe - 1; 
##   matrix[N, K_pe - 1] Xc_pe;  // centered version of X_pe 
##   vector[K_pe - 1] means_X_pe;  // column means of X_pe before centering 
##   int Kc_angio = K_angio - 1; 
##   matrix[N, K_angio - 1] Xc_angio;  // centered version of X_angio 
##   vector[K_angio - 1] means_X_angio;  // column means of X_angio before centering 
##   int Kc_ddimer = K_ddimer - 1; 
##   matrix[N, K_ddimer - 1] Xc_ddimer;  // centered version of X_ddimer 
##   vector[K_ddimer - 1] means_X_ddimer;  // column means of X_ddimer before centering 
##   int Kc_treat = K_treat - 1; 
##   matrix[N, K_treat - 1] Xc_treat;  // centered version of X_treat 
##   vector[K_treat - 1] means_X_treat;  // column means of X_treat before centering 
##   int Kc_death = K_death - 1; 
##   matrix[N, K_death - 1] Xc_death;  // centered version of X_death 
##   vector[K_death - 1] means_X_death;  // column means of X_death before centering 
##   for (i in 2:K_pe) { 
##     means_X_pe[i - 1] = mean(X_pe[, i]); 
##     Xc_pe[, i - 1] = X_pe[, i] - means_X_pe[i - 1]; 
##   } 
##   for (i in 2:K_angio) { 
##     means_X_angio[i - 1] = mean(X_angio[, i]); 
##     Xc_angio[, i - 1] = X_angio[, i] - means_X_angio[i - 1]; 
##   } 
##   for (i in 2:K_ddimer) { 
##     means_X_ddimer[i - 1] = mean(X_ddimer[, i]); 
##     Xc_ddimer[, i - 1] = X_ddimer[, i] - means_X_ddimer[i - 1]; 
##   } 
##   for (i in 2:K_treat) { 
##     means_X_treat[i - 1] = mean(X_treat[, i]); 
##     Xc_treat[, i - 1] = X_treat[, i] - means_X_treat[i - 1]; 
##   } 
##   for (i in 2:K_death) { 
##     means_X_death[i - 1] = mean(X_death[, i]); 
##     Xc_death[, i - 1] = X_death[, i] - means_X_death[i - 1]; 
##   } 
## } 
## parameters { 
##   vector[Kc_pe] b_pe;  // population-level effects 
##   real temp_pe_Intercept;  // temporary intercept 
##   vector[Kc_angio] b_angio;  // population-level effects 
##   real temp_angio_Intercept;  // temporary intercept 
##   vector[Kc_ddimer] b_ddimer;  // population-level effects 
##   real temp_ddimer_Intercept;  // temporary intercept 
##   real<lower=0> sigma_ddimer;  // residual SD 
##   vector[Kc_treat] b_treat;  // population-level effects 
##   real temp_treat_Intercept;  // temporary intercept 
##   vector[Kc_death] b_death;  // population-level effects 
##   real temp_death_Intercept;  // temporary intercept 
## } 
## transformed parameters { 
## } 
## model { 
##   vector[N] mu_pe = Xc_pe * b_pe + temp_pe_Intercept; 
##   vector[N] mu_angio = Xc_angio * b_angio + temp_angio_Intercept; 
##   vector[N] mu_ddimer = Xc_ddimer * b_ddimer + temp_ddimer_Intercept; 
##   vector[N] mu_treat = Xc_treat * b_treat + temp_treat_Intercept; 
##   vector[N] mu_death = Xc_death * b_death + temp_death_Intercept; 
##   // priors including all constants 
##   target += student_t_lpdf(temp_pe_Intercept | 3, 0, 10); 
##   target += student_t_lpdf(temp_angio_Intercept | 3, 0, 10); 
##   target += student_t_lpdf(temp_ddimer_Intercept | 3, 221, 39); 
##   target += student_t_lpdf(sigma_ddimer | 3, 0, 39)
##     - 1 * student_t_lccdf(0 | 3, 0, 39); 
##   target += student_t_lpdf(temp_treat_Intercept | 3, 0, 10); 
##   target += student_t_lpdf(temp_death_Intercept | 3, 0, 10); 
##   // likelihood including all constants 
##   if (!prior_only) { 
##     target += bernoulli_logit_lpmf(Y_pe | mu_pe); 
##     target += bernoulli_logit_lpmf(Y_angio | mu_angio); 
##     target += normal_lpdf(Y_ddimer | mu_ddimer, sigma_ddimer); 
##     target += bernoulli_logit_lpmf(Y_treat | mu_treat); 
##     target += bernoulli_logit_lpmf(Y_death | mu_death); 
##   } 
## } 
## generated quantities { 
##   // actual population-level intercept 
##   real b_pe_Intercept = temp_pe_Intercept - dot_product(means_X_pe, b_pe); 
##   // actual population-level intercept 
##   real b_angio_Intercept = temp_angio_Intercept - dot_product(means_X_angio, b_angio); 
##   // actual population-level intercept 
##   real b_ddimer_Intercept = temp_ddimer_Intercept - dot_product(means_X_ddimer, b_ddimer); 
##   // actual population-level intercept 
##   real b_treat_Intercept = temp_treat_Intercept - dot_product(means_X_treat, b_treat); 
##   // actual population-level intercept 
##   real b_death_Intercept = temp_death_Intercept - dot_product(means_X_death, b_death); 
## }

Kết quả thô của mạng mô hình đa biến được tóm tắt bằng hàm summary như sau: Dù sử dụng quy luật phân bố hơi khác (bernoulli thay vì binomial), các tham số trong mô hình có thể xem như tương đương với kết quả của hàm glm cho từng mô hình đơn lẻ.

tên các effects trong kết quả được trình bày theo quy ước: (tên mô hình_hiệu ứng), thí dụ death_peYes có nghĩa là tham số cho hiệu ứng Pe=Yes trong mô hình tiên lượng nguy cơ tử vong.

summary(k_fit_brms)

##  Family: MV(bernoulli, bernoulli, gaussian, bernoulli, bernoulli) 
##   Links: mu = logit
##          mu = logit
##          mu = identity; sigma = identity
##          mu = logit
##          mu = logit 
## Formula: pe ~ wells 
##          angio ~ pe 
##          d.dimer ~ pregnant + pe 
##          treat ~ angio + d.dimer 
##          death ~ pe + treat 
##    Data: PE (Number of observations: 10000) 
## Samples: 1 chains, each with iter = 3000; warmup = 1000; thin = 2;
##          total post-warmup samples = 1000
## 
## Population-Level Effects: 
##                     Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## pe_Intercept           -3.90      0.08    -4.06    -3.74       1000 1.00
## angio_Intercept        -2.22      0.04    -2.30    -2.16        902 1.00
## ddimer_Intercept      210.25      0.34   209.59   210.87       1000 1.00
## treat_Intercept        -5.90      0.16    -6.21    -5.59       1000 1.00
## death_Intercept        -4.19      0.10    -4.39    -4.00        910 1.00
## pe_wells                0.58      0.02     0.54     0.61       1000 1.00
## angio_peYes             3.28      0.07     3.15     3.42       1000 1.00
## ddimer_pregnantYes     29.28      0.99    27.36    31.16        795 1.01
## ddimer_peYes           68.37      0.73    66.98    69.77        986 1.00
## treat_angioPositive     1.73      0.06     1.61     1.85        990 1.00
## treat_d.dimer           0.02      0.00     0.02     0.02       1000 1.00
## death_peYes             5.48      0.14     5.22     5.77        947 1.00
## death_treatYes         -1.94      0.12    -2.17    -1.71        932 1.00
## 
## Family Specific Parameters: 
##              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma_ddimer    29.88      0.20    29.47    30.26        859 1.00
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Nhi trích xuất các chuỗi MCMC cho một số hiệu ứng cần khảo sát, và tính Odds ratio cho các hiệu ứng thuộc mô hình logistic bằng hàm exp(), sau đó vẽ biểu đồ phân phối hậu nghiệm và diễn tiến chuỗi MCMC như sau:

postdf=k_fit_brms$fit%>%as.data.frame()%>%
  as_tibble()%>%
  mutate(Iteration=as.numeric(c(1:1000)))

postdf[,c(1,2,4:7,10:13)]<-exp(postdf[,c(1,2,4:7,10:13)])

longpost=postdf%>%
  gather(b_pe_wells:b_death_treatYes,
         key="Variable",value="effect")

library(ggridges)

ggplot(longpost)+
  geom_density(aes(x=effect,fill=Variable),col="black",show.legend = F,alpha=0.5)+
  facet_wrap(~Variable,ncol=2,scales = "free")+
  scale_fill_manual(values=pals::tol(8))+theme_bw()

ggplot(longpost)+
  geom_path(aes(y=effect,x=Iteration,col=Variable),
            show.legend = F,alpha=0.8)+
  facet_wrap(~Variable,ncol=2,scales = "free")+
  scale_color_manual(values=pals::tol(8))+theme_bw()

Phân phối hậu nghiệm cho các điều kiện/hiệu ứng này hoàn toàn hợp lý, và có thể diễn giải như thường lệ, thí dụ: bệnh nhân bị thuyên tắc phổi có ddimer tăng trung bình 68 đơn vị, nguy cơ tử vong cao đến 240 lần. Việc điều trị làm giảm nguy cơ tử vong với Odds-ratio = 0.14, …

longpost%>%group_by(Variable)%>%
  summarise_at("effect",funs(Median=median,
                           SD=sd,
                           LL=quantile(.,probs=0.025),
                           UL=quantile(.,probs=0.975)))

Trước khi kết thúc, Nhi làm thêm một thủ thuật nhỏ nữa , để trình bày nội dung mô hình mạng dưới dạng một sơ đồ DAG có kèm Odds-ratio và hiệu ứng của mỗi mô hình:

5 Kết luận

Bài thực hành đến đây là hết, chúng ta vừa thử nghiệm tính năng mô hình cấu trúc đa biến của package brms. Tính năng này khá lợi hại, vì nó cho phép chúng ta lai ghép giữa thống kê Bayes và những phương pháp mô hình cấu trúc khác như SEM, DAG, Latent var. analysis, network … với tính linh hoạt rất lớn, thậm chí hơn cả những package sử dụng JAGs sampler.

Về giá trị thực tiễn, phương pháp mạng mô hình Bayes (Bayesian network)mở ra nhiều ứng dụng tiềm năng cho nghiên cứu Y học lâm sàng, vì cho phép các bạn kể một câu chuyện hấp dẫn và hợp lý với dữ liệu hiện có, suy diễn nhân quả, biện luận quy trình đưa ra quyết định chẩn đoán, phác đồ điều trị…bằng kết quả phân tích đa biến,dựa vào lý thuyết y học và tất cả những ưu thế mà trường phái Bayes có thể mang lại.

Chúc các bạn thành công.