Kiểm định t theo Bayes

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ứ hai trong dự án Bayes for Vietnam, được reboot lại từ project cùng tên năm 2016. Trong bài trước, Nhi đã giới thiệu một số thông tin cơ bản về trường phái Bayes, một cách tiếp cận khác cho những phân tích và suy diễn thống kê. 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 định.

Bản thân trường phái frequentist và null hypothesis testing (NHST) không có gì sai cả, nhưng nó không đủ khả năng để giải quyết những khủng hoảng và bất cập trong suy diễn thống kê ứng dụng vào nghiên cứu. Khái niệm ý nghĩa thống kê và ngưỡng p_value < 0.05 đã sinh ra nhiều hệ quả tai hại, gồm những cách gian lận rất đáng trách, thí dụ :

• Bịa đặt giả thuyết H0 và H1 SAU KHI nhìn thấy p_value

• Thực hiện hàng loạt kiểm định theo kiểu ăn may (tra tấn dữ liệu), dấu kết quả p>0.05 và khoe kết quả p<0.05

• Hack trị số p (và cả CI) để có ý nghĩa thống kê bằng mọi giá

Như vậy NHST không giải quyết mọi vấn đề :thí dụ ý nghĩa thống kê hoàn toàn khác với ý nghĩa lâm sàng, mà nó còn gây ra vấn đế mới. Nó cũng gây ngộ nhận vì không có nhiều người hiểu bản chất p-value là gì.

Bayes có 2 ưu điểm mà không thể có trong trường phái cũ (frequentist), đó là:

1. Hòa hợp giữa giả định/kết quả có trước và dữ liệu thực tế,

2. Rất linh hoạt: Bayes không còn bị ràng buộc bởi các giả định điều kiện như phân phối chuẩn, outliers, cỡ mẫu, phương sai đồng nhất…; và khi có phân phối hậu định trong tay, có nhiều cách để diễn giải, kiểm định giả thuyết nghiên cứu và mỗi cách đều khắc phục các bất cập của trị số p.

2 Ngôn ngữ STAN

Có nhiều cách để thực hành Bayes. Đơn giản nhất là dùng những chương trình (package) viết sẵn, thí dụ BayesFactor, MCMC hoặc các package chuyên biệt cho từng loại model. Các công cụ này đưa bạn thẳng tới kết quả cuối cùng là Bayes factor của tham số, nhưng không cho biết cơ chế hoạt độngbên trong và bạn cũng không tùy chỉnh được cơ chế này.

Phức tạp hơn, chúng ta có những ngôn ngữ lập trình cho Gibs sampler như BUGs và JAGs, chúng đều có giao thức trong R và các package ứng dụng. Cuối cùng, chúng ta có ngôn ngữ STAN với cơ chế hoàn toàn khác, và STAN cũng có giao thức và package ứng dụng trong R như rstan, rstanarm, brms…

Năm 2013, Nhi bắt đầu học Bayes, vừa đọc lý thuyết, vừa viết code bằng WinBUGs, sau đó là JAGs. Cuốn sách của J . Kruschke và Richard McElreath đã cải đạo cho Nhi từ tín đồ của phái Frequentist thành tín đồ của phái Bayes. Sau đó STAN ra đời, cùng với 2nd edition của cuốn sách của Kruschke và hàng loạt các giao thức STAN trong R đã cho phép Nhi thực hành Bayes một cách thoải mái hơn. Do STAN có nhiều ưu điểm và được hỗ trợ đầy đủ nên dự án Bayes for Vietnam chọn STAN làm ngôn ngữ cho tutorial.

STAN là một ngôn ngữ lập trình thống kê xác suất, được phát triển cho R, Python và Matlab.Cơ chế tạo ra chuỗi Markov Monte Carlo của STAN hoàn toàn khác biệt với JAGs và BUGs là 2 sampler có trước. STAN sử dụng phương pháp Hamiltonian Monte Carlo (HMC) sampling dựa vào lý thuyết của Metropolis-Hastings (1953, 1970). Bạn có thể install STAN từ đây:

https://github.com/stan-dev

3 T test

Kể từ lúc William Sealy Gosset (1908) phát minh ra kiểm định t và phân phối t dưới bút danh là Student, t test đã trở thành công cụ được sử dụng phổ biến nhất trong hầu hết các nghiên cứu y học, bao gồm luận văn của các sinh viên. Mục đích giản dị của t test là để so sánh giá trị trung bình giữa 2 phân nhóm dựa vào null hypothesis testing trên trị số t. Thực ra test t chính là một mô hình tuyến tính mô tả biến kết quả Y theo 2 biến số giả chỉ phân nhóm. Trong R, ta có thể thực hiện test t dựa vào công cụ trực tiếp hoặc qua mô hình glm.

Mục tiêu của bài thực hành này là giới thiệu một cách tiếp cận hoàn toàn mới theo trường phái Bayes để so sánh 2 phân nhóm, thay thế cho test t.

Do bài toán này rất phổ biến, đã có nhiều công cụ đã được tạo ra trong R cho phép làm trực tiếp t-test Bayes, thí dụ package BayesFactor, BEST… Tuy nhiên, Nhi sẽ dẫn các bạn đi theo con đường dài hơn, làm thủ công mọi thứ từ code mô hình cho đến diễn giải kết quả, nhưng đổi lại, các bạn sẽ hiểu rõ hơn về cơ chế của mô hình, cấu trúc STAN code, về prior, và kiểm soát mọi thứ, kể cả vẽ biểu đồ theo ý thích.

Nhi cũng sẽ giới thiệu với các bạn cả 3 cách suy diễn Bayes, dùng Bayes Factor, ROPE và ngưỡng so sánh. Khi kết thúc bài thực hành, các bạn sẽ có trong tay một quy trình hoàn chỉnh thay thế hoàn toàn t test cổ điển bằng t test Bayes.

4 Từ t test đến t test 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).

Nếu các bạn tiếp cận được với những tài liệu Thống kê độc đáo, như giáo trình SPSS và R của Andy Field hay quyển Statistical Rethinking của Richard McElreath, bạn sẽ nhận ra rằng quy trình so sánh một biến định lượng giữa các phân nhóm thực chất là một mô hình hồi quy tuyến tính. Suy rộng ra, gần như tất cả những câu hỏi nghiên cứu đều có thể diễn đạt bằng ngôn ngữ mô hình. Do đó, học về mô hình là một cách đầu tư khôn ngoan vì nó trang bị cho các bạn một công cụ phổ quát để tiếp cận vấn đề.

Trong bài toán so sánh đại lượng Y giữa các phân nhóm, mô hình của chúng ta sẽ gồm những yếu tố sau đây :

Một đại lượng Y là biến liên tục mà ta muốn khảo sát, ta gọi nó là biến kết quả (outcome). Giá trị outcome Y này được mô tả bằng một kiểu phân phối xác định.

Theo tinh thần của t test của Gosset, outcome là giá trị trung bình của biến Y, và khác biệt trung bình giữa 2 phân nhóm G1 và G2. Khác biệt trung bình này được mô tả bằng 1 phân phối t. Như ta biết, phân phối t có 3 tham số là mu (location), sigma (scale) và kiểu hình (hay độ tự do) Nu.

Mục tiêu của t test Bayes là tìm phân phối hậu định của 3 tham số này.Để làm việc đó, ta phải dựng mô hình trong STAN

5 Dữ liệu minh họa

Bài toán minh họa của chúng ta là dataset ICU, đây là một nghiên cứu mô tả cắt ngang, khảo sát huyết áp và mạch của 200 bệnh nhân shock nhiễm khuẩn trong khoa ICU. Câu hỏi nghiên cứu giả định của chúng ta là : So sánh huyết áp tâm thu (SysBP) giữa 2 phân nhóm bệnh nhân Tử vong và Sống sót.

# The Data and Story Library (DASL), http://lib.stat.cmu.edu/DASL/Datafiles/ICU.html

df=read.csv("https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/Stat2Data/ICU.csv")%>%as_tibble()

df$Outcome<-recode_factor(df$Survive,
                         `0`="Dead",
                         `1`="Survived")

df$Survive<-recode_factor(df$Survive,
                          `0`="1",
                          `1`="2")

psych::describeBy(df$SysBP,df$Survive)

## 
##  Descriptive statistics by group 
## group: 1
##    vars  n   mean    sd median trimmed   mad min max range skew kurtosis
## X1    1 40 118.83 41.08    126  117.22 32.62  36 256   220  0.6      1.4
##     se
## X1 6.5
## -------------------------------------------------------- 
## group: 2
##    vars   n   mean   sd median trimmed   mad min max range skew kurtosis
## X1    1 160 135.64 29.8    132  133.97 29.65  48 224   176 0.42     0.34
##      se
## X1 2.36

df%>%ggplot(aes(x=SysBP,fill=Outcome))+
  geom_density(alpha=.6,col="black")+
  geom_histogram(aes(y=..density..,col=Outcome),alpha=0.5)+
  facet_wrap(~Outcome,ncol=1,scales="free_y")+
  theme_bw()

6 t-test theo trường phái frequentist

Đầu tiên, Nhi sẽ làm t-test theo trường phái cổ điển

t.test(df$SysBP~df$Outcome,alternative="less")

## 
##  Welch Two Sample t-test
## 
## data:  df$SysBP by df$Outcome
## t = -2.4341, df = 49.726, p-value = 0.009281
## alternative hypothesis: true difference in means is less than 0
## 95 percent confidence interval:
##       -Inf -5.237818
## sample estimates:
##     mean in group Dead mean in group Survived 
##               118.8250               135.6438

effectsize=function(outcome,factor,data){
  d=data
  n2=Hmisc::describe(d[factor])%>%.[1]%>%.[[factor]]%>%.$values%>%.$frequency%>%.[1]
  n1=Hmisc::describe(d[factor])%>%.[1]%>%.[[factor]]%>%.$values%>%.$frequency%>%.[2]
  var=d[outcome]%>%split(d[factor])%>%map(~var(.))
  v2=var[[1]]%>%.[1]
  v1=var[[2]]%>%.[1]
  m2=d[outcome]%>%split(d[factor])%>%map(~psych::describe(.))%>%.[1]%>%.[[1]]%>%.$mean
  m1=d[outcome]%>%split(d[factor])%>%map(~psych::describe(.))%>%.[2]%>%.[[1]]%>%.$mean
  s=sqrt(((n1-1)*v1+(n2-1)*v2)/(n1+n2-2))
  cohend=((m1-m2)/s)
  hedgeG=((m1-m2)/s)*(1- (3/(4*(n1+n2)-9)))
  cat("Effect size for Independent t test","\n")
  cat("Cohen's d value=",cohend,"\n")
  cat("Hedge's G=",hedgeG,"\n")
}

effectsize(outcome="SysBP",factor="Outcome",data=df)

## Effect size for Independent t test 
## Cohen's d value= 0.5201267 
## Hedge's G= 0.5181541

Kết quả cho thấy huyết áp giảm một cách có ý nghĩa thống kê ở nhóm bệnh nhân tử vong so với nhóm sống sót (t=-2.43, p=0.0093, Cohen’s d=0.52).

7 Viết code STAN cho t-test BAYES

Trước khi tiếp cận dữ liệu (data), ta đã có thể đặt giả thuyết về mô hình. Trạng thái sơ khai này của mô hình được xác định bằng các priors (phân phối tiền nghiệm, hay tiền định) của Outcome và tất cả các Tham số trong mô hình

library(rstan)

stanmodelcode ="

data {
int<lower=1> N;                               // Co mau (lower=1: de kiem tra null data ?)
int<lower=2> Ng;                              // So phan nhom (lower =2: it nhat phai co 2 phan nhom)
real<lower=5> n1;                              // Co mau phan nhom 1: it nhat phai co 5 cases
real<lower=5> n2;                              // Co mau phan nhom 2: it nhat phai co 5 cases
vector[N] y;                                  // Bien ket qua Y la 1 vector co do dai=N
int<lower=1, upper=2> groupID[N];            //  group ID la 1 so nguyen nhan 2 gia tri =1 hoac =2
}

transformed data{
real meany;                                   // Bien meany la trung binh cua Y la 1 so thuc; prior cua no la mu
meany = mean(y);                              // Xac dinh gia tri meany= ham mean
}

parameters {
vector[2] mu;                                 // Khai bao tham so cua phan phoi Student-t: mu, sigma va nu; mu la 1 vector 2 gia tri
vector<lower=0>[2] sigma;                     // sigma la 1 vector, gioi han duoi =0, nhan 2 gia tri
real<lower=0, upper=100> nu;                  // nu la do tu do hay kieu hinh cua phan phoi t, gioi han 0:100
}

transformed parameters {                        // Phan hoan chuyen tham so cua mo hinh stan: khong can thiet
}

model {                                         // Khai bao noi dung mo hinh 
// khai bao tien dinh (prior)
mu ~ normal(meany, 10);                       // Tien dinh (prior) cua tham so mu la phan phoi Gaussian, tb=meany, sd=10
sigma ~ cauchy(0, 5);                         // Tien dinh cua tham so Sigma la pp Cauchy(0,5)
nu ~ exponential(1.0/29);                     // Tien dinh cua tham so Nu la pp Exponential voi rate=29


// khai bao likelihood = 1 vong lap

for (n in 1:N){
y[n] ~ student_t(nu, mu[groupID[n]], sigma[groupID[n]]);
}

}

generated quantities {                        // Tinh toan cac tham so trong t_test 
vector[N] Res;                               // Phan phoi hau dinh cua phan vi gia tri du bao Y theo mo hinh, yRep la 1 vector do dai N
real muDiff;                                  // Khac biet trung binh : muDiff la 1 so thuc
real s;                                       // pooled sd
real Cohensd;                                 // Cohen's d effect size (Jacob Cohen) la 1 so thuc
real HedgesG;                                 // Hedges's G effect size ( Larry Hedges, 1981), la so thuc

for (n in 1:N){
  Res[n] = student_t_rng(nu, mu[groupID[n]], sigma[groupID[n]]);
}

muDiff = mu[2] - mu[1];
s=sqrt(((n1-1)*sigma[1]^2+(n2-1)*sigma[2]^2)/(n1+n2-2));      //Tinh pooled sd
Cohensd = muDiff / s;                                      // Tinh Cohen's d
HedgesG = (muDiff / s)*(1- (3/(4*(n1+n2)-9)));                                    // Tinh Hedges's G
}
"

Một chương trình (program) Bayes trong ngôn ngữ STAN có cấu trúc nhiều phần, với danh pháp như sau:

1. Phần data: có nội dung là khai báo các tham số về dữ liệu, mỗi tham số phải được định danh, và định dạng. Trong thí dụ này, có 6 tham số data là: cỡ mẫu N (là 1 số nguyên >1), số phân nhóm Ng (số nguyên dương > 2), cỡ mẫu mỗi phân nhóm (n1,n2) – lẽ ra cũng là số nguyên, nhưng được định dạng là số thực (real), vì Nhi muốn dùng chúng để tính effect size (Cohen’s d), và trong STAN các phép tính số nguyên sẽ cho ra kết quả là số nguyên, làm tròn, là điều ta không mong muốn. Biến kết quả y là 1 vector kích thước N, Tên phân nhóm groupID là 1 số nguyên, nhận 2 giá trị 1, 2.

2. Phần transformed data để thực hiện hoán chuyển dữ liệu (thường là biến kết quả) trước khi đưa vào hàm likelihood. Việc hoán chuyển này có bản chất là áp dụng 1 link function lên Y. Trong thí dụ này, outcome của mô hình t-test là 1 giá trị trung bình của Y, nên ta khai báo biến outcome sau hoán chuyển là meany (1 số thực), và công thức hoán chuyển : mean y = mean(y).

3. Phần tham số (parameters): Đây là những tham số mà ta cần tìm phân phối hậu định, kết quả của mô hình Bayes. Trong thí dụ này, ta có 3 tham số, mu, sigma và nu tương ứng với location, scale và shape của một phân phối Student t. mu là 1 vector, sigma là 1 vector, có 2 nhánh cho 2 phân nhóm (như vậy kết quả sẽ có mu1, mu2, sigma1, sigma2); nu là 1 số thực giới hạn từ 0-100.

4. Phần tham số hoán chuyển (transformed parameters): không dùng đến ở đây nên được để trống

5. Phần model : nội dung của mô hình Bayes: Đầu tiên, ta khai báo Priors cho mỗi tham số. Nhắc lại : Prior là một hàm mật độ xác suất mà ta giả định (tin) rằng hợp lý nhất cho tham số. Prior là thông tin có trước khi nhìn thấy kết quả, nó có thể là kết quả trong quá khứ (y văn), kinh nghiệm cá nhân, hoặc giả định lý thuyết. Ở đây, ta lần lượt đặt 3 prior cho 3 tham số mu, sigma và nu như sau:

Mu là location của phân phối t, của biến kết quả Mean(Y). Vì mu là 1 số thực, và sự biến thiên ngẫu nhiên của nó có bản chất là kết quả của một quá trình thêm/bớt những cá thể có giá trị y cao/thấp một cách ngẫu nhiên, do đó phân phối hợp lý nhất của mu sẽ là Gaussian (normal). Thật vậy, Gaussian là phân phối hợp lý nhất cho 1 biến liên tục khi mà ta không có khả năng đưa ra bất cứ giả định nào khác, nó là quy luật của mọi hiện tượng sinh lý, tự nhiên. Phân phối Gaussian (normal) này có location= trung bình của y và scale=10.

Prior của sigma và nu khó giải thích hơn (Nhi không có giả định nào nên buộc phải dùng phân phối cauchy cho sigma và exponential với rate =29 cho nu, 2 prior này đã được kiểm chứng trong package BEST và tài liệu của Kruschke). Để dễ hình dung, ta có thể thử vẽ density plot của 2 prior này. Ghi chú là bạn có thể dùng phân phối Gamma thay cho exponential cũng được.

p1=rexp(5000,1/29)%>%
  as.data.frame()%>%
  ggplot(aes(x=.))+
  geom_density(alpha=0.5,col="black",fill="red")+
  theme_bw()+ggtitle("Prior cho nu")

p2=rnorm(5000,mean(df$SysBP),10)%>%
  as.data.frame()%>%
  ggplot(aes(x=.))+
  geom_density(alpha=0.5,col="black",fill="red")+
  theme_bw()+ggtitle("Prior cho Mu")

p3=rcauchy(5000,location = 0,scale =5)%>%
  as.data.frame()%>%
  ggplot(aes(x=.))+xlim(c(0,max(df$SysBP)))+
  geom_density(alpha=0.5,col="black",fill="red")+
  theme_bw()+ggtitle("Prior cho Sigma")

p4=df%>%ggplot(aes(x=SysBP))+
  geom_density(alpha=0.5,col="black",fill="red")+
  theme_bw()+ggtitle("Phân phối thực tế của Y")

gridExtra::grid.arrange(p1,p2,p3,p4,ncol=2)

## Warning: Removed 2486 rows containing non-finite values (stat_density).

Nếu không có prior, STAN sẽ dùng flat prior, nhưng kết quả sẽ không chính xác

Sau khi khai báo prior, ta viết 1 vòng lặp cho hàm likelihood : nội dung của nó là ước lượng giá trị của outcome y với quy luật phân phối student-t được xác định bằng 3 tham số mu, sigma, nu, tùy theo phân nhóm groupID.

Đến đây là hết nội dung model.

6. Phần cuối cùng của chương trình STAN có tên là generated quantities, đây là một chức năng rất hữu dụng, vì nó cho phép chúng ta tạo ra bất kì trị số thống kê mà ta thích dựa vào mô hình. Đây chính là giao điểm giữa trường phái Bayes và thống kê cổ điển, ta có thể rút ra phân phối hậu định của những trị số thống kê cho t-test mà ta từng biết, như :

Phân phối hậu định của giá trị Y cho mỗi cá thể trong mẫu : Res

Phân phối hậu định cho khác biệt trung bình : muDiff = mu phân nhóm 2 – mu phân nhóm 1

Phân phối hậu định cho Effect -size của t test : Cohen’s d và Hedge’s G

(nếu thích, bạn có thể tính cả p_value cho 1 null hypothesis testing: thí dụ H1=muDiff > 0, chẳng hạn, nhưng Nhi không muốn đi lạc quá xa về phái frequentist ở đây).

8 Thi hành mô hình STAN

Tác giả John Kruschke đã chia sẻ kinh nghiệm của ông khi viết code cho mô hình STAN : đi ngược từ likelihood model lên trên, và di chuyển con trỏ chuột giữa các block, để đảm bảo không sót, không sai (chỉ cần thiếu một dấu ; hay {, chương trình sẽ không thể compile được và báo lỗi, dĩ nhiên đây là điều tốt, vì quy trình sampling MCMC rất tốn kém thời gian, nên STAN kiểm tra lỗi cú pháp trước khi thi hành.

Sau khi đã đóng gói mô hình, ta khai báo data cho mô hình trong 1 list, đây chính là các tham số trong block data. Thí dụ ở đây : N=200, Ng=2, n1=40, n2=160, groupID là một vector số nguyên từ column Survive trong dataframe, y là vector của column SysBP trong data frame.

Bây giờ ta đưa đoạn code và data vào hàm stan để bắt đầu thi hành mô hình. Trước tiên, STAN sẽ kiểm tra cú pháp của đoạn code, tương hợp giữa các tham số trong đoạn code và list data. Nếu tất cả đều hợp chuẩn, đoạn code và data sẽ được compile thành 1 program C++, và thi hành. Program này kích hoạt một sampler để bắt đầu lấy mẫu các chuỗi MCMC theo cơ chế mà ta đã đưa ra, thí dụ : 10.000 lượt iteration, khởi động 1000 lượt, chạy song song trên 4 cores, và rút gọn MCMC (thin =10).

groupID=as.numeric(df$Survive)

data= list(N=200, Ng=2, n1=40, n2=160, groupID=groupID, y=df$SysBP)

set.seed(1234)

fit = stan(model_code=stanmodelcode, data=data, iter=10000, warmup=1000, 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 file1f5c17726c4c.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"

Mô hình chỉ mất vài phút để converge, bạn sẽ có 1 object khoảng 6.5 Mb, bên trong nó là các chuỗi MCMC cho tất cả những tham số, trị số mà bạn đã khai báo. Cụ thể, nó chứa 200 chuỗi MCMC cho giá trị dự báo SysBP của từng cá thể (biến Res), 2 chuỗi MCMC cho mu1, mu2, 2 chuỗi cho sigma1, sigma2, 1 chuỗi cho muDiff (kết quả chính của t test), 1 chuỗi cho Cohen’s d và 1 chuỗi cho Hedge’s G

9 Tham dò kết quả phân phối hậu nghiệm

Đầu tiên, chúng ta xem qua phân phối hậu định của mu, sigma, muDiff (kết quả chính), nu và effect size Cohen’s d, Hedges’ G

print(fit, digits=3, pars=c('mu', 'sigma', 'muDiff', 'Cohensd', 'HedgesG','nu'))

## Inference for Stan model: 999c243cf5747c34c22dbcd0b41da6a4.
## 4 chains, each with iter=10000; warmup=1000; thin=10; 
## post-warmup draws per chain=900, total post-warmup draws=3600.
## 
##             mean se_mean     sd    2.5%     25%     50%     75%   97.5%
## mu[1]    122.336   0.091  5.443 111.888 118.691 122.137 125.941 133.372
## mu[2]    134.670   0.039  2.292 130.225 133.129 134.613 136.191 139.233
## sigma[1]  37.172   0.088  5.233  27.840  33.537  36.974  40.365  48.359
## sigma[2]  27.477   0.035  2.123  23.214  26.096  27.556  28.908  31.607
## muDiff    12.333   0.098  5.901   0.391   8.419  12.435  16.325  23.448
## Cohensd    0.418   0.003  0.203   0.013   0.282   0.419   0.555   0.809
## HedgesG    0.417   0.003  0.202   0.013   0.281   0.417   0.553   0.806
## nu        20.898   0.295 17.375   4.576   9.101  14.625  26.419  71.901
##          n_eff  Rhat
## mu[1]     3600 0.999
## mu[2]     3530 1.002
## sigma[1]  3565 1.000
## sigma[2]  3600 1.001
## muDiff    3600 1.000
## Cohensd   3600 1.000
## HedgesG   3600 1.000
## nu        3472 1.000
## 
## Samples were drawn using NUTS(diag_e) at Sun Sep 17 12:34:30 2017.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at 
## convergence, Rhat=1).

Chi tiết quan trọng nhất, đó là phân phối hậu định của muDiff (khác biệt trung bình) không chứa giá trị zero, và trung vị = 12.435, do đó ta có rất nhiều khả năng giả thuyết: SysBP thấp hơn ở nhóm tử vong là đúng. Hai giá trị effect-size theo Cohen và Hedges thấp hơn 1 chút so với t-test cổ điển, nhưng vẫn có thể diễn giải là “hiệu ứng mạnh”.

Ta có thể vẽ phân phối hậu nghiệm của MuDiff

muDiff = extract(fit, pars='muDiff',inc_warmup=F)$muDiff%>%as_tibble()

muDiff$group=factor(rep(c(1:100), e=3600/100))

muDiff%>%ggplot(aes(x=value))+ 
  geom_density(alpha=0.5,col="red4",size=1) +
  geom_point(y=0, alpha=.01, size=2,col="red3") +
  geom_line(aes(group=group, color=group), stat='density', alpha=.2,show.legend = F) +scale_x_continuous("Khác biệt trung bình")+
  theme_bw()

Tiếp theo, Nhi sẽ khảo sát phân phối hậu định của 200 giá trị dự báo của SysBP từ mô hình Bayes, ta muốn biết là prior có hợp lý không ? phân phối hậu định của kết quả có phù hợp với quan sát thực tế hay không ? (vì phân phối hậu định dựa vào prior). Kết quả cho thấy phân phối hậu định của Y (2 hình chuông màu xanh, đỏ) khá phù hợp với quan sát thực tế (2 đường cong màu đen)

resdf=extract(fit, pars='Res',inc_warmup=F)$Res%>%reshape2::melt(.)%>%as_tibble()
colnames(resdf) = c('iteration', 'case', 'SysBP' )
resdf$groupID=factor(rep(groupID, e=3600))

resdf$groupID=recode_factor(resdf$groupID,
                          `1`="Dead",
                          `2`="Survived")

resdf$case=as.factor(resdf$case)

resdf%>%ggplot(aes(x=SysBP)) + 
  geom_density(aes(group=groupID, fill=groupID), color=NA, alpha=.25) +
  geom_line(aes(group=case, col=groupID), stat='density', alpha=.05) +
  geom_point(aes(x=df$SysBP, y=0, color=factor(Outcome)), alpha=.15, size=2, data=df)+
  xlim(c(min(df$SysBP),max(df$SysBP)))+
  geom_density(aes(group=Outcome, linetype=Outcome, x=df$SysBP), alpha=.05, data=df)+
  theme_bw()

## Warning: Removed 4768 rows containing non-finite values (stat_density).

## Warning: Removed 4768 rows containing non-finite values (stat_density).

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) :

p1=muDiff%>%mutate(.,Iteration=as.numeric(rownames(.)))%>%ggplot(aes(x=Iteration,y=value))+
  geom_path(col="orangered",alpha=0.7)+
  theme_bw()+ggtitle("Khác biệt trung bình")

p2=muDiff%>%ggplot(aes(x=value))+ 
  geom_density(alpha=0.5,fill="orangered",size=1) +
  scale_x_continuous("Khác biệt trung bình")+
  theme_bw()

cohend = extract(fit, pars='Cohensd',inc_warmup=F)$Cohensd%>%as_tibble()

p3=cohend%>%mutate(.,Iteration=as.numeric(rownames(.)))%>%
  ggplot(aes(x=Iteration,y=value))+
  geom_path(col="red",alpha=0.7)+
  theme_bw()+ggtitle("Effect-size")

p4=cohend%>%ggplot(aes(x=value))+ 
  geom_density(alpha=0.5,fill="red",size=1) +
  scale_x_continuous("Cohen's d")+
  theme_bw()

nu = extract(fit, pars='nu',inc_warmup=F)$nu%>%as_tibble()

p5=nu%>%mutate(.,Iteration=as.numeric(rownames(.)))%>%
  ggplot(aes(x=Iteration,y=value))+
  geom_path(col="purple",alpha=0.7)+
  theme_bw()+ggtitle("Nu")

p6=nu%>%ggplot(aes(x=value))+ 
  geom_density(alpha=0.5,fill="purple",size=1) +
  scale_x_continuous("Cohen's d")+
  theme_bw()

gridExtra::grid.arrange(p1,p2,p3,p4,p5,p6,ncol=2)

10 Suy diễn Bayes

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) :

11 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 (thí dụ H1 : muDiff > 0, hoặc Cohen’s d > 0.4 chẳng hạn). 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.

Thí dụ: Ta muốn tính Bayes Factor cho giả thuyết H1: muDiff >3 so với giả thuyết H0: muDiff < hoặc = 3:

Kết quả cho ra BF = 15.29, tỉ trọng chứng cứ yếu

brms::hypothesis(fit,"muDiff>3",alpha=0.05)

## Hypothesis Tests for class :
##                  Estimate Est.Error l-95% CI u-95% CI Evid.Ratio Star
## (muDiff)-(3) > 0     9.33       5.9    -0.75      Inf      15.29     
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

Bây giờ, ta muốn tính Bayes Factor cho giả thuyết H1: muDiff >0 so với giả thuyết H0: muDiff < hoặc =0:

Kết quả cho ra BF = 45.75, tỉ trọng chứng cứ cao hơn 30 nên mức độ xác tín cũng chắc chắn hơn

brms::hypothesis(fit,"muDiff>0",alpha=0.05)

## Hypothesis Tests for class :
##              Estimate Est.Error l-95% CI u-95% CI Evid.Ratio Star
## (muDiff) > 0    12.33       5.9     2.25      Inf      45.75    *
## ---
## '*': The expected value under the hypothesis lies outside the 95% CI.

Bạn có thể nêu ra cùng lúc hàng loạt cặp giả thuyết H0/H1, và biểu diễn cùng lúc nhiều giá trị BF tương ứng cho mỗi giả thuyết, như trong hình dưới đây :

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

for(i in (1:21)){
  thr=i-1
  threshold[i]=thr
  hyp=paste("muDiff>",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_point(show.legend = F,size=5,shape=21,col="black")+
  geom_text(aes(label=round(BayesFactor,2)),col="black",show.legend = F,nudge_y = 5,nudge_x=1,angle = 45)+
  theme_bw()+
  geom_hline(yintercept = c(10,30),linetype=2,col="blue")+
  scale_fill_gradient(low="gold",high="red")

12 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 (thí dụ: huyết áp giảm 0-3 mmHg được xem như không có ý nghĩa lâm sàng). 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.

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 = muDiff$value,compVal=0,rope=c(1,5),credMass=0.95)

##                Estimated
## Mean          12.3331004
## Median        12.4352939
## Mode          13.4199996
## HDIlevel       0.9500000
## LL             0.3912786
## UL            23.4514280
## CompVal        0.0000000
## PcntGtCompVal 97.8611111
## ROPElow        1.0000000
## ROPEhigh       5.0000000
## PcntLtROPE     3.2500000
## PcntInROPE     7.0277778
## PcntGtROPE    89.7222222

Trong thí dụ trên, Nhi đặt ngưỡng so sánh CompVal = 0.0, như vậy nó cho phép kiểm tra giả thuyết H0: muDiff < hay =0 và H1: muDiff >0; đồng thời Nhi đặt khoảng ROPE từ 0-5, như vậy giả thuyết ở đây là: nếu MuDiff rơi vào trong khoảng 0-5 (ROPE), xem như không có ý nghĩa lâm sàng, và 2 giả thuyết khác là: MuDiff <0 và MuDiff >5.

Kết quả cho thấy có 97.86 % mật độ phân phối hậu định của MuDiff nằm bên phải CompVal (>0), chỉ có 3.25 % mật độ phân phối hậu định MuDiff nằm ngoài và dưới ngưỡng 0 của ROPE, 7.03 % phân phối hậu định rơi vào bên trong ROPE,nhưng có đến 89.72 % mật độ nằm bên ngoài và cao hơn ngưỡng trên = 5 của ROPE.

13 Tổng kết

Bài thực hành đến đây là chấm dứt. Qua bài này, các bạn đã nắm được một quy trình hoàn chỉnh để thực hiện t-test BAYES, thay thế cho t-test cổ điển, bao gồm:

• Cấu trúc chương trình và cú pháp của STAN code
• Diễn đạt t-test theo ngôn ngữ mô hình
• Cách chọn Prior
• Khai thác phân phối hậu định
• Suy diễn Bayes bằng 3 phương pháp

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