Mảnh ghép cuối cùng cho suy diễn Bayes
Lê Ngọc Khả Nhi
13 Tháng 6 2018
1 Giới thiệu
Khi thực hiện project Bayes for Vietnam, Nhi có hy vọng khuyến khích anh chị em đồng nghiệp chuyển sang dùng Thống kê theo trường phái Bayes thay cho những kiểm định rời rạc khác. Qua những bài thực hành trước đây, chúng ta đã thấy khả năng mạnh mẽ và linh hoạt của phương pháp Bayes. Tuy nhiên mục tiêu phổ cập phương pháp này cho số đông sẽ khó thành công nếu buộc người học phải viết code thủ công hoàn toàn cho nội dung mô hình và quy trình khai thác kết quả. Vì lý do này, Nhi sẽ dần chuyển từ cách làm thủ công với chính ngôn ngữ STAN sang sử dụng những công cụ đơn giản hơn, thí dụ package brms.
Như các bạn đã biết, suy diễn thống kê theo Bayes có hai bước, đầu tiên ta phải có một mô hình (bạn có thể dựng mô hình bằng rstan, stanarm hay brms), sau đó từ mô hình này ta có thể khảo sát phân bố hậu nghiệm từ các chuỗi MCMC để kiểm định một giả thuyết tùy chọn.
Một số phương pháp đã được giới thiệu cho bước khai thác mô hình, bao gồm mô tả trực quan phân bố hậu nghiệm, tóm tắt khoảng mật độ cao nhất (Highest density Interval, HDI), kiểm tra giả thuyết vô hiệu với ngưỡng đơn lẻ (Compval) hay khoảng vô hiệu thực dụng (ROPE), giá trị p_value hay Bayes factor. Tuy nhiên, tất cả những quy trình này cần phải viết hàm thủ công.
Trong phiên bản mới nhất của package sjstats, tác giả Daniel Lüdecke đã bổ sung một loạt hàm cho phép suy diễn Bayes trực tiếp từ object mô hình brms và rstanarm. Bài thực hành hôm nay sẽ giới thiệu với các bạn về các hàm này.
library(tidyverse)
library(rstan)
library(brms)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())Trước hết, ta sẽ dựng 2 mô hình hồi quy Bayes bằng brms:
Mô hình thứ nhất là một thiết kế MANCOVA có nội dung như sau: Đây là một thí nghiệm sinh lý trên động vật vào năm 1994 của GS. J. Ludbrook. mục tiêu nhằm khảo sát hiệu ứng tương tác của liều hoạt chất phenylbiguanide (kích thích co mạch) và 1 chất đối vận serotonin (có tác dụng giảm co mạch) gây thay đổi huyết áp ở 5 con thỏ.
Biến kết quả là BPchange=thay đổi huyết áp, yếu tố thứ nhất là Treatment chỉ phân nhóm sử dụng MDL đối kháng thụ thể serotonin và Dose chỉ liều phenylbiguanid. Do thí nghiệm lặp lại 5 lần ở mỗi con vật nên ta có thểm random effect trong mô hình.
mdl_df<-read.csv("http://vincentarelbundock.github.io/Rdatasets/csv/MASS/Rabbit.csv")
head(mdl_df)%>%knitr::kable()| X | BPchange | Dose | Run | Treatment | Animal |
|---|---|---|---|---|---|
| 1 | 0.5 | 6.25 | C1 | Control | R1 |
| 2 | 4.5 | 12.50 | C1 | Control | R1 |
| 3 | 10.0 | 25.00 | C1 | Control | R1 |
| 4 | 26.0 | 50.00 | C1 | Control | R1 |
| 5 | 37.0 | 100.00 | C1 | Control | R1 |
| 6 | 32.0 | 200.00 | C1 | Control | R1 |
mod<-brm(BPchange~Treatment*Dose+(1|Animal),
data=mdl_df,
family = brmsfamily("gaussian", link_sigma = "log"),
seed = 1234,iter=1500,chains = 1
)##
## SAMPLING FOR MODEL 'gaussian brms-model' NOW (CHAIN 1).
##
## Gradient evaluation took 0.02 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 200 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 1500 [ 0%] (Warmup)
## Iteration: 150 / 1500 [ 10%] (Warmup)
## Iteration: 300 / 1500 [ 20%] (Warmup)
## Iteration: 450 / 1500 [ 30%] (Warmup)
## Iteration: 600 / 1500 [ 40%] (Warmup)
## Iteration: 750 / 1500 [ 50%] (Warmup)
## Iteration: 751 / 1500 [ 50%] (Sampling)
## Iteration: 900 / 1500 [ 60%] (Sampling)
## Iteration: 1050 / 1500 [ 70%] (Sampling)
## Iteration: 1200 / 1500 [ 80%] (Sampling)
## Iteration: 1350 / 1500 [ 90%] (Sampling)
## Iteration: 1500 / 1500 [100%] (Sampling)
##
## Elapsed Time: 0.384 seconds (Warm-up)
## 0.12 seconds (Sampling)
## 0.504 seconds (Total)summary(mod)## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: BPchange ~ Treatment * Dose + (1 | Animal)
## Data: mdl_df (Number of observations: 60)
## Samples: 1 chains, each with iter = 1500; warmup = 750; thin = 1;
## total post-warmup samples = 750
##
## Group-Level Effects:
## ~Animal (Number of levels: 5)
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sd(Intercept) 2.83 1.99 0.34 8.84 86 1.00
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 4.26 2.41 0.28 11.44 76 1.03
## TreatmentMDL -4.22 2.18 -8.36 0.47 330 1.00
## Dose 0.14 0.02 0.11 0.18 274 1.00
## TreatmentMDL:Dose -0.01 0.02 -0.05 0.04 376 1.00
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 5.89 0.57 4.92 7.10 473 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).Thứ hai là một mô hình logistic với phân bố bernoulli (link function là hàm logit), để khảo sát liên hệ giữa các yếu tố: Giới tính, bệnh Ung thư, Suy thận, Nhiễm trùng, Gãy xương, kỹ thuật CPR đối với nguy cơ tử vong của bệnh nhân cấp cứu tại khoa ICU.
icu_df<-read.csv("https://www.openml.org/data/get_csv/53980/ICU.csv")%>%select(STA,SEX,CAN,CRN,INF,FRA)
icu_df$CAN<-icu_df$CAN%>%recode_factor(.,"1"="NO","2"="YES")
icu_df$CRN<-icu_df$CRN%>%recode_factor(.,"1"="NO","2"="YES")
icu_df$INF<-icu_df$INF%>%recode_factor(.,"1"="NO","2"="YES")
icu_df$FRA<-icu_df$FRA%>%recode_factor(.,"1"="NO","2"="YES")
icu_df$SEX<-icu_df$SEX%>%recode_factor(.,"1"="M","2"="F")
head(icu_df)%>%knitr::kable()| STA | SEX | CAN | CRN | INF | FRA |
|---|---|---|---|---|---|
| 0 | F | NO | NO | YES | NO |
| 0 | M | NO | NO | NO | NO |
| 0 | M | NO | NO | NO | NO |
| 0 | M | NO | NO | YES | YES |
| 0 | F | NO | NO | YES | NO |
| 0 | M | NO | NO | YES | NO |
logmod=brm(STA~CAN+CRN+INF+FRA+SEX,
family=bernoulli(link = "logit"),data=icu_df,
seed = 1234,iter=1500,chains = 1
)##
## SAMPLING FOR MODEL 'bernoulli brms-model' NOW (CHAIN 1).
##
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 1500 [ 0%] (Warmup)
## Iteration: 150 / 1500 [ 10%] (Warmup)
## Iteration: 300 / 1500 [ 20%] (Warmup)
## Iteration: 450 / 1500 [ 30%] (Warmup)
## Iteration: 600 / 1500 [ 40%] (Warmup)
## Iteration: 750 / 1500 [ 50%] (Warmup)
## Iteration: 751 / 1500 [ 50%] (Sampling)
## Iteration: 900 / 1500 [ 60%] (Sampling)
## Iteration: 1050 / 1500 [ 70%] (Sampling)
## Iteration: 1200 / 1500 [ 80%] (Sampling)
## Iteration: 1350 / 1500 [ 90%] (Sampling)
## Iteration: 1500 / 1500 [100%] (Sampling)
##
## Elapsed Time: 0.262 seconds (Warm-up)
## 0.161 seconds (Sampling)
## 0.423 seconds (Total)summary(logmod)## Family: bernoulli
## Links: mu = logit
## Formula: STA ~ CAN + CRN + INF + FRA + SEX
## Data: icu_df (Number of observations: 200)
## Samples: 1 chains, each with iter = 1500; warmup = 750; thin = 1;
## total post-warmup samples = 750
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept -2.02 0.34 -2.71 -1.38 750 1.00
## CANYES 0.07 0.61 -1.22 1.24 750 1.00
## CRNYES 1.13 0.52 0.13 2.09 750 1.00
## INFYES 0.89 0.39 0.17 1.64 750 1.00
## FRAYES -0.05 0.74 -1.50 1.25 750 1.00
## SEXF 0.03 0.38 -0.72 0.72 750 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).Tiếp theo, các bạn install 2 package sjstats và sjPlot vào R, đây là 2 công cụ rất tiện dụng cho thực hành thống kê. Chúng cung cấp các hàm cho phép đơn giản hóa nhiều quy trình khai thác kết quả của mô hình trong R.
library(sjstats)
library(sjPlot)2 Hàm hdi(): khảo sát phân phối hậu nghiệm
Hàm hdi() có công dụng tóm tắt thông tin về phân phối hậu nghiệm của các tham số trong mô hình (hay hiệu ứng ta cần khảo sát). HDI là một thuật ngữ của John Kruschke, với ý nghĩa gần giống như Confidence Interval trong thống kê cổ điển. Kruschke khuyến cáo dùng thuật ngữ Highest density interval để mô tả đúng bản chất mật độ một phân phối liên tục thực sự của hiệu ứng từ chuỗi MCMC, thay cho ý nghĩa “khoảng tin cậy” là kết quả ước lượng từ sự lấy mẫu ngẫu nhiên lặp lại trong phái Frequentist.
Ta có thể tùy chọn ngưỡng cho HDI: thí dụ 90% hay 95%
hdi(mod,prob=c(0.90, 0.95, 0.975),type="fixed")%>%print(digits=3)##
## # Highest Density Interval
##
## HDI(90%) HDI(95%) HDI(98%)
## b_Intercept [ 0.441 7.257] [-0.687 7.907] [ 0.087 12.591]
## b_TreatmentMDL [-7.997 -0.740] [-8.690 -0.080] [-8.566 1.245]
## b_Dose [ 0.116 0.169] [ 0.114 0.175] [ 0.106 0.178]
## b_TreatmentMDL.Dose [-0.048 0.029] [-0.055 0.038] [-0.057 0.051]
## sigma [ 4.952 6.775] [ 4.759 6.875] [ 4.711 7.161]Cho mô hình hồi quy logistic, ta có thể áp dụng một hàm hoán chuyển exp( )để có HDI cho Odds-ratio:
hdi(logmod,prob=c(0.8, 0.9, 0.95),type="fixed",trans="exp")%>%print(digits=3)##
## # Highest Density Interval
##
## HDI(80%) HDI(90%) HDI(95%)
## b_Intercept [0.086 0.203] [0.074 0.223] [0.067 0.253]
## b_CANYES [0.543 2.340] [0.394 2.850] [0.311 3.648]
## b_CRNYES [1.556 6.014] [1.461 8.099] [1.165 8.301]
## b_INFYES [1.529 3.977] [1.362 4.792] [1.207 5.329]
## b_FRAYES [0.406 2.560] [0.327 3.410] [0.223 3.505]
## b_SEXF [0.665 1.712] [0.586 1.978] [0.484 2.054]Nhận xét: hàm hdi() rất tiện dụng cho mục tiêu khảo sát nhanh ý nghĩa thống kê của phân phối hậu nghiệm. So với những hàm tóm tắt chuỗi MCMC của brms, broom hay coda, hàm hdi linh hoạt hơn khi cho phép người dùng ấn định nhiều ngưỡng probability; tuy nhiên nó không cung cấp giá trị trung tâm (trung bình, trung vị) của phân bố này.
3 Hàm rope(): khoảng vô hiệu thực dụng
Hàm rope() cho phép áp dụng một khoảng tùy chọn [-b,b] với -b và b là hai ngưỡng trên và dưới, với giả định là phần phân phối hậu nghiệm rơi vào bên trong khoảng này sẽ được xem như vô hiệu (hay không có ý nghĩa) về mặt thực dụng. Thí dụ nếu Odds-ratio nằm trong khoảng [-1.5,1.5 ]cho biết ta chỉ chấp nhận OR có ý nghĩa nếu giá trị của nó > 1.5 hay <-1.5.
Thí dụ cho mô hình MANCOVA, ta đặt rope = -2 đến 1, kết quả như sau
rope(mod,rope=c(-2,1),type="fixed")%>%print(digits=3)##
## # Proportions of samples inside and outside the ROPE
##
## inside outside
## b_Intercept 6.667% 93.333%
## b_TreatmentMDL 13.467% 86.533%
## b_Dose 100.000% 0.000%
## b_TreatmentMDL.Dose 100.000% 0.000%
## sigma 0.000% 100.000%Kết quả này được diễn giải là: Hoạt chất MDL làm giảm huyết áp với hơn 86 % phân phối hậu nghiệm nằm ngoài ngưỡng -2
Nhận xét: So với code thủ công, hàm rope của package sjstats không xác định được tỉ lệ mật độ phân bố nằm trên và dưới khoảng rope, nhưng chỉ mới tính được bao nhiêu % nằm trong hay ngoài ROPE.
Tương tự, ta áp dụng ROPE=-2 đến 2 cho OR của mô hình logistic như sau:
rope(logmod,rope=c(-2,2),trans="exp")%>%print(digits=3)##
## # Proportions of samples inside and outside the ROPE
##
## inside outside
## b_Intercept 100.000% 0.000%
## b_CANYES 84.267% 15.733%
## b_CRNYES 20.267% 79.733%
## b_INFYES 31.600% 68.400%
## b_FRAYES 84.800% 15.200%
## b_SEXF 97.067% 2.933%Theo kết quả này, chỉ có yếu tố Bệnh lý Thận mạn tính là có OR đáng kể nhất (80% nằm ngoài ngưỡng -2:2).
4 Hàm Equivalence test: kiểm định giả thuyết vô hiệu
Hàm Equivalence test thực hiện một kiểm định để suy diễn thống kê tự động cho mỗi hiệu ứng trong mô hình Bayes, dựa vào 1 khoảng vô hiệu thực dụng chung. Giả thuyết cần kiểm định là liệu phân bố hậu nghiệm có nằm ngoài ROPE hay không ? Kết quả có thể là: Chấp nhận giả thuyết, bác bỏ giả thuyết hoặc Hoài nghi, không thể quyết định.
equi_test(mod,out="plot",rope=c(-2,-1),eff_size = 0.1)
equi_test(mod,out="txt",rope=c(-2,-1),eff_size = 0.1)##
## # Test for Practical Equivalence of Model Predictors
##
## Effect Size: 0.10
## ROPE: [-2.00 -1.00]
## Samples: 750
##
## H0 %inROPE HDI(95%)
## b_Intercept (*) reject 0.00 [-0.69 7.91]
## b_TreatmentMDL (*) undecided 7.73 [-8.69 -0.08]
## b_Dose (*) reject 0.00 [ 0.11 0.18]
## b_TreatmentMDL.Dose (*) reject 0.00 [-0.06 0.04]
## sigma reject 0.00 [ 4.76 6.87]Tương tự, ta có thể suy diễn thống kê cho OR của mô hình logistic
equi_test(logmod,out="plot",rope=c(-1,1),eff_size = 0.1,trans="exp")
equi_test(logmod,out="txt",rope=c(-1,1),eff_size = 0.1,trans="exp")##
## # Test for Practical Equivalence of Model Predictors
##
## Effect Size: 0.10
## ROPE: [-1.00 1.00]
## Samples: 750
##
## H0 %inROPE HDI(95%)
## b_Intercept reject 0.00 [-2.71 -1.38]
## b_CANYES undecided 89.87 [-1.17 1.29]
## b_CRNYES undecided 40.67 [ 0.15 2.12]
## b_INFYES undecided 60.27 [ 0.19 1.67]
## b_FRAYES undecided 84.00 [-1.50 1.25]
## b_SEXF accept 99.07 [-0.73 0.72]Nhận xét: Đây là một ý tưởng thú vị, vì kết quả của kiểm định có thể giúp đưa ra quyết định về “ý nghĩa thống kê” của hiệu ứng trong mô hình, dựa vào ROPE. Tuy nhiên, nhược điểm của cách làm này đó là chỉ có 1 ngưỡng vô hiệu chung cho mọi hiệu ứng. Phương pháp suy diễn bằng Bayes factor linh hoạt hơn nhiều. Tuy biểu đồ phân bố hậu nghiệm được cung cấp, nó có phẩm chất không tốt bằng các biểu đồ được vẽ thủ công từ ggplot2 và dữ liệu MCMC thô.
5 Hàm tidy_stan: Tóm tắt mô hình
Hàm tidy_stan cho phép tóm tắt dễ dàng nội dung mô hình , bản thân nó bao gồm cả công dụng của hàm hdi nhưng bổ sung thêm giá trị Trung vị của phân bố hậu nghiệm, và 3 trị số kiểm tra phẩm chất chuỗi MCMC:
tidy_stan(mod,prob=0.95)##
## # Summary Statistics of Stan-Model
##
## estimate std.error HDI(95%) neff_ratio Rhat mcse
## Intercept 4.23 2.02 [-0.69 7.91] 0.10 1.03 0.28
## TreatmentMDL -4.22 2.06 [-8.69 -0.08] 0.44 1.00 0.12
## Dose 0.14 0.02 [ 0.11 0.18] 0.36 1.00 0.00
## TreatmentMDL.Dose -0.01 0.02 [-0.06 0.04] 0.50 1.00 0.00
## sigma 5.87 0.55 [ 4.76 6.87] 0.63 1.00 0.03Ngoài hiệu ứng chính ta còn có thể tóm tắt random effects
tidy_stan(mod,prob=0.95,type="random")##
## # Summary Statistics of Stan-Model
##
## ## Random effect (Intercept)
##
## estimate std.error HDI(95%) neff_ratio Rhat mcse
## Animal.R1 2.31 2.09 [-0.93 8.10] 0.12 1.02 0.26
## Animal.R2 -0.41 1.58 [-6.15 3.64] 0.11 1.01 0.25
## Animal.R3 -0.12 1.50 [-4.83 3.89] 0.09 1.03 0.26
## Animal.R4 0.11 1.61 [-4.43 4.76] 0.09 1.03 0.28
## Animal.R5 -1.66 1.85 [-7.40 1.69] 0.09 1.02 0.28Với mô hình logistic, ta có thể áp dụng hàm exp() để tính OR trực tiếp và tóm tắt kết quả này:
tidy_stan(logmod,prob=0.95,trans="exp")##
## # Summary Statistics of Stan-Model
##
## estimate std.error HDI(95%) neff_ratio Rhat mcse
## Intercept 0.13 0.35 [0.07 0.25] 1 1 0.01
## CANYES 1.11 0.59 [0.31 3.65] 1 1 0.02
## CRNYES 3.11 0.55 [1.17 8.30] 1 1 0.02
## INFYES 2.42 0.40 [1.21 5.33] 1 1 0.01
## FRAYES 0.98 0.71 [0.22 3.50] 1 1 0.03
## SEXF 1.03 0.38 [0.48 2.05] 1 1 0.01Nhận xét: Hàm tidy_stan có nhiều ưu điểm, như dễ sử dụng, chỉ giữ lại những kết quả chính yếu nhất của nội dung mô hình, khả năng hoán chuyển dữ liệu trực tiếp cho phép tính Incidence rate ratio (IRR) cho count data hay Odds-ratio cho logistic model.
Hàm tidyMCMC của package broom cũng có khả năng tóm tắt nội dung mô hình Bayes nhưng cú pháp của nó phức tạp hơn nhiều
broom::tidyMCMC(mod,conf.int = TRUE,
estimate.method = "median",
conf.method = "HPDinterval",
pars=c("b_Intercept",
"b_TreatmentMDL",
"b_Dose",
"b_TreatmentMDL:Dose",
"sigma"))## term estimate std.error conf.low conf.high
## 1 b_Intercept 4.225200073 2.41111699 -0.48124134 8.03849857
## 2 b_TreatmentMDL -4.217310163 2.17556023 -8.56608659 -0.07976323
## 3 b_Dose 0.143562406 0.01654160 0.11445559 0.17534723
## 4 b_TreatmentMDL:Dose -0.007982953 0.02375522 -0.05529648 0.03753917
## 5 sigma 5.869718744 0.57264430 4.75905840 6.871253796 Hàm plot_model: khảo sát trực quan
Ngoài khả năng tóm tắt nội dung mô hình, package sjPlot còn có một hàm rất tiện dụng khác là plot_model, hàm này cho phép mô tả trực quan tất cả các loại mô hình thông dụng trong R, bao gồm mô hình STAN và brms. Nó sẽ vẽ một biểu đồ boxplot mô tả phân bố hậu nghiệm cho mỗi hiệu ứng trong mô hình:
plot_model(mod,
vline.color = "red",
sort.est = TRUE,
show.values = TRUE,
value.offset = .3)+theme_bw()
plot_model(mod,type="re",
vline.color = "red",
sort.est = TRUE,
show.values = TRUE,
value.offset = .3)+theme_bw()
Khi áp dụng cho mô hình logistic, ta chỉ cần thêm hàm hoán chuyển exp trong tùy chỉnh trans để vẽ biểu đồ cho OR:
plot_model(logmod,trans="exp",
vline.color = "red",
sort.est = TRUE,
show.values = TRUE,
value.offset = .3,
show.p=TRUE)+theme_bw()
Nhận xét: Biểu đồ trực quan cho mô hình Bayes sử dụng Boxplot thay vì point range plot vì ta đang khảo sát 1 phân bố liên tục thực sự của mỗi tham số trong mô hình. Các biểu đồ này giản dị và vừa đủ cho việc công bố kết quả. Tuy nhiên, phẩm chất mỹ thuật của chúng không thể so sánh được với các biểu đồ được dựng thủ công, nhất là density plot của package ggridges.
7 Kết luận
Vào năm 2016 có người khẳng định lạc quan rằng “Sớm muộn gì tất cả người dùng R đều sẽ theo trường phái Bayes”. Điều này có thể tin được nếu ta biết rằng hiện nay phương pháp suy diễn Bayes đã trở nên dễ dàng và phổ biến hơn bao giờ hết trong R. Ngôn ngữ STAN trở thành nền tảng chung cho thống kê Bayes, ít nhất 3 package trong R cho phép kết nối với STAN, trong đó brms là công cụ đơn giản và mạnh nhất. Những packages như sjstats và sjPlot đã cung cấp cho người dùng mảnh ghép cuối cùng còn thiếu - khai thác mô hình Bayes. Với những hàm được giới thiệu trong bài, ta hoàn toàn có thể hy vọng về một tương lai rất rộng mở cho trường phái Bayes, từng bước thay thế cho những phương pháp kiểm định cổ điển khác.