BAYES thay thế kiểm định Chi2

Sử dụng ngôn ngữ STAN

1 Giới thiệu

Thân chào các bạn, đây là bài thực hành thứ ba trong dự án Bayes for Vietnam, được reboot lại từ project cùng tên năm 2016. Mục tiêu của dự án này là phổ biến phương pháp thống kê theo trường phái Bayes cho các bạn bác sĩ và sinh viên y khoa, nhằm thay thế cho những công cụ thống kê truyền thống.

Qua 2 bài trước, chắc hẳn các bạn đã được làm quen với quy trình phân tích và suy diễn thống kê theo trường phái Bayes và nhận ra những điểm khác biệt của Bayes so với quy trình cổ điển. Nguyên lý của phương pháp này dựa vào định lý Bayes, được diễn đạt bằng ngôn ngữ mô hình, xác suất điều kiện. Suy diễn thống kê theo Bayes liên kết giữa 4 phần: Dữ liệu quan sát thực tế, giả thuyết (phân phối) tiền định Prior, hàm likelihood và kết quả phân phối hậu địn h. Quy trình Bayes đòi hỏi người thực hành nhìn thấu suốt vấn đề mà họ phải đối diện đến từng chi tiết nhỏ nhất: Họ phải nhận diện ra những biến số trong bài toán, đâu là biến kết quả, đâu là tham số. Tiếp theo, người dùng phải chuyển được bài toán sang ngôn ngữ mô hình và chọn một mô hình phù hợp. Từ mô hình này, ta lại phải xác định các tham số, hàm likelihood, phân phối tiền định cho tham số. Khi đã phác thảo xong mô hình, ta viết code để đưa mô hình vào một Sampler để lấy mẫu chuỗi MCMC cho tham số và các trị số thống kê khác. Các chuỗi MCMC này chứa phân phối hậu định, cũng là kết quả của phân tích. Cuối cùng, ta suy diễn thống kê trên phân phối hậu định, bằng nhiều cách khác nhau bao gồm Bayes factor, khoảng vô nghĩa thực dụng (ROPE) hay ngưỡng so sánh (Comp.Val).

Tại một số vị trí nào đó, Bayes và thống kê cổ điển có thể giao nhau, thí dụ từ mô hình Bayes, ta có thể rút ra phân phối hậu định cho các trị số thống kê, bao gồm cả Effect-size. Như vậy khác biệt quan trọng nhất giữa Bayes và các kiểm định thống kê quy ước, đó là có sự hòa hợp giữa prior và posterior, đó là xác suất điều kiện, likelihood và mô hình thay cho những công thức và trị số thống kê mang tính đại diện và xấp xỉ (như Chi2 cho Pearson Chi2 test, t cho Student t-test, F cho ANOVA , r cho phân tích tương quan vv…). Khác biệt thứ hai nằm ở cách diễn giải kết quả. Thực ra, P_vale hoàn toàn có thể được tính dựa vào phân phối hậu định, nhưng bản thân phân phối hậu định đã chứa nhiều thông tin đến mức ta không cần phải dùng đến phản nghiệm, Null hypothesis testings và p_value nữa.

Sau khi thay thế Student-t test trong bài thứ 2, bài này Nhi sẽ hướng dẫn các bạn dùng Bayes để thay thế cho một công cụ phổ biến khác đó là Chi2 test trong phân tích tần số, tỉ lệ của bảng chéo

Phân tích bảng chéo thực ra là một bài toán xác suất rất đơn giản, đến mức sinh viên Y khoa năm thứ nhất có thể giải quyết dễ dàng; tại sao ta phải dùng đến Bayes ? Nhi hy vọng sau bài thực hành bạn sẽ có câu trả lời cho mình.

2 Bài toán bảng chéo và 3 cách tư duy

Bộ số liệu này có nguồn gốc từ 1 nghiên cứu thử nghiệm lâm sàng của Edwards Koch năm 1988, với mục tiêu khảo sát hiệu quả điều trị bệnh Viêm khớp của một loại thuốc (nhóm Treatment), so với nhóm sử dụng giả dược (placebo). Hiệu quả điều trị được đánh giá như một biến định danh với 3 bậc giá trị : Cải thiện rõ rệt (marked), cải thiện tương đối (Some) và không có cải thiện (None) triệu chứng của bệnh.

library(tidyverse)

df=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/vcd/Arthritis.csv")

crosstab=xtabs(~Treatment+Improved,data=df)

crosstab

##          Improved
## Treatment Marked None Some
##   Placebo      7   29    7
##   Treated     21   13    7

library(vcd)

mosaic(crosstab,shade=TRUE,legend=TRUE,gp = shading_max)

Như vậy bài toán này có 2 biến số (yếu tố) : Treatment và Outcome, cả 2 đều là biến định danh với 2 hoặc 3 bậc giá trị. Với dữ liệu như vậy, ta có thể phát biểu câu hỏi nghiên cứu (giả thuyết H0/H1) theo nhiều cách khác nhau; và chú ý rằng mỗi cách như vậy sẽ gợi ý cho chúng ta sử dụng một hình thức/công cụ thống kê khác nhau:

Cách thứ 1

Giả thuyết cần chứng minh là : Có mối liên hệ giữa yếu tố phân nhóm điều trị và mức độ cải thiện triệu chứng; Ghi chú: liên hệ ở đây được hiểu ngầm – là về tần số / tỉ lệ phân bố. Cách tư duy này sẽ dẫn chúng ta đến Chisquared test của Pearson với H0: 2 yếu tố Điều trị và Cải thiện triệu chứng là độc lập với nhau (không có liên hệ giữa chúng), độ tự do cho Chi2= 2. Chúng ta sẽ bàn về Chi2test trong phần tiếp theo.

Cách thứ 2

Giả thuyết cần chứng minh sẽ là: Có sự khác biệt về tỉ lệ phân bố của bệnh nhân Có cải thiện triệu chứng giữa 2 phân nhóm : Dùng thuốc và dùng Giả dược. Cách tư duy này sẽ dẫn chúng ta đến bài toán so sánh tỉ lệ điều trị thành công giữa 2 phân nhóm (Treatment vs Placebo), ghi chú là tỉ lệ này chính là tỉ số giữa số bệnh nhân có cải thiện chia cho tổng số bệnh nhân ở từng phân nhóm điều trị. Bài toán này dẫn chúng ta đến một mô hình : số bệnh nhân điều trị thành công là biến kết quả, được mô tả bằng phân phối Binomial (theta, ntrial), với ntrial = tổng số bệnh nhân ở mỗi phân nhóm điều trị, và theta là tỉ lệ điều trị thành công, tham số của mô hình. Như vậy bài toán này gần giống với bài toán so sánh trung bình giữa 2 phân nhóm (kiểm định t), chỉ thay phân phối Student t bằng phân phối binomial và tham số Mu bằng tham số Theta. Ghi chú: để thực hiện phân tích, ta phải gộp 2 phân loại: Cải thiện đáng kể và cải thiện tương đối thành một.

Cách thứ 3:

Giả thuyết được phát biểu như sau: 1 bệnh nhân bất kì nếu được điều trị bằng thuốc X sẽ có xác suất cải thiện triệu chứng cao hơn so với khi dùng placebo. Outcome của bài toán này là 1 xác suất cho 1 cá thể. Chắc các bạn cũng nhận ra, đây là 1 mô hình logistic hoặc mô hình hồi quy đa giá trị (tùy theo biến kết quả được gộp thành 1 biến nhị giá :Có/không có cải thiện hay giữ nguyên 3 bậc giá trị. Cách thứ 3 này hoàn toàn khác với 2 cách trên, và không có liên quan đến phân tích bảng chéo/ tần số, nên chúng ta sẽ không xét về nó trong bài.

Tuy cách 1 và 2 đều biểu diễn dữ liệu bằng bảng chéo, nhưng bản chất câu hỏi và cách giải quyết rất khác nhau. Như đã trình bày, kiểm định Chi2 là giải pháp cho câu hỏi chứng minh về tính độc lập/liên hệ giữa 2 yếu tố (biến định danh). Do đó, cần phải nêu rõ 1 điều ở đây: công dụng của Chisquare test không phải là để so sánh về tỉ lệ phân phối giữa các phân nhóm. Rất nhiều bạn sinh viên nhầm lẫn ở điểm này, nên các bạn thường viết trong báo cáo/luận văn như sau : tỉ lệ bệnh nhân giữa 2 phân nhóm được SO SÁNH bằng test Chi2. Đây là 1 phát biểu sai vì khi làm Chi2 test ta không so sánh cái gì cả. Mặt khác, bạn sẽ thấy rằng: ngay cả khi chúng ta trả lời thành công cho câu hỏi thứ 1 (test chi2 cho ra kết quả ý nghĩa, p<0.05), chứng minh được có mối liên hệ giữa yếu tố Nhóm điều trị và Yếu tố Cải thiện triệu chứng, ta vẫn chưa thể kết luận gì về Hiệu ứng của Thuốc so với Placebo cả ! (giả sử ta hoán chuyển tên 2 phân nhóm Placebo và điều trị, bạn vẫn sẽ có 1 kết quả test Chi2 dương tính như thường nhưng ý nghĩa lâm sàng khác hẳn !); do đó Cách đặt vấn đề thứ 2 và thứ 3 tốt hơn nhiều so với cách thứ 1, vì bạn có trong tay bằng chứng định lượng (xác suất khỏi bệnh, khác biệt tần số bệnh nhân khỏi bệnh, odss-ratio) giữa 2 phân nhóm.

3 Nhắc lại về Chi2 test cổ điển

Chisquared test là gì ? Đây là một kiểm định rất cổ xưa, do Karl Pearson tạo ra từ năm 1900, công dụng của nó là kiểm tra mối liên hệ giữa Hàng và Cột trong một bảng chéo.

Giả thuyết H0 của Chi2 test là: hàng và cột trong bảng chéo độc lập với nhau = không có liên hệ giữa 2 yếu tố ta đang xét. Giả thuyết H1 là: có mối liên hệ giữa hàng và cột trong bảng chéo

Chisquared test có 2 giả định, thứ nhất, mọi quan sát trong bảng chéo là độc lập với nhau (1 cá thể chỉ xuất hiện 1 lần ở hàng/cột); thứ 2, không quá 20% tần số giả định có giá trị < 5. Giả định thứ 2 này khá phiền toái, có thể hiểu là với 1 bảng 2x2, và tần số ở 1 ô nào đó <5, ta không thể dùng test Chi2 được nữa, vì khi đó đường cong mật độ phân phối của Chi2 sẽ có hình dạng bất đối xứng, không chính xác để tính xác suất cho trị số p. Giả định này có thể được kiểm tra đơn giản bằng cách xem kết quả tần số giả định trong bảng chéo.

Trị số thống kê của test có bản chất là tổng bình phương sai số (residual) giữa Tần số quan sát được trên thực tế trong bảng chéo (observed) và tần số giả định (expected) nếu giả thuyết H0 là đúng.

Nếu giả thuyết H0 đúng, trị số thống kê này tuân theo quy luật phân phối Chisquared với độ tự do =df; do đó nó được đặt tên là Chisquared. Thực ra phân phối Chi2 rất phổ biến trong kiểm định thống kê; có vài chục loại kiểm định mà trị số thống kê được mô tả bằng phân phối Chi2.

Tần số giả định:

\[expected = \frac{rowsum * colsum}{total}\]

Độ tự do df của Chi2 test được tính như sau:

\[df = (nrow - 1)(ncol - 1)\]

Pearson residual được tính như sau:

\[residual = \frac{observed - expected}{\sqrt{expected}}\]

\[\chi^2 = \sum{\frac{(observed - expected)^2}{expected}} = \sum{residuals}\]

Chisquare là một phân phối liên tục, nó có 1 tham số kiểu hình là Nu=df. Minh họa: hình dáng của phân phối Chisquared ở độ tự do từ 2 đến 15; khi df tăng dần, Chi2 trở nên đối xứng và gần với phân phối Gaussian.

# Chi2 

sim=data.frame(Val=rep(NA,500),nu=rep(NA,500))
temp=sim

for (i in (2:15)){
  temp$Val=rchisq(n=500,df=i)
  temp$nu=i
  sim=rbind(temp,sim)
}

sim$nu=as.factor(sim$nu)

sim%>%na.omit%>%ggplot(aes(x=Val,fill=nu,col=nu))+
  geom_density(alpha=0.1)+
  theme_bw()+ggtitle("Chisquared distribution")

Chisquared test có một số biến thể, bao gồm test Exact của Fisher, phương pháp mô phỏng Monte Carlo. Effect-size của chi2 test là V (theo Cramer). Hàm sau đây sẽ thực hiện cùng lúc cả 3 loại test và tính Cramer’s V cho 1 bảng chéo:

chisquare=function(crosstab,B,level){
  test1=chisq.test(crosstab)
  test2=chisq.test(crosstab,simulate.p.value = T,B=B)
  test3=test3=fisher.test(crosstab)
  p=Monte_Carlo_pval=test2$p.value[[1]]
  print(rbind(Chisqr=test1$statistic[[1]],
              df=test1$parameter[[1]],
              Chi2p_val=test1$p.value[[1]],
              Monte_Carlo_pval=p,
              Fisher_exact_pval=test3$p.value[[1]],
              Cramer_V=lsr::cramersV(crosstab)))
        }

chisquare(crosstab,10000)

##                           [,1]
## Chisqr            13.055019853
## df                 2.000000000
## Chi2p_val          0.001462643
## Monte_Carlo_pval   0.001699830
## Fisher_exact_pval  0.001393195
## Cramer_V           0.394229505

Khi dùng hàm summary cho 1 bảng chéo, R tự động làm test Chi2

crosstab%>%summary()

## Call: xtabs(formula = ~Treatment + Improved, data = df)
## Number of cases in table: 84 
## Number of factors: 2 
## Test for independence of all factors:
##  Chisq = 13.055, df = 2, p-value = 0.001463

Hàm chisq.test cho phép làm Chi2 test và trích xuất các trị số bên trong:

test1=chisq.test(crosstab)

test1

## 
##  Pearson's Chi-squared test
## 
## data:  crosstab
## X-squared = 13.055, df = 2, p-value = 0.001463

Thí dụ: Observed values:

test1$observed

##          Improved
## Treatment Marked None Some
##   Placebo      7   29    7
##   Treated     21   13    7

Expected values:

test1$expected

##          Improved
## Treatment   Marked None     Some
##   Placebo 14.33333 21.5 7.166667
##   Treated 13.66667 20.5 6.833333

Pearson residual

test1$residuals

##          Improved
## Treatment      Marked        None        Some
##   Placebo -1.93699199  1.61749160 -0.06225728
##   Treated  1.98367320 -1.65647289  0.06375767

4 Chisquared test BAYES

Ta có thể mô phỏng Chi2 test theo phương pháp BAYES như sau:

Ta muốn xác định phân phối hậu định của trị số Chisquared và Effect-szie của nó là Cramer’s V bằng một mô hình Bayes có likelihood function ước tính Chisquared theo phân phối Chi2 với tham số là Nu. Tiền định của mô hình là thal số Nu (độ tự do của phân phối Chi2) có phân phối uniform (0.01, 100). Dữ liệu đầu vào của mô hình là 2 vector chứa tần số quan sát được (observed) và tần số giả định (expect) trong bảng chéo.

expect=test1$expected%>%as.vector()
obs=test1$observed%>%as.vector()

Đây là nội dung mô hình STAN

# Bayesian Chi2 test

library(rstan)

stanmodelcode ="

data{
int<lower=4> ncell;
vector[ncell] observed;
vector[ncell] expect;
}

transformed data {

real residual[ncell];
real chisqr;
real total;

for (i in 1:ncell) {
residual[i]=(observed[i]-expect[i])^2/expect[i];}

chisqr=sum(residual);
total = sum(observed);

}

parameters {
real nu;
}

model {

// khai bao tien dinh (prior)

nu ~ uniform(0.01,100);

// likelihood

chisqr ~ chi_square(nu);

}

generated quantities {
real cramersV;
real chi2;

chi2=chi_square_rng(nu);
cramersV=sqrt(chi2/total);
}
"

Mô tả nội dung model STAN:

1. Block data: cần 3 thành phần: ncell= số ô trong bảng chéo, 2 vector observed và expect.

2. Block transformed data: khai báo vector residual là 1 số thực, chisqr là 1 số thực, total là 1 số thực. residual được tính cho từng ô trong vector, dựa vào expect và observed vector. Có residual, ta dễ dàng tính được Chi2 = sum(residual). Total = tổng tất cả giá trị observed.

3. Model Bayes có likelihood chisqr ~ chi_square(nu), nu là một số thực, chỉ độ tự do của Chi2, nu có prior là phân phối uniform.

4. Block generated quantities: ta lấy mẫu ngẫu nhiên cho biến chi2 là 1 số thực,dựa vào phân phối hậu định của nu, rồi từ đó tính Cramer’s V là effectsize, cũng là 1 số thực.

Khi thi hành mô hình BAYES này, chúng ta sẽ có phân phối hậu định của Nu, Chi2 và Cramers’V.

data=list(ncell=6,observed=obs,expected=expect)

fit = stan(model_code=stanmodelcode, data=data, iter=1000, warmup=200, chains=1,cores=4, thin=10)

## In file included from C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from filef081bc235d1.cpp:8:
## C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
##  #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##  ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
## 
## SAMPLING FOR MODEL '16c3f04e4fe9d979aca9ed660448df9f' NOW (CHAIN 1).
## Rejecting initial value:
##   Error evaluating the log probability at the initial value.
## Exception: chi_square_lpdf: Degrees of freedom parameter is -0.80088, but must be > 0!  (in 'modelf0894ab16_16c3f04e4fe9d979aca9ed660448df9f' at line 35)
## 
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 201 / 1000 [ 20%]  (Sampling)
## Iteration: 300 / 1000 [ 30%]  (Sampling)
## Iteration: 400 / 1000 [ 40%]  (Sampling)
## Iteration: 500 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.01 seconds (Warm-up)
##                0.01 seconds (Sampling)
##                0.02 seconds (Total)

summary(fit)

## $summary
##                mean    se_mean        sd      2.5%       25%        50%
## nu       15.2526184 0.68121374 5.4843252 7.5471615 11.373353 14.3500838
## cramersV  0.4071654 0.01341651 0.1105121 0.1922191  0.349046  0.4089721
## chi2     14.9388936 0.94685718 7.6189792 3.1081116 10.237335 14.0497266
## lp__      3.4585111 0.11744526 0.8783617 0.4998120  3.376739  3.7322409
##                 75%      97.5%    n_eff     Rhat
## nu       18.8198123 29.6969801 64.81561 1.007013
## cramersV  0.4624003  0.6452856 67.84858 1.028366
## chi2     17.9605729 34.9771457 64.74774 1.039745
## lp__      3.9522331  3.9831383 55.93399 1.021946
## 
## $c_summary
## , , chains = chain:1
## 
##           stats
## parameter        mean        sd      2.5%       25%        50%        75%
##   nu       15.2526184 5.4843252 7.5471615 11.373353 14.3500838 18.8198123
##   cramersV  0.4071654 0.1105121 0.1922191  0.349046  0.4089721  0.4624003
##   chi2     14.9388936 7.6189792 3.1081116 10.237335 14.0497266 17.9605729
##   lp__      3.4585111 0.8783617 0.4998120  3.376739  3.7322409  3.9522331
##           stats
## parameter       97.5%
##   nu       29.6969801
##   cramersV  0.6452856
##   chi2     34.9771457
##   lp__      3.9831383

Ta có thể thấy rằng kết quả Median của Chi2 và Cramer’s V tương đương với test Chi2 frequentist

Ta có thể diễn giải kết quả trực tiếp dựa vào Cramers’V, hoặc dùng phân phối hậu định của Chi2 để ước tính một giá trị p_value nhằm kiểm tra giả thuyết H0. Giá trị p này = 0.0007, rất gần với p_value theo phương pháp mô phỏng Monte Carlomaé ta tìm ra ở trên.

chi2=extract(fit,pars="chi2")%>%.[[1]]

p=rep(NA,length(chi2))

for (i in 1:length(chi2)){
  p[i]=pchisq(chi2[i],2,lower.tail=F)}

median(p)

## [1] 0.0008897486

Như vậy, cách làm này không cho ta biết thêm thông tin nào mới so với test Chi2 cổ điển.

5 So sánh tỉ lệ bằng mô hình Binomial BAYES

Một cách tổng quát, suy diễn Bayes dựa trên định lý Bayes như sau :

\[ p(\theta |y,x)\propto p(y|\theta ,x) p(\theta ,x) \] theo đó, phân phối hậu định của một tham số Theta (khi ta có trong tay dữ liệu biến kết quả Y và hằng số X là một matrix các predictors trong mô hình) tỉ lệ với tích của hàm likelihood (xác suất có điều kiện của y khi có theta và x) và phân phối tiền định (prior= giả thuyết về phân phối của theta).

Ta có thể biến cách tiếp cận thứ 2 (So 2 sánh tỉ lệ thành công giữa 2 phân nhóm Điều trị vs Placebo) thành một mô hình Bayes như sau:

Dữ liệu đầu vào gồm có:

1 vector Proportion có 2 giá trị là tần số điều trị thành công tương ứng cho 2 nhóm

1 vector Total với 2 giá trị tổng số bệnh nhân ở mỗi phân nhóm

Ở mỗi phân nhóm, biến kết quả của mô hình có thể xem như là số lần thành công (Proportion[i]) rút ra được sau n[i] lần thử từ một hộp kín chứa n[i] phần tử (=Total[i]).

Như vậy quy luật đằng sau thử nghiệm này chính là phân phối Binomial. Như ta biết, phân phối Binomial có 2 tham số là số lần thử (n trial) và tỉ lệ thành công (Theta). Số lần thử ở đây ta đã biết, chính là tổng số bệnh nhân ở mỗi phân nhóm, ta chỉ còn phải xác định phân phối hậu định của tham số Theta.

# Binomial

sim=data.frame(Val=rep(NA,100),theta=rep(NA,100))
temp=sim
theta=c(0.05,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)

for (i in (1:10)){
  temp$Val=rbinom(n=100,prob=theta[i],size=41)
  temp$theta=theta[i]
  sim=rbind(temp,sim)
}

sim$theta=as.factor(sim$theta)

sim%>%na.omit%>%ggplot(aes(x=Val,fill=theta,col=theta))+
  geom_density(alpha=0.2)+
  theme_bw()+ggtitle("Binomial distribution / 41 trials")

Giả thuyết tiền định cho tham số Theta là gì ?: Do theta là 1 xác suất có giá trị từ 0.0 đến 1.0, sẽ hợp lý nếu ta giả định rằng theta được mô tả bằng phân phối Beta. Nếu bạn không có thông tin nào về prior, bạn có thể chọn giải pháp an toàn khi đặt phân phối Beta(shape1, shape2) có hình dạng “phẳng” tức gần giống 1 phân phối Uniform(0,1), bằng cách dùng 2 tham số shape bằng nhau = 1.

Ở đây Nhi dùng 1 prior do Jeffrey đề nghị là Beta (0.5, 0.5). Một prior khác bạn có thể thử là prior của Haldane : Beta (0.0001, 0.0001).

#Beta Distribution

sim=data.frame(Val=rep(NA,1000),shape=rep(NA,500))
temp=sim
para=c(0.0001,0.5,1,2,3,4)

for (i in (1:6)){
  temp$Val=rbeta(n=1000,shape1=para[i],shape2=para[i])
  temp$shape=para[i]
  sim=rbind(temp,sim)
}

sim$shape=as.factor(sim$shape)

sim%>%na.omit%>%ggplot(aes(x=Val,fill=shape,col=shape))+
  geom_density(alpha=0.2)+
  theme_bw()+ggtitle("Beta distribution")

Sau đây là nội dung model STAN

# Proportional test BAYES


stanmodelcode ="
data {
int <lower=2> ng;
int proportion[ng];
int total[ng];
}

transformed data {
}

parameters {
real theta[ng];
}

model {

// khai bao tien dinh (prior)

theta ~ beta(0.5,0.5);

// likelihood
for (i in 1:ng) {
proportion[i] ~ binomial(total[i],theta[i]);
}

}

generated quantities {

int pred[ng];

for (i in 1:ng){
pred[i] = binomial_rng(total[i],theta[i]);
}



}

"

1. Block data gồm có 3 thành phần: ng là 1 số nguyên >2, chỉ số phân nhóm cần so sánh, trong thí dụ này ng=2; proportion là 1 vector số nguyên có độ dài ng, chứa tần số bệnh nhân có đáp ứng điều trị ở 2 phân nhóm placebo và treatment; total cũng là 1 vector số nguyên, chứa tổng số bệnh nhân ở mỗi phân nhóm placebo và treatmen

2. Block transformed data không cần thiết

3. block parameter khai báo 1 tham số theta là 1 vector số thực, có độ dài ng

4. Block model gồm có: prior của theta là beta(0.5, 0.5);

hàm likelihood của biến kết quả là proportion[i] ~ binomial(total[i],theta[i]);

5. Block generated quantities:

Tính giá trị dự báo pred1 và pred2 theo mô hình binomial : pred[i] = binomial_rng(total[i],theta[i]);

Ta tạo 2 vector: positive = số bệnh nhân có cải thiện triệu chứng, và rowsum = tổng số bệnh nhân ở 2 phân nhóm điều trị vs placebo

marked=crosstab[,1]%>%as.vector()
some=crosstab[,3]%>%as.vector()
none=crosstab[,2]%>%as.vector()

positive=marked+some
rowsum=marked+some+none

Và ta thi hành mô hình STAN

data=list(ng=2,proportion=positive,total=rowsum)

fit = stan(model_code=stanmodelcode, data=data, iter=1000, warmup=200, chains=1,cores=4, thin=1)

## In file included from C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.4/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from filef081f962d39.cpp:8:
## C:/Users/Admin/Documents/R/win-library/3.4/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
##  #  define BOOST_NO_CXX11_RVALUE_REFERENCES
##  ^
## <command-line>:0:0: note: this is the location of the previous definition
## cc1plus.exe: warning: unrecognized command line option "-Wno-ignored-attributes"
## 
## SAMPLING FOR MODEL '04fe3b5c2c1476993474804c7e83a5a1' 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 / 1000 [  0%]  (Warmup)
## Iteration: 100 / 1000 [ 10%]  (Warmup)
## Iteration: 200 / 1000 [ 20%]  (Warmup)
## Iteration: 201 / 1000 [ 20%]  (Sampling)
## Iteration: 300 / 1000 [ 30%]  (Sampling)
## Iteration: 400 / 1000 [ 40%]  (Sampling)
## Iteration: 500 / 1000 [ 50%]  (Sampling)
## Iteration: 600 / 1000 [ 60%]  (Sampling)
## Iteration: 700 / 1000 [ 70%]  (Sampling)
## Iteration: 800 / 1000 [ 80%]  (Sampling)
## Iteration: 900 / 1000 [ 90%]  (Sampling)
## Iteration: 1000 / 1000 [100%]  (Sampling)
## 
##  Elapsed Time: 0.01 seconds (Warm-up)
##                0.02 seconds (Sampling)
##                0.03 seconds (Total)

Sau khi thi hành mô hình BAYES này trong STAN, ta sẽ có phân phối hậu định của 2 tham số Theta1, tương ứng với tỉ lệ điều trị thành công ở mỗi phân nhóm (theta quyết định đặc tính của phân phối binomial, và từ đó quyết định tần số bệnh nhân tại mỗi ô trong bảng 2x2). Theta1 cho nhóm Placebo, Theta2 cho nhóm Treatment. Ta cũng có phân phối hậu định của tần số bệnh nhân (có cải thiện triệu chứng) được ước tính : Pred1 và Pred2.

summary(fit)

## $summary
##                 mean     se_mean         sd        2.5%         25%
## theta[1]   0.3309912 0.003078324 0.06609964   0.2102120   0.2858369
## theta[2]   0.6808666 0.002788962 0.07180417   0.5412657   0.6333668
## pred[1]   14.2150000 0.173538683 4.21704925   7.0000000  11.0000000
## pred[2]   27.9112500 0.156339165 4.23713524  19.0000000  25.0000000
## lp__     -52.1387069 0.051342673 0.92616717 -54.7483711 -52.4387371
##                  50%         75%       97.5%    n_eff      Rhat
## theta[1]   0.3305032   0.3724852   0.4772503 461.0727 0.9988131
## theta[2]   0.6833480   0.7316253   0.8138947 662.8483 0.9993492
## pred[1]   14.0000000  17.0000000  23.0000000 590.5066 0.9999734
## pred[2]   28.0000000  31.0000000  36.0000000 734.5292 1.0000668
## lp__     -51.8454698 -51.4643321 -51.2403454 325.4032 1.0001991
## 
## $c_summary
## , , chains = chain:1
## 
##           stats
## parameter         mean         sd        2.5%         25%         50%
##   theta[1]   0.3309912 0.06609964   0.2102120   0.2858369   0.3305032
##   theta[2]   0.6808666 0.07180417   0.5412657   0.6333668   0.6833480
##   pred[1]   14.2150000 4.21704925   7.0000000  11.0000000  14.0000000
##   pred[2]   27.9112500 4.23713524  19.0000000  25.0000000  28.0000000
##   lp__     -52.1387069 0.92616717 -54.7483711 -52.4387371 -51.8454698
##           stats
## parameter          75%       97.5%
##   theta[1]   0.3724852   0.4772503
##   theta[2]   0.7316253   0.8138947
##   pred[1]   17.0000000  23.0000000
##   pred[2]   31.0000000  36.0000000
##   lp__     -51.4643321 -51.2403454

Kết quả sơ bộ cho thấy tỉ lệ bệnh nhân có cải thiện triệu chứng cao hơn hẳn ở nhóm trị liệu so với nhóm placebo (theta2=0.68 so với theta1=0.328).

Từ 2 chuỗi MCMC chứa phân phối hậu định của theta1, theta2, ta có thể tạo ra một vector chứa phân phối hậu định của Khác biệt xác suất điều trị thành công giữa 2 phân nhóm: Difprop = Theta2 – Theta1. Suy diễn thống kê sẽ dựa vào vector Difprop này, sử dụng các phương pháp mà Nhi đã trình bày trong bài trước là Comp.Val, ROPE và Bayes Factor.

s=as.data.frame(fit)%>%.[,-5]

colnames(s)=c("Theta_Placebo","Theta_Treated","Pred_Placebo","Pred_Treated")

s$Difprop=s$Theta_Treated-s$Theta_Placebo

Hmisc::describe(s$Difprop)

## s$Difprop 
##        n  missing distinct     Info     Mean      Gmd      .05      .10 
##      800        0      793        1   0.3499   0.1092   0.1813   0.2184 
##      .25      .50      .75      .90      .95 
##   0.2868   0.3536   0.4170   0.4727   0.5071 
## 
## lowest : 0.06469176 0.08309505 0.09891270 0.09959744 0.10642911
## highest: 0.59264643 0.59302557 0.59726109 0.62291389 0.62578527

s$pseudoGroup=factor(rep(c(1:20),e=nrow(s)/20))
s$Iter=as.numeric(rownames(s))

Phân phối hậu định của khác biệt về tỉ lệ đáp ứng điều trị là 0.173 -0.519, median khác biệt là 0.3554

Ta có thể vẽ hình phân phối hậu định cho Theta1,Theta2,Difprop như sau:

p1=s%>%ggplot(aes(x=Theta_Treated))+
  geom_density(alpha=0.5,col="red4",size=1) +
  geom_point(y=0, alpha=.01, size=2,col="red3") +
  geom_line(aes(group=pseudoGroup, color=pseudoGroup), stat='density', alpha=.5,show.legend = F)+
  scale_x_continuous("Theta_Treated")+
  theme_bw()


p2=s%>%ggplot(aes(x=Theta_Placebo))+
  geom_density(alpha=0.5,col="blue4",size=1) +
  geom_point(y=0, alpha=.01, size=2,col="blue3") +
  geom_line(aes(group=pseudoGroup, color=pseudoGroup), stat='density', alpha=.5,show.legend = F)+
  scale_x_continuous("Theta_Placebo")+
  theme_bw()

p3=s%>%ggplot(aes(x=Difprop))+
  geom_density(alpha=0.5,col="purple",size=1) +
  geom_point(y=0, alpha=.01, size=2,col="purple") +
  geom_line(aes(group=pseudoGroup, color=pseudoGroup), stat='density', alpha=.5,show.legend = F)+
  scale_x_continuous("Proportion Difference")+
  theme_bw()

p11=s%>%ggplot(aes(x=Iter,y=Theta_Treated))+
  geom_path(alpha=.7,show.legend = F,col="red")+
  scale_y_continuous("Theta_Treated")+
  theme_bw()

p21=s%>%ggplot(aes(x=Iter,y=Theta_Placebo))+
  geom_path(alpha=.7,show.legend = F,col="blue")+
  scale_y_continuous("Theta_Placebo")+
  theme_bw()

p31=s%>%ggplot(aes(x=Iter,y=Theta_Placebo))+
  geom_path(alpha=.7,show.legend = F,col="violet")+
  scale_y_continuous("Difprop")+
  theme_bw()

gridExtra::grid.arrange(p11,p1,p21,p2,p31,p3,ncol=2)

Tiếp theo, Nhi sẽ giới thiệu với các bạn cả 3 cách suy diễn thống kê theo Bayes, đó là Bayes Factor, Ngưỡng so sánh (CompVal) và khoảng vô nghĩa thực dụng (ROPE) :

6 Bayes Factor

Bản chất của Bayes Factor là tỉ trọng chứng cứ ủng hộ cho một giả thuyết H1 so với giả thuyết H0. Giả thuyết này được xác định bằng 1 ngưỡng giá trị đặc biệt. Bayes factor cho phép kiểm chứng (xác nhận hay loại trừ) giả thuyết H1 có liên quan tới ngưỡng giá trị này.

Ranh giới của H0 và H1 là BF=1, BF càng cao (>1), thì độ khả tín của kết luận (H1) càng lớn, và ngược lại, BF càng thấp (<1) thì chứng cứ ủng hộ cho H0 càng cao. Thông thường khi BF > 10 đã là « strong evidence », BF >30 thì tỉ trọng chứng cứ đã là rất mạnh ; còn khi BF > 100 thì có thể gần như khẳng định chắc chắn H1 là đúng.

BF có thể tính trực tiếp nhờ hàm hypothesis trong package brms, nó nhận 1 input là mô hình STAN (bất kì mô hình Bayes nào do STAN tạo ra, dù là brms, rstan hay rstanarm), và xuất kết quả là Evidence ratio chính là BF.

7 ROPE và CompVal

Ngoài BF, tác giả John Kruschke đề xuất diễn giải kết quả phân tích Bayes dựa vào tỉ trọng phân phối hậu nghiệm sau khi phân chia bằng 1 khoảng thực dụng (ROPE) hoặc 1 ngưỡng (CompVal). Mục tiêu là xác định mật độ (%) của phân phối hậu nghiệm nằm trong, hay ngoài (cao hơn, thấp hơn) ROPE, hoặc cao/thấp hơn ngưỡng CompVal.

Ngưỡng so sánh : CompVal có bản chất là 1 ngưỡng giá trị nhất định, cho phép chia phân phối hậu nghiệm thành 2 vùng: Cao hơn/Thấp hơn. Công dụng của nó là khẳng định 1 giả thuyết H1 hoặc loại trừ giả thuyết H0 có hàm ý so sánh với 1 ngưỡng ý nghĩa đặc biệt.

ROPE (Region of practice equivalence) có bản chất là 1 khoảng giá trị trên thang đo, giới hạn bởi 2 ngưỡng trên/dưới. ROPE cho phép chia phân phối hậu nghiệm ra 3 vùng:

Vùng trong ROPE xem như vô nghĩa

Ngoài ROPE : gồm Vùng cao / Vùng Thấp

ROPE cho phép khẳng định 1 giả thuyết H1 hoặc loại trừ 1 giả thuyết H0 liên quan tới 1,2 hoặc 3 vùng với ý nghĩa thực dụng tùy chọn. Do đó ROPE tiện ích hơn nhiều so với BF và Null hypothesis testing, vì nó cho phép trả lời cả 3 câu hỏi : kết quả cao hay thấp so với ngưỡng so sánh ? kết quả có ý nghĩa lâm sàng hay không ? Mức độ khả tín là bao nhiêu ?

Phương pháp J.Kruschke áp dụng trên output thuộc lớp MCMC và dựa vào 1 vài hàm do tác giả viết, Nhi chỉ cải biên lại đôi chút.

#Hypothesis testing

# Kruschke codes

HDIF= function( sampleVec,credMass=0.975 ) {
  sortedPts = sort( sampleVec )
  ciIdxInc = ceiling( credMass * length( sortedPts ) )
  nCIs = length( sortedPts ) - ciIdxInc
  ciWidth = rep( 0 , nCIs )
  for ( i in 1:nCIs ) {
    ciWidth[ i ] = sortedPts[ i + ciIdxInc ] - sortedPts[ i ]
  }
  HDImin = sortedPts[ which.min( ciWidth ) ]
  HDImax = sortedPts[ which.min( ciWidth ) + ciIdxInc ]
  HDIlim = c( HDImin , HDImax )
  return( HDIlim )
}

SUMK=function(paramSampleVec,compVal=NULL , ROPE=NULL , credMass=0.975) {
  meanParam = mean( paramSampleVec )
  medianParam = median( paramSampleVec )
  dres = density( paramSampleVec )
  modeParam = dres$x[which.max(dres$y)]
  hdiLim = HDIF( paramSampleVec , credMass=credMass )
  if ( !is.null(compVal) ) {
    pcgtCompVal = ( 100 * sum( paramSampleVec > compVal ) 
                    / length( paramSampleVec ) )
  } else {
    compVal=NA
    pcgtCompVal=NA
  }
  if ( !is.null(ROPE) ) {
    pcltRope = ( 100 * sum( paramSampleVec < ROPE[1] ) 
                 / length( paramSampleVec ) )
    pcgtRope = ( 100 * sum( paramSampleVec > ROPE[2] ) 
                 / length( paramSampleVec ) )
    pcinRope = 100-(pcltRope+pcgtRope)
  } else { 
    ROPE = c(NA,NA)
    pcltRope=NA 
    pcgtRope=NA 
    pcinRope=NA 
  }  
  return( c( Mean=meanParam , Median=medianParam , Mode=modeParam , 
             HDIlevel=credMass , LL=hdiLim[1] , UL=hdiLim[2] , 
             CompVal=compVal , PcntGtCompVal=pcgtCompVal , 
             ROPElow=ROPE[1] , ROPEhigh=ROPE[2] ,
             PcntLtROPE=pcltRope , PcntInROPE=pcinRope , PcntGtROPE=pcgtRope ) )
}

summaryKruschke=function(MCMC,compVal=NULL, rope=NULL,credMass=NULL){
  summaryInfo = NULL
  summaryInfo = cbind(summaryInfo, "Estimated"= SUMK(MCMC,
                                                     compVal=compVal,
                                                     ROPE=rope,credMass=credMass))
  return(summaryInfo)
}

summaryKruschke(MCMC = s$Difprop,compVal=0.05,rope=c(0.1,0.2),credMass=0.975)

##                 Estimated
## Mean            0.3498755
## Median          0.3536088
## Mode            0.3650493
## HDIlevel        0.9750000
## LL              0.1256309
## UL              0.5498100
## CompVal         0.0500000
## PcntGtCompVal 100.0000000
## ROPElow         0.1000000
## ROPEhigh        0.2000000
## PcntLtROPE      0.5000000
## PcntInROPE      5.2500000
## PcntGtROPE     94.2500000

#Bayes factor

brms::hypothesis(s,"Difprop>0.05",alpha=0.01)

## Hypothesis Tests for class :
##                      Estimate Est.Error l-99% CI u-99% CI Evid.Ratio Star
## (Difprop)-(0.05) > 0      0.3       0.1     0.07      Inf        Inf    *
## ---
## '*': The expected value under the hypothesis lies outside the 99% CI.

Trong thí dụ trên, mục tiêu diễn giải là phân phối hậu định của khác biệt tỉ lệ cải thiện triệu chứng (theta) giữa phân nhóm điều trị và placebo.

Nhi đặt ngưỡng so sánh CompVal = 0.05, như vậy nó cho phép kiểm tra giả thuyết H0: Difprop < hay = 5% và H1:Difprop > 5%; đồng thời Nhi đặt khoảng ROPE từ 10-20%, như vậy giả thuyết ở đây là: nếu khác biệt rơi vào trong khoảng 0.1 - 0.2 (ROPE), xem như không có ý nghĩa lâm sàng, và 2 giả thuyết khác là: Dif < 5% và Dif > 5%.

Kết quả cho thấy có 99.875 % mật độ phân phối hậu định của khác biệt nằm bên phải CompVal (>5%), chỉ có 1.125 % mật độ phân phối hậu định khác biệt tỉ lệ nằm ngoài và dưới ngưỡng bên trái=10% của ROPE, 6.375 % phân phối hậu định rơi vào bên trong ROPE,nhưng có đến 92.5% mật độ nằm bên ngoài và cao hơn ngưỡng trên = 20% của ROPE.

Nếu dùng BF với ngưỡng H1 là 5% khác biệt, ta sẽ có giá trị BF=799, rất lớn; cho thấy độ xác tín rất cao của khác biệt về tỉ lệ cải thiện triệu chứng giữa 2 phân nhóm.

Khi tăng dần ngưỡng H1 từ 5% lên 100%, BF giảm dần như sau:

threshold=rep(NA,20)
BayesFactor=rep(NA,20)

for(i in (1:20)){
  thr=i/20
  threshold[i]=thr
  hyp=paste("(`theta[2]`-`theta[1]`)>",thr,sep="")
  bf=brms::hypothesis(fit,hyp,alpha=0.05)
  BayesFactor[i]=bf$hypothesis$Evid.Ratio
}

bfdf=cbind(threshold,BayesFactor)%>%as_tibble()

bfdf%>%ggplot(aes(x=threshold,y=BayesFactor,fill=BayesFactor))+
  geom_path()+
  geom_point(show.legend = F,size=5,shape=21,col="black")+
  geom_text(aes(label=round(BayesFactor,2)),col="black",show.legend = F,angle = 60,nudge_y=50,nudge_x=0.05,size=4)+
  theme_bw()+
  geom_hline(yintercept = c(30,100),linetype=2,col="blue")+
  scale_fill_gradient(low="gold",high="red")

Kết quả của phân tích Bayes có thể tóm tắt như sau: Tỉ lệ bệnh nhân có cải thiện triệu chứng được so sánh giữa 2 phân nhóm Dùng thuốc và Placebo thông qua một mô hình phân phối Binomial theo Bayes. Phân phối hậu định của khác biệt cho thấy nhóm điều trị có tỉ lệ cải thiện cao hơn nhóm placebo đến 35.5 % (khoảng mật độ lớn nhất từ 17.3% đến 51.9%).

8 Tổng kết

Bài thực hành đến đây là chấm dứt. Ngoài một số kỹ năng viết code cho mô hình trong STAN, Nhi đã chuyển đến cho các bạn một số thông điệp chính như sau:

Cùng một bài toán về bảng chéo và phân tích tần số; chúng ta có thể tiếp cận theo nhiều cách. Mỗi cách tương ứng với một phát biểu câu hỏi/giả thuyết khác nhau và sẽ dẫn ta đến những phương pháp thống kê khác nhau . Tuy Chi2 test là một bài toán rất đơn giản, nhưng khi giải quyết nó bằng BAYES, chúng ta học được rất nhiều điều thú vị, và điều này đúng cho mọi quy trình Bayes. Ít nhất, ta đã có khái niệm về 3 loại phân phối khi viết mô hình Bayes cho Chi2 test, đó là phân phối Chi2,phân phối binomial, và phân phối beta, ta có thể diễn đạt bảng chéo thành likelihood và model.

Chi2 test không phải là cách tiếp cận tối ưu, so với mô hình logistic, log linear và so sánh tỉ lệ bằng mô hình binomial.

Chi2 test không có công dụng so sánh tỉ lệ. Nếu bạn muốn so sánh, nên dùng binomial test.

Hẹn gặp lại các bạn vào 1 dịp khác. Nếu các bạn có hứng thú tham gia reboot cho dự án Bayes for Vietnam, xin liên lạc với nhóm chúng tôi.

Xin cảm ơn và hẹn gặp lại