BAYES-ANOVA cho thí nghiệm lặp lại

Sử dụng ngôn ngữ STAN

Ghi chú:

Bs. Lê Ngọc Khả Nhi là thành viên trong Coreteam của nhóm BAV (Bayesian analysis for Vietnam). Trong bài viết có sử dụng 1 số ý tưởng, R function và code của anh Nhật Nam là team leader. Cảm ơn bạn Đinh Tiến Tài trong việc kiểm tra STAN code.

1 Giới thiệu

Thân chào các bạn, đây là bài thực hành thứ 6 trong dự án Bayes for Vietnam. Mục tiêu của chúng tôi là phổ cập về phương pháp thống kê theo trường phái Bayes nhằm thay thế hoàn toàn những công cụ truyền thống. Nối tiếp bài trước với chủ đề ANOVA đơn biến, bài này sẽ giới thiệu với các bạn thiết kế Thí nghiệm lặp lại (Repeated measure) và cách phân tích dữ liệu loại này. Nhi giả định là các bạn đã quen thuộc với cơ chế của suy diễn Bayes, cấu trúc ngôn ngữ STAN và quy trình khai thác phân phối hậu định, nên một số phần sẽ được giản lược.

Trong bài này chúng ta sẽ thay thế phân tích phương sai cho thí nghiệm lặp lại (Repeated measure ANOVA) bằng mô hình hỗn hợp và hồi quy Bayes. Như thường lệ, bài giảng sẽ đi theo 3 bước như sau:

Trước hết, chúng tôi sẽ ôn lại lý thuyết về ANOVA repeated measure, đưa ra 1 bài toán và minh họa quy trình ANOVA cổ điển theo phái frequentist.

Sau đó chúng tôi sẽ chuyển bài toán này thành mô hình Bayes, và hướng dẫn các bạn viết STAN code cho mô hình.

Cuối cùng, chúng ta sẽ khai thác phân phối hậu định cho các tham số cần quan tâm, và suy diễn Bayes.

2 Ôn tập về ANOVA thí nghiệm lặp lại

Bài toán minh họa

Dữ liệu minh họa trong bài này là một thí nghiệm sinh lý trên động vật vào năm 1994 của GS. J. Ludbrook. Trong thí nghiệm nguyên thủy, mục tiêu nhằm khảo sát hiệu ứng tương tác của 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ỏ. Để đơn giản vấn đề, chúng tôi loại bỏ yếu tố can thiệp với đối vận serotonin, và chỉ giữ lại nhóm chứng và hiệu ứng của phenylbiguanide.

Câu hỏi nghiên cứu giả định trong bài toán hiện tại, do đó sẽ là: khảo sát hiệu ứng của liều thuốc phenylbiguanide lên huyết áp của thỏ.

Quy trình thí nghiệm là như sau: 5 con thỏ sẽ được tiêm tĩnh mạch lần lượt 6 liều phenylbiguanide từ D1 (thấp nhất) đến D6 (cao nhất) với khoảng cách 10 phút. Huyết áp thỏ được đo cùng lúc (như vậy mỗi con thỏ được đo huyết áp 6 lần).

library(tidyverse)

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

# Within only

df$step=rep(c("D1","D2","D3","D4","D5","D6"),5)

df=subset(df,Treatment=="Control")%>%.[,c(2:3,6,7)]%>%as_tibble()%>%mutate(logDose=log(.$Dose))

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

BPchange	Dose	Animal	step	logDose
0.5	6.25	R1	D1	1.832582
4.5	12.50	R1	D2	2.525729
10.0	25.00	R1	D3	3.218876
26.0	50.00	R1	D4	3.912023
37.0	100.00	R1	D5	4.605170
32.0	200.00	R1	D6	5.298317

Bảng thống kê mô tả:

# Descriptive analysis

d=psych::describeBy(df$BPchange,df$step)

t1=rbind(d$D1[,c(2:13)], 
         d$D2[,c(2:13)], 
         d$D3[,c(2:13)], 
         d$D4[,c(2:13)],
         d$D5[,c(2:13)],
         d$D6[,c(2:13)])%>%as.data.frame()%>%round(.,2)

row.names(t1)=c("D1","D2","D3","D4","D5","D6")

knitr::kable(t1)

	n	mean	sd	median	trimmed	mad	min	max	range	skew	kurtosis	se
D1	5	1.00	0.40	1.0	1.00	0.37	0.50	1.5	1.00	0.00	-1.91	0.18
D2	5	2.35	1.39	1.5	2.35	0.37	1.25	4.5	3.25	0.57	-1.70	0.62
D3	5	5.60	2.70	5.0	5.60	1.48	3.00	10.0	7.00	0.64	-1.40	1.21
D4	5	17.40	5.46	16.0	17.40	4.45	12.00	26.0	14.00	0.54	-1.54	2.44
D5	5	27.80	7.19	27.0	27.80	8.90	20.00	37.0	17.00	0.13	-2.05	3.22
D6	5	27.20	6.22	29.0	27.20	5.93	18.00	33.0	15.00	-0.42	-1.81	2.78

Bây giờ chúng ta sẽ phân tích bài toán này.

Kết quả là BPchange, 1 biến định lượng liên tục. Mỗi con thỏ đều tham gia vào cả 6 lần đo tương ứng D1 đến D6. Khảo sát trực quan cho thấy khuynh hướng chung là BPchange tăng dần từ D1 đến D6 với dao động ngẫu nhiên tùy theo cá thể thỏ;

# Exploration

library(viridis)

## Loading required package: viridisLite

df%>%ggplot(aes(x=step,y=BPchange,fill=Animal,group=1))+
  geom_line(aes(color=Animal),size=1)+
  geom_point(shape=21,size=3,col="black")+
  theme_bw()+
  scale_color_viridis(discrete=T,option="C",direction=-1)+
  scale_fill_viridis(discrete=T,option="C",direction=-1)+
  facet_wrap(~Animal,scale="free",ncol=2)

Thí nghiệm lặp lại là một thiết kế nghiên cứu trong đó kết quả định lượng Yij sẽ được khảo sát nhiều lần tại nhiều thời điểm (Tj) khác nhau hoặc dưới những điều kiện/phương pháp (Fj) khác nhau tên cùng đối tượng (Si).

Thiết kế này có rất nhiều ứng dụng trong nghiên cứu y học, thí dụ:

• Đánh giá hiệu quả điều trị của phương pháp phẫu thuật: mỗi bệnh nhân được khảo sát 1 chỉ số lâm sàng Y trước và sau khi phẫu thuật (lưu ý: Repeated measure ANOVA là dạng tổng quát của paired samples Student t-test).

• Khảo sát tính tương hợp của 3 máy đo đường huyết: Cùng một bệnh nhân sẽ được đo đường huyết bởi 3 máy đo, theo thứ tự ngẫu nhiên. Một thiết kế tương tự: khảo sát sự tương hợp chẩn đoán của 3 bác sĩ chẩn đoán hình ảnh khi cho họ đọc kết quả của cùng 1 phim CT-scan.

• Nghiên cứu theo dõi kéo dài (longitudinal study): Một nhóm bệnh nhân được theo dõi diễn tiến của marker lâm sàng Y mỗi tháng một lần trong suốt quá trình điều trị bằng loại thuốc X, trong thời gian 1 năm.

…

Có 2 chi tiết quan trọng cần lưu ý:

1. Cần phân biệt giữa 2 thiết kế: Khảo sát lặp lại (repeated measure) và nghiên cứu dài hạn (longitudinal study):

Repeated measure không bắt buộc có biến số thời gian. Đa số trường hợp, việc lặp lại thí nghiệm đồng nghĩa với việc khảo sát tại nhiều thời điểm khác nhau và theo trình tự thời gian; thí dụ trong trường hợp này, 6 liều thuốc khác nhau sẽ được đưa vào cơ thể mỗi con thỏ theo thứ tự tăng dần, bắt đầu từ D1 và kết thúc ở D6; khi đó yếu tố Thời gian đã hòa nhập làm một với yếu tố Can thiệp và Điều kiện thí nghiệm, và bị ẩn đi.

Trong thí dụ thứ 2, thậm chí người ta còn muốn triệt tiêu vai trò của thời điểm/thời gian, bằng cách ngẫu nhiên hóa thứ tự của thí nghiệm; thí dụ cá thể A được đo bằng 3 phương pháp theo trình tự M1,M2,M3, nhưng cá thể B thì đổi sang trình tự khác: M2,M1,M3. Một số trường hợp khác, thậm chí không có yếu tố thời gian vì các thí nghiệm được thực hiện đồng thời, thí dụ thí nghiệm In-vitro trong đó một đoạn mạch máu được cắt ra thành 3 phần và mỗi phần cho tiếp xúc với 1 hoạt chất co mạch khác nhau.

Trong khi đó, Thời gian là một biến số quan trọng, tuyến tính và biểu kiến trong longitudinal study. Thời gian này có thể rất dài (đúng như tên gọi : longitudinal), qua nhiều năm, tháng; nhưng cũng có thể rất ngắn, thí dụ bạn khảo sát sự biến thiên của 1 đoạn tín hiệu Y với tần số lấy mẫu 10 Hz, kéo dài 30 giây. Thời gian nên được khảo sát như 1 biến liên tục trong longitudinal study.

2. Dạng chính tắc của ANOVA thường xem xét 2 hiệu ứng: between và within effect. Thí dụ thí nghiệm gốc trên thỏ này, between effect do hiệu ứng của đối vận serotonin, còn within effect do liều thuốc phenylbiguanide,

Trong giới hạn của bài này, chúng ta chỉ bàn về Repeated measure trong một trường hợp đơn giản nhất và đặc biệt, đó là cả 3 hiệu ứng Thời gian, Can thiệp, điều kiện khảo sát đã hòa nhập làm một trong biến số Dose step, tức là chỉ còn duy nhất within effect.

Tuy nhiên nguyên tắc của ANOVA là không đổi. Ta vẫn định nghĩa được F từ Mean square (MS) của mô hình, so với MS của Residual

Phương sai của kết quả Y có thể được giải thích qua một mô hình đơn giản như sau :

\[Y_{ij} \sim \mu + S_{i} + F_{j} + \varepsilon _{ij}\]

Mô hình này ước lượng giá trị Y cho cá thể i tại điều kiện thí nghiệm j từ : trung bình của Y trong quần thể là Mu, hiệu ứng ngẫu nhiên Si (random effect, S = subject, phần biến thiên ngẫu nhiên ở từng cá thể) và hiệu ứng chính do yếu tố F (Fixed effect, hay within subject variance, F ở đây là điều kiện thí nghiệm).

Phần sai số tồn lưu là epsilon, mô hình không thể giải thích được. Như vậy, việc đưa thêm random effect Si vào mô hình cho phép giảm thiểu sai số và phân lập được yếu tố ngẫu nhiên cá thể.

Mô hình này còn có tên là random intercept, vì Si chỉ đại diện cho trạng thái cơ bản của mỗi cá thể trước khi chịu can thiệp bởi thí nghiệm F, và ta giả định là mỗi cá thể sẽ đáp ứng như nhau trong thí nghiệm F.

Một mô hình khác phức tạp hơn, random slope còn xét cả sự tương tác ngẫu nhiên giữa Si và Fj. Tuy nhiên, trong trường hợp này ta chọn mô hình Random Intercept vì cỡ mẫu thấp và mỗi cá thể (i) chỉ được khảo sát 1 lần duy nhất tại mỗi điều kiện Fj.

Đồ thị của 5 mô hình tuyến tính có thể được hình dung như sau:

df%>%ggplot(aes(x=logDose,y=BPchange,fill=Animal))+
  geom_point(alpha=0.7,aes(color=step),show.legend = F,size=2)+
  geom_smooth(aes(col=Animal),
              size=1,
              method="lm",se=F,
              show.legend = F)+
  theme_bw()+
  scale_color_viridis(discrete=T,option="C",direction=-1)+
  scale_fill_viridis(discrete=T,option="C",direction=-1)

Như đã nói trong bài trước, các sinh viên Y khoa thường được dạy xử lý vấn đề theo hướng so sánh trung bình giữa các phân nhóm. Cách tư duy này có thể đúng cho hầu hết trường hợp, như thí dụ này quả thực ta đang muốn so sánh BPchange trung bình giữa 6 phân nhóm (bậc) D1:D6;

df%>%ggplot()+
  geom_boxplot(aes(x=step,y=BPchange,fill=step),
               alpha=0.7,show.legend = F)+
  geom_hline(yintercept = 0,linetype=2,col="red4")+
  theme_bw()+
  scale_color_viridis(discrete=T,option="C",direction=-1)+
  scale_fill_viridis(discrete=T,option="C",direction=-1)

Nhưng chúng tôi khuyến khích cách đặt mục tiêu : khảo sát hiệu ứng của yếu tố liều thuốc lên kết quả Y. Cách suy nghĩ này tổng quát, linh hoạt hơn và tự nhiên dẫn dắt bạn đến việc dùng mô hình.

ANOVA không được giảng dạy một cách hệ thống trong chương trình thống kê, nhất là ANOVA cho thí nghiệm lặp lại. Do đó, khi thực hiện nghiên cứu lâm sàng, các bác sĩ không biết đến phương pháp này. Thói quen thường gặp của họ là chia nhỏ bài toán ra thành nhiều cặp so sánh rồi dùng paired sample t-test hay Wilcoxon sign rank test để giải quyết. Việc này hoàn toàn sai.

Bài toán ANOVA cho phép đo lặp lại có thể được giải quyết bằng 2 cách: Thủ công/cổ điển hoặc Mô hình Mixed model:

Trước hết, ta xem xét các giả định

Mọi phân tích ANOVA đều yêu cầu những giả định như sau :

1. Phân phối chuẩn :

Yj tại mỗi điều kiện thí nghiệm (phân nhóm) phải phân phối chuẩn chung quanh giá trị trung bình Muj.

df%>%
  split(.$step)%>%
  map(~fBasics::shapiroTest(.$BPchange))

## $D1
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.9868
##   P VALUE:
##     0.9672 
## 
## Description:
##  Sat Oct 28 21:31:21 2017 by user: Admin
## 
## 
## $D2
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.8261
##   P VALUE:
##     0.13 
## 
## Description:
##  Sat Oct 28 21:31:21 2017 by user: Admin
## 
## 
## $D3
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.9031
##   P VALUE:
##     0.4272 
## 
## Description:
##  Sat Oct 28 21:31:21 2017 by user: Admin
## 
## 
## $D4
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.927
##   P VALUE:
##     0.576 
## 
## Description:
##  Sat Oct 28 21:31:21 2017 by user: Admin
## 
## 
## $D5
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.9415
##   P VALUE:
##     0.6768 
## 
## Description:
##  Sat Oct 28 21:31:21 2017 by user: Admin
## 
## 
## $D6
## 
## Title:
##  Shapiro - Wilk Normality Test
## 
## Test Results:
##   STATISTIC:
##     W: 0.9149
##   P VALUE:
##     0.4974 
## 
## Description:
##  Sat Oct 28 21:31:21 2017 by user: Admin

Giả định này trên thực tế, hầu như không bao giờ được thỏa mãn, và những kiểm định thống kê để kiểm tra nó đều có nhược điểm nên cách xử trí tốt nhất là nhận xét trực quan và quyết định riêng của bạn (tính đồng dạng, outliers, …).

library(ggjoy)

## Warning: package 'ggjoy' was built under R version 3.4.2

## Loading required package: ggridges

## Warning: package 'ggridges' was built under R version 3.4.2

## The ggjoy package has been deprecated. Please switch over to the
## ggridges package, which provides the same functionality. Porting
## guidelines can be found here:
## https://github.com/clauswilke/ggjoy/blob/master/README.md

df%>%ggplot(aes(y=step,x=BPchange,
                fill=step,
                col=step
))+
  geom_joy(alpha=0.7,scale=1,show.legend = F)+
  geom_rug(alpha=0.6,show.legend = F)+
  geom_vline(xintercept = 0,linetype=2,col="red4")+
  coord_flip()+theme_bw()+
  scale_color_viridis(discrete=T,option="C",direction=-1)+
  scale_fill_viridis(discrete=T,option="C",direction=-1)

## Picking joint bandwidth of 2.16

Đồ thị cho thấy: Có sự vi phạm rõ ràng giả định phân phối chuẩn từ D3:D6, Shapiro - Wilk test đã không thể phát hiện. Mặt khác, cỡ mẫu quá thấp (n=5) cho mỗi phân nhóm - lẽ ra với cỡ mẫu này thậm chí ta phải dùng Non parametric test chứ không thể nghĩ đến một mô hình ANOVA. Vì lý do minh họa, ta sẽ thử đi tiếp:

2. Tính độc lập của quan sát : điều này cho phép giả định về sự độc lập và ngẫu nhiên của sai số eij. Những cá thể trong từng phân nhóm phải độc lập với nhau. Thí nghiệm lặp lại là một trường hợp đặc biệt mà giả định này bị vi phạm một cách có hệ thống và chủ ý, vì mỗi cá thể tham gia vào tất cả bậc j của điều kiện thí nghiệm (phân nhóm), ta cũng đã hiệu chỉnh điều này bằng cách phân lập random effect Si và between subject effect Fj nên không phải bận tâm về giả định này.

3. Phương sai đồng nhất : Mỗi phân nhóm (điều kiện thí nghiệm, bậc của F) phải có phương sai như nhau. Với thí nghiệm Repeated measure, ta có giả định quan trọng hơn là « Sphericity », hay tính đồng nhất về phương sai của các giá trị khác biệt ở mỗi cá thể (trong 1 covariance matrix). Thí dụ, nếu ta có K điều kiện thí nghiệm (lần lặp lại), thì sẽ có k(k-1)/2 cặp khác biệt, giả định sphericity là chúng đều có phương sai như nhau.

Giả định Sphericity có thể kiểm tra bằng test Mauchly, hoặc quan sát covariance matrix. Tuy nhiên bước này có thể bỏ qua nếu ta hiệu chỉnh df của within subject factor (F) bằng hệ số epsilon theo phương pháp Greenhouse-Geisser (GG) hoặc Huynh-Feildt (HF). Epsilon = 1 khi không có vi phạm giả định

Tuy nhiên, giả định sphericity là không cần thiết (không thể bị vi phạm), nếu K<3 (chỉ lặp lại thí nghiệm 2 lần). Sphericity cũng không cần thiết nếu trong mô hình không có tương tác giữa F và S (random intercept model, hay mỗi cá thể chỉ đo 1 lần tại mỗi điều kiện K) vì variance của khác biệt sẽ =0 như nhau cho mọi cá thể. Như vậy, trong trường hợp hiện tại, việc hiệu chỉnh df và kiểm tra sphericity là không cần thiết.

Phương pháp cổ điển: hàm aov

summary(aov(BPchange~step+Error(Animal), df))

## 
## Error: Animal
##           Df Sum Sq Mean Sq F value Pr(>F)
## Residuals  4  268.8   67.19               
## 
## Error: Within
##           Df Sum Sq Mean Sq F value   Pr(>F)    
## step       5   3752   750.4   60.13 2.32e-11 ***
## Residuals 20    250    12.5                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hàm aov cho phép thực hiện F test đơn giản cho mô hình BPchange~step+Error(Animal)

Giải pháp Mixed model

Những lợi ích khi tiếp cận ANOVA cho thí nghiệm lặp lại bằng mô hình hỗn hợp (Mixed model):

Khả năng tùy biến cao hơn nhiều: bạn có thể mở rộng cấu trúc mô hình thoải mái, linh hoạt để thích nghi với độ phức tạp của thiết kế thí nghiệm. Thí dụ để khảo sát đồng thời nhiều Phương pháp trị liệu, tương tác với yếu tố Thời gian, những biến số khác, phân cụm bệnh nhân… Thậm chí, kết quả của bạn cũng không bắt buộc phải là biến định lượng (bạn có thể dựng mixed logistic model luôn).

Mô hình cho phép bạn xử lý yếu tố Thời gian như một biến số định lượng liên tục thực sự. ANOVA cổ điển chỉ nhìn thời điểm như các thứ bậc của biến định tính. Mixed model có thể mô hình hóa kết quả theo thời gian một cách liên tục chứ không chỉ ước lượng trung bình tại mỗi thời điểm. Điều này rất hữu ích nếu thí nghiệm được lặp lại định kỳ tại nhiều, rất nhiều thời điểm. Điều này đồng nghĩa với việc Thời gian có thể làm covariate (hiệp biến số) trong tương tác với một yếu tố khác.

Phương pháp mô hình cho phép dung nạp dữ liệu trống (missing value) và cỡ mẫu bất cân xứng giữa các phân nhóm. Trong quy trình ANOVA cổ điển, nếu dữ liệu thiếu sót tại bất cứ thời điểm nào thì toàn bộ dữ liệu còn lại của cá thể đó sẽ bị loại bỏ. Khi dùng mixed model, chỉ thời điểm có missing value bị loại bỏ, các bậc giá trị khác của Fj hay Tj vẫn được phân tích. Mixed model xử lý dữ liệu dưới dạng “long format” chứ không phải đa biến (wide format) nên số lượng không đồng đều giữa các phân nhóm không còn là vấn đề nữa.

Đơn giản hóa việc thực hiện suy diễn thống kê, thí dụ F test (vì lúc này Mean of squares chỉ còn hai bậc: Model và Residual). Từ mô hình cũng dễ dàng kiểm định ý nghĩa thống kê của từng bộ phận của phương sai, bao gồm cả random effects và covariance matrix. Cuối cùng, tư duy bằng ngôn ngữ mô hình là điều kiện quan trọng để bạn bước qua trường phái Bayes.

Ta có thể thực hiện Mixed model bằng 2 packages trong R: nlme và lme4

nlme package

# nlme method

library(nlme)
mod1=lme(BPchange~step, random=~1|Animal, df,method='REML')
summary(mod1)

## Linear mixed-effects model fit by REML
##  Data: df 
##        AIC      BIC    logLik
##   161.0754 170.4998 -72.53768
## 
## Random effects:
##  Formula: ~1 | Animal
##         (Intercept) Residual
## StdDev:    3.019727 3.532439
## 
## Fixed effects: BPchange ~ step 
##             Value Std.Error DF   t-value p-value
## (Intercept)  1.00  2.078311 20  0.481160  0.6356
## stepD2       1.35  2.234110 20  0.604267  0.5525
## stepD3       4.60  2.234110 20  2.058985  0.0528
## stepD4      16.40  2.234110 20  7.340729  0.0000
## stepD5      26.80  2.234110 20 11.995825  0.0000
## stepD6      26.20  2.234110 20 11.727262  0.0000
##  Correlation: 
##        (Intr) stepD2 stepD3 stepD4 stepD5
## stepD2 -0.537                            
## stepD3 -0.537  0.500                     
## stepD4 -0.537  0.500  0.500              
## stepD5 -0.537  0.500  0.500  0.500       
## stepD6 -0.537  0.500  0.500  0.500  0.500
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -1.86101507 -0.40461779  0.01142447  0.55587032  1.50371349 
## 
## Number of Observations: 30
## Number of Groups: 5

anova(mod1)

##             numDF denDF  F-value p-value
## (Intercept)     1    20 82.07769  <.0001
## step            5    20 60.13340  <.0001

intervals(mod1)

## Approximate 95% confidence intervals
## 
##  Fixed effects:
##                   lower  est.     upper
## (Intercept) -3.33528009  1.00  5.335280
## stepD2      -3.31027291  1.35  6.010273
## stepD3      -0.06027291  4.60  9.260273
## stepD4      11.73972709 16.40 21.060273
## stepD5      22.13972709 26.80 31.460273
## stepD6      21.53972709 26.20 30.860273
## attr(,"label")
## [1] "Fixed effects:"
## 
##  Random Effects:
##   Level: Animal 
##                    lower     est.    upper
## sd((Intercept)) 1.285609 3.019727 7.092946
## 
##  Within-group standard error:
##    lower     est.    upper 
## 2.591129 3.532439 4.815710

VarCorr(mod1)

## Animal = pdLogChol(1) 
##             Variance  StdDev  
## (Intercept)  9.118751 3.019727
## Residual    12.478125 3.532439

lme4 package

#lme4 method

library(lme4)

mod3=lmer(BPchange~step + (1|Animal), df,REML=TRUE)

library(lmerTest)

mod3=update(mod3)
summary(mod3)

## Linear mixed model fit by REML t-tests use Satterthwaite approximations
##   to degrees of freedom [lmerMod]
## Formula: BPchange ~ step + (1 | Animal)
##    Data: df
## 
## REML criterion at convergence: 145.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.86102 -0.40462  0.01142  0.55587  1.50371 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Animal   (Intercept)  9.119   3.020   
##  Residual             12.478   3.532   
## Number of obs: 30, groups:  Animal, 5
## 
## Fixed effects:
##             Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)    1.000      2.078 12.689   0.481   0.6386    
## stepD2         1.350      2.234 20.000   0.604   0.5525    
## stepD3         4.600      2.234 20.000   2.059   0.0528 .  
## stepD4        16.400      2.234 20.000   7.341 4.29e-07 ***
## stepD5        26.800      2.234 20.000  11.996 1.37e-10 ***
## stepD6        26.200      2.234 20.000  11.727 2.04e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##        (Intr) stepD2 stepD3 stepD4 stepD5
## stepD2 -0.537                            
## stepD3 -0.537  0.500                     
## stepD4 -0.537  0.500  0.500              
## stepD5 -0.537  0.500  0.500  0.500       
## stepD6 -0.537  0.500  0.500  0.500  0.500

# Kenward-Roger denominator df method

anova(mod3,ddf="Kenward-Roger")

## Analysis of Variance Table of type III  with  Kenward-Roger 
## approximation for degrees of freedom
##      Sum Sq Mean Sq NumDF DenDF F.value    Pr(>F)    
## step 3751.8  750.35     5    20  60.133 2.315e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta có thể thấy cả 3 kết quả là tương đương. Theo đó, F test cho thấy Liều thuốc có hiệu ứng ý nghĩa làm tăng BPChange. (F(5,20)=60.133, p_value < 0.0001).

Giải pháp lai và hoàn chỉnh: package afex

afex là một package mới được tạo ra trong năm 2017, và đây là một công cụ mạnh nhất mà Nhi từng biết để thực hiện ANOVA. afex sử dụng cú pháp đơn giản và hỗ trợ tất cả thiết kế ANOVA với between, within và random effects, cho phép hiệu chỉnh vi phạm giả định sphericity bằng GG hay HF, cho phép làm posthoc testvới hiệu chỉnh ngưỡng ý nghĩa (thí dụ Bonferroni) và contrast analysis để cho ra kết quả phong phú không thua kém những software thương mại khác. Với repeated mesure ANOVA, afex dùng giải pháp lai, tức là phần lõi của nó dựa vào mixed model (package lme4), nhưng vẫn xuất ra kết quả theo truyền thống.

library(afex)

aovt=aov_ez(id="Animal",
            dv="BPchange",
            data=df,
            within=c("step"),
            anova_table = list(correction="GG",
                               es="pes",
                               p_adjust_method="bonferroni"))


summary(aovt)

## 
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
## 
##                 SS num Df Error SS den Df      F    Pr(>F)    
## (Intercept) 5514.9      1   268.76      4 82.078 0.0008227 ***
## step        3751.8      5   249.56     20 60.133 2.315e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

nice(aovt)

## Anova Table (Type 3 tests, bonferroni-adjusted)
## 
## Response: BPchange
##   Effect    df   MSE         F pes p.value
## 1   step 5, 20 12.48 60.13 *** .94  <.0001
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1
## 
## Sphericity correction method: GG

Khi kết hợp afex và package lsmean, ta có một giải pháp rất mạnh để làm post-hoc test và contrast analysis:

library(lsmeans)

lsm=lsmeans(aovt,specs = "step")

lsm

##  step lsmean       SE    df  lower.CL  upper.CL
##  D1     1.00 2.078311 12.69 -3.501122  5.501122
##  D2     2.35 2.078311 12.69 -2.151122  6.851122
##  D3     5.60 2.078311 12.69  1.098878 10.101122
##  D4    17.40 2.078311 12.69 12.898878 21.901122
##  D5    27.80 2.078311 12.69 23.298878 32.301122
##  D6    27.20 2.078311 12.69 22.698878 31.701122
## 
## Confidence level used: 0.95

lsmdf=lsm%>%summary()%>%as_tibble()


lsmdf%>%ggplot(aes(x=step, 
                   y=lsmean,
                   fill=step,
                   group=1))+
  geom_errorbar(aes(ymin=lsmean-SE, 
                    ymax=lsmean+SE), 
                width=0.2,size=1) +
  geom_line(size=1,col="red4")+
  geom_point(size=5,shape=21,col="black")+
  scale_y_continuous("BPchange")+
  theme_bw()

Thí dụ, một phân tích tương phản với 2 giả thuyết H0 : Trung bình BPchange của 3 liều D123 tương đương với của D456, và: Trung bình BPChange là như nhau giữa 3 nhóm liều: D12, D34 và D56:

Ta chỉ cần đưa trọng số contrast vào matrix và dùng hàm contrast trên object lsm:

contrast(lsm, 
         list("D123/D456"=c(-1,-1,-1,1,1,1)/3,
              "D12/D34/D56"=c(-2,-2,-1,-1,1,1)/3),
         adjust="bonferroni")

##  contrast     estimate       SE    df t.ratio p.value
##  D123/D456   21.150000 1.289864 20.00  16.397  <.0001
##  D12/D34/D56  8.433333 2.563147 10.04   3.290  0.0162
## 
## P value adjustment: bonferroni method for 2 tests

Kết quả bác bỏ cả 2 giả thuyết H0 nêu trên, như vậy có tương phản ý nghĩa giữa nhóm liều thấp và liều cao, giữa 3 nhóm liều Thấp, Trung bình và Cao.

Ta cũng có thể làm test Contrast hàm đa thức bậc cao:

contrast(lsm,method="poly")

##  contrast  estimate        SE df t.ratio p.value
##  linear      219.15 13.217176 20  16.581  <.0001
##  quadratic    18.85 14.478691 20   1.302  0.2077
##  cubic       -94.35 21.194634 20  -4.452  0.0002
##  quartic     -16.25  8.359276 20  -1.944  0.0661
##  degree 5     16.95 25.077829 20   0.676  0.5068

pcdf=contrast(lsm,method="poly")%>%coef()%>%as_tibble()

colnames(pcdf)=c("Step","linear","Quadratic","Cubic","Quartic","Deg.5")

pcdf%>%gather(linear:Deg.5,key="Poly",value="Coef")%>%
  ggplot(aes(x=Step,y=Poly,fill=Coef))+
  geom_tile(col="black")+scale_fill_gradient2(low="blue",high="red",mid="white")+
  theme_bw()

pcdf%>%gather(linear:Deg.5,key="Poly",value="Coef")%>%
  ggplot(aes(x=Step,y=Coef,col=Poly,group=1))+
  geom_path()+facet_wrap(~Poly,ncol=2)+
  theme_bw()

Kết quả Polynomial contrast analysis cho thấy khuynh hướng thay đổi của BPchange có thể được xấp xỉ bằng 1 hàm bậc 3 (cubic) nghịch biến, hoặc bậc 1 (tuyến tính).Điều này có nghĩa gì: Huyết áp tăng dần theo liều thuốc (hàm linear); tuy nhiên ở những liều thấp nhất, huyết áp chỉ tăng nhẹ, tại một liều trung gian nào đó, sẽ tăng đột biến, sau đó sẽ duy trì pha bình nguyên ở những liều cao hơn (hàm cubic)

Một loại contrast analysis khác có thể thử, đó là reversed pairwise (bắt cặp tuần tự nghịch chiều), tức là so sánh giữa liều sau và liều thấp hơn::

contrast(lsm,"revpairwise",adjust="bonferroni")

##  contrast estimate       SE df t.ratio p.value
##  D2 - D1      1.35 2.234111 20   0.604  1.0000
##  D3 - D1      4.60 2.234111 20   2.059  0.7914
##  D3 - D2      3.25 2.234111 20   1.455  1.0000
##  D4 - D1     16.40 2.234111 20   7.341  <.0001
##  D4 - D2     15.05 2.234111 20   6.736  <.0001
##  D4 - D3     11.80 2.234111 20   5.282  0.0005
##  D5 - D1     26.80 2.234111 20  11.996  <.0001
##  D5 - D2     25.45 2.234111 20  11.392  <.0001
##  D5 - D3     22.20 2.234111 20   9.937  <.0001
##  D5 - D4     10.40 2.234111 20   4.655  0.0023
##  D6 - D1     26.20 2.234111 20  11.727  <.0001
##  D6 - D2     24.85 2.234111 20  11.123  <.0001
##  D6 - D3     21.60 2.234111 20   9.668  <.0001
##  D6 - D4      9.80 2.234111 20   4.387  0.0043
##  D6 - D5     -0.60 2.234111 20  -0.269  1.0000
## 
## P value adjustment: bonferroni method for 15 tests

cdf=contrast(lsm,"revpairwise",adjust="bonferroni")%>%
  coef()%>%
  as_tibble()

cdf%>%gather(c.1:c.15,key="Pairs",value="Coef")%>%
  ggplot(aes(x=step,y=Pairs,fill=as.factor(Coef)))+
  geom_tile(col="black")+scale_fill_manual(values=c("blue","white","red"))+
  theme_bw()

Kết quả cho thấy D4 chính là liều kích hoạt sự tăng huyết áp đột biến, so với 3 liều trước đó, sau đó huyết áp tăng dần từ D4-D5 nhưng bắt đầu giảm từ D5 đến D6.

3 Giải pháp Bayes cho thí nghiệm lặp lại

Trong mô hình STAN, chúng ta phân tích bài toán hiện thời theo thứ bậc như sau:

• Giá trị của Y (BPchange) được mô tả bằng phân phối Gaussian (normal) với tham số : vị trí trung tâm Mu và scale Sigma, đây là hàm likelihood của STAN model

• Giá trị Mu[ij] được ước lượng bằng một mô hình :

\[\mu _{ij} = \beta _{j}.X_{j} + \gamma_{i}\]

4 Mô hình Bayes trong STAN

Trong đó, Fixed effect – hiệu ứng thí nghiệm được ước tính như tích của 6 tham số beta(j) với j =1,2,3,4,5,6 và design matrix Xj với 6 predictors (D1,D2,D3,D4,D5,D6); Random effect – biến thiên ngẫu nhiên cá thể được ước lượng như 5 tham số gamma(i) tương ứng với 5 con thỏ trong thí nghiệm

• Mô hình tường minh có nội dung:

\[\mu = \beta _{1}D1 + \beta _{2}D2 + \beta _{3}D3 + \beta _{4}D4 + \beta _{5}D5 + \beta _{6}D6 + \gamma_{i}\]

• Mặt khác, ta giả định rằng gamma có phân phối chuẩn với 1 tham số scale là Sigma_B

• Mục tiêu của chúng ta là tìm phân phối hậu định của : 6 tham số beta, tham số sigma, 5 tham số gamma và sigma_B

• Về Prior (tiền định), ta có thể giả định như sau: beta có phân normal (0,100); gamma có phân phối normal(0,sigma_B), sigma_B và sigma có phân phối half Cauchy

• Dữ liệu đầu vào của mô hình STAN gồm có:

1. Cỡ mẫu n= độ dài matrix long data (n=30)

2. Số dummy variables trong design matrix, hay bậc k của yếu tố can thiệp trong thí nghiệm, k=6 vì có 6 liều thuốc (mô hình này không chứa Intercept, một phần Intercept đã nằm trong tham số gammas).

3. nB = số lượng cá thể, ở đây nB=5

4. Vector biến kết quả Y, chính là BPchange

5. Design matrix X chứa các dummy variable tạo ra từ 6 bậc yếu tố can thiệp là dose step

6. Vector số nguyên B, là cột mã số định danh cho mỗi cá thể trong long data, lưu ý : vector này có độ dài = n chứ không phải nB

Ta viết cấu trúc mô hình này vào STAN code :

# STAN

rstanString="
data{
int n; //co mau
int nX; //predictors
int nB; // so block
vector [n] y;  // outcome
matrix [n,nX] X; // design matrix
int B[n]; //block vector
}

parameters{
vector [nX] beta; // fixed effect
real sigma; 
vector [nB] gamma; // random effect
real sigma_B;
}

transformed parameters {
vector[n] mu;  // fixed effect
mu = X*beta;
for (i in 1:n) {
mu[i]= mu[i] + gamma[B[i]];  //mean predicted
}

} 
model{
// Priors
beta ~ normal(0,100);
gamma ~ normal(0,sigma_B);
sigma_B ~ cauchy(0,25);
sigma ~ cauchy(0,25);

//likelihood
y ~ normal( mu , sigma );
}
"

Sau khi code mô hình, ta sẽ tạo ra dữ liệu đầu vào cho phân tích Bayes, dưới dạng 1 list, trong đó ta khai báo các thành phần mà block data đã mô tả.

Xmat=model.matrix(~step, data=df)

datalist=with(df, list(y=BPchange, X=Xmat, nX=ncol(Xmat),
                       B=as.integer(Animal),
                       n=nrow(df), 
                       nB=length(levels(Animal))))

Ta đưa data vào mô hình và kích hoạt quy trình lấy mẫu MCMC với các tùy chỉnh Ta sẽ chạy 4 chuỗi song song: Trước hết sampler sẽ chạy 1000 lượt khởi động, khi ổn định thì kết quả MCMC bắt đầu được ghi lại 2000 lượt cho mỗi chuỗi sau đó rút gọn 1/2 để cuối cùng chỉ còn 1000 giá trị x 4 chuỗi = 4000 giá trị cho mỗi tham số, và chạy song song trên tối đa 4 cores của máy tính.

library(rstan)

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

set.seed(123)

fit <- stan(data=datalist,
                           model_code=rstanString,
                           chains=4,
                           iter=3000,
                           warmup=1000,
                           thin=2,
            control=list(adapt_delta=0.99))

## 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 file3743d06efc.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"

Quy trình converged cực nhanh vì dữ liệu có n thấp. Bây giờ ta có 1 object tương đối lớn, chứa 4 chuỗi MCMCx1000 lượt cho 6 tham số beta, sigma, sigma_B, 5 tham số gamma, và 30 giá trị ước lượng cho Mu.

5 Phân phối hậu định cho Fixed và Random effects

Vì ta đang có trong tay phân phối hậu định của tất cả tham số (thành phần) của mô hình, ta có thể lần lượt khai thác chúng, bắt đầu bằng Fixed effect (hiệu ứng chính)

Fixed effect

Phân phối hậu định của 6 tham số beta (chính là trung bình của BPchange cho 6 bậc của dosestep, hay fixed effect của dose step sau khi hiệu chỉnh cho random effect) có thể được khảo sát như sau :

postdf=as.data.frame(fit)%>%as_tibble()

betadf=postdf[,c(1:6)]%>%mutate(.,Iteration=as.numeric(rep(c(1:1000),4)),
                                  Chain=as.factor(rep(c(1:4),each=1000)))

colnames(betadf)=c("D1","D2","D3","D4","D5","D6","Iteration","Chain")

head(betadf)%>%.[,c(1:6)]%>%knitr::kable()

D1	D2	D3	D4	D5	D6
-3.600280	-0.3920784	2.712608	17.35059	25.99524	26.35438
-2.990803	-1.2436725	2.734570	13.21955	25.03401	24.06605
-5.783254	1.1506698	2.456025	17.42084	27.39022	27.65323
-5.046266	5.7077594	4.030842	19.15349	30.08561	27.24923
-2.275975	-2.9764782	3.723750	12.18430	23.75162	23.44466
-1.372470	4.2986273	6.874671	18.53192	25.72957	31.57401

p1=betadf%>%gather(D1:D6,key="Step",value="Fixed_effect")%>%
  ggplot(aes(x=Fixed_effect,fill=Chain,col=Chain))+
  geom_density(alpha=0.3,show.legend = F)+
  facet_wrap(~Step,ncol=1,scales = "free_y")+
  theme_bw(6)
  
p2=betadf%>%gather(D1:D6,key="Step",value="Fixed_effect")%>%
  ggplot(aes(y=Fixed_effect,x=Iteration,col=Chain))+
  geom_path(alpha=0.5,show.legend = F)+
  facet_wrap(~Step,ncol=1,scales = "free")+
  theme_bw(6)

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

# Plot

betadf$pseudoGroup=factor(rep(c(1:20),e=nrow(betadf)/20))

betadf%>%gather(D1:D6,key="Step",value="Fixed_effect")%>%
ggplot()+
  geom_joy(aes(y=Step,
              x=Fixed_effect,
              col=pseudoGroup,
              fill=Step),
           alpha=0.03,show.legend = F,scale=1)+
  geom_vline(xintercept=median(df$BPchange),col="red3",linetype=2)+
  scale_x_continuous(breaks=c(-20,-10,0,10,20,30,40,50))+
  coord_flip()+
  scale_colour_brewer(palette = "PuRd")+
  scale_fill_brewer(palette="Reds")+
  theme_bw()

Tiếp theo, ta sẽ khảo sát Random effects, hay 5 tham số gammas cho 5 con thỏ:

Gammas cho biết sai biệt ngẫu nhiên, cơ bản và cá thể của từng con thỏ bên ngoài hiệu ứng chính của 6 liều thuốc

Random effects

# Random effects (gammas)

gammadf=postdf[,c(8:12)]%>%mutate(.,Iteration=as.numeric(rep(c(1:1000),4)),
                                Chain=as.factor(rep(c(1:4),each=1000)))

colnames(gammadf)=c("R1","R2","R3","R4","R5","Iteration","Chain")

gammadf$pseudoGroup=factor(rep(c(1:20),e=nrow(gammadf)/20))

gammadf%>%gather(R1:R5,key="Blocks",value="Random_effect")%>%
  ggplot()+
  geom_joy(aes(y=Blocks,
               x=Random_effect,
               col=pseudoGroup,
               fill=Blocks),
           alpha=0.03,show.legend = F,scale=1)+
  geom_vline(xintercept=0,col="blue3",linetype=2)+
  coord_flip()+
  scale_colour_brewer(palette = "Blues")+
  scale_fill_brewer(palette="Blues")+
  theme_bw()

Ta thậm chí có thể khảo sát phân phối hậu nghiệm cho Mu (Yhat) của 5 con thỏ tại từng bậc của dosestep:

# Y_hat

yhatdf=postdf[,c(14:43)]%>%mutate(.,Iteration=as.numeric(rep(c(1:1000),4)),
                                  Chain=as.factor(rep(c(1:4),each=1000)))

yhatdf$pseudoGroup=factor(rep(c(1:20),e=nrow(yhatdf)/20))

yhatdf%>%gather(`mu[1]`:`mu[30]`,key="Observation",value="Y_hat")%>%
  mutate(.,Step=rep(c("D1","D2","D3","D4","D5","D6"),each=20000))%>%
  ggplot()+
  geom_joy(aes(y=Observation,
               x=Y_hat,
               col=pseudoGroup,
               fill=Step),
           alpha=0.02,show.legend = F,scale=1)+
  geom_vline(xintercept=median(df$BPchange),col="blue3",linetype=2)+
  scale_fill_viridis(discrete=T,option="C",direction = -1)+
  scale_color_brewer(palette="RdPu")+
  theme_bw()+
  facet_wrap(~Step,ncol=2,scales="free_y")

6 Bảng ANOVA từ mô hình BAYES

Nội dung của bước này là tương tự như cách ta làm để dựng bảng ANOVA cổ điển, chỉ khác đó là MSM và MSR và các sum of squares không được xác định từ mô hình tuyến tính dùng algorithm Least square hay REML, mà dùng chính mô hình Bayes với 6 tham số Beta (fiex effect) và 5 tham số gamma (random effect). Để đơn giản, ta lấy Trung vị của phân phối hậu định của mỗi tham số để tạo ra 1 mô hình đại diện.

Đầu tiên, ta tính Y_hat từ mô hình đại diện như sau:

coefs=as.matrix(fit)%>%.[,c(1:6)]%>%apply(.,2, median)

gams=as.matrix(fit)%>%.[,c(8:12)]%>%apply(.,2, median)
gam=c(rep(gams[[1]],6),
      rep(gams[[2]],6),
      rep(gams[[3]],6),
      rep(gams[[4]],6),
      rep(gams[[5]],6)
      )

predDF=data.frame(x=df$step,truth=df$BPchange)

Yhat=(coefs %*% t(model.matrix(~predDF$x))+gam)%>%as.vector()

Sau đó, từ mô hình đại diện Bayes này, ta tính các SS, MSM, MSR, F values cũng như các effect-sizes.

Ghi chú: Bậc tự do của các thành phần được tính như sau:

\[df_{total}= (N-1)\]

\[df_{fix}= (k-1)\] \[df_{random}= (s-1)\]

\[df_{residual} = (N-k-s+1)\] Ở đây, Dose step là k=6 do có 6 bậc giá trị của Dosestep, s=5 cho 5 con thỏ, N=30.

Ta bắt đầu tính thủ công msm, msr, F, p_value cho F test và effectsizes là eta2 và omega2:

df1=5
df2=30-5-6+1
n=30
k=6

ssm=sum( (Yhat-mean(predDF$truth))^2 )
ssr=sum( (predDF$truth-Yhat)^2 )
sst=sum( (predDF$truth-mean(predDF$truth))^2 )
msm=ssm/df1
msr=ssr/df2
f=msm/msr

eta2=ssm/sst
r=sqrt(ssm/sst)
R2=r^2
omega2=(ssm-df1*msr)/(sst+msr) # or ((k-1)*(F-1))/((k-1)*(F-1)*(n*k))
#2-sides F test
p2sides=2*(1-(pf(f,df1,df2,F)))
#1-side F test
p1side=(1-(pf(f,df1,df2,F)))

rbind(ssm,ssr,sst,msm,msr,f,df1,df2,p1side,p2sides,eta2,omega2,r,R2)%>%round(.,3)%>%as.matrix()%>%.[,1]

##      ssm      ssr      sst      msm      msr        f      df1      df2 
## 3906.821  265.877 4270.085  781.364   13.294   58.776    5.000   20.000 
##   p1side  p2sides     eta2   omega2        r       R2 
##    0.000    0.000    0.915    0.897    0.957    0.915

Như vậy, đây là một giải pháp “lai” giữa hồi quy Bayes và null hypothesis testing. Thậm chí ta có cả p values.

Kết quả trên đây có thể xem như tương đương với kết quả của quy trình ANOVA theo trường phái Frequentist. Giá trị F gần như nhau (60 so với 58.8), p_value <0.001

7 Post-hoc test theo BAYES

Quy trình này thực sự thay thế post-hoc test theo phái frequentist bằng phương pháp Bayes. Nội dung của post hoc test theo Bayes rất đơn giản: Ta sẽ tạo ra 15 chuỗi MCMC chứa phân phối hậu định của “khác biệt trung bình” giữa 2 phân nhóm được bắt cặp tuần tự (Ta có 6 phân nhóm, nên có 15 cặp so sánh ). Suy diễn Bayes sẽ được thực hiện trực tiếp trên mỗi chuỗi MCMC, không cần dùng p_value hay hiệu chỉnh.

posthocdf=data.frame(
  D1D2=rep(NA,4000),
  D1D3=rep(NA,4000),
  D1D4=rep(NA,4000),
  D1D5=rep(NA,4000),
  D1D6=rep(NA,4000),
  D2D3=rep(NA,4000),
  D2D4=rep(NA,4000),
  D2D5=rep(NA,4000),
  D2D6=rep(NA,4000),
  D3D4=rep(NA,4000),
  D3D5=rep(NA,4000),
  D3D6=rep(NA,4000),
  D4D5=rep(NA,4000),
  D4D6=rep(NA,4000),
  D5D6=rep(NA,4000)
)

for (i in 1:4000) {
  posthocdf$D1D2[i]=betadf$D2[i]-betadf$D1[i]
  posthocdf$D1D3[i]=betadf$D3[i]-betadf$D1[i]
  posthocdf$D1D4[i]=betadf$D4[i]-betadf$D1[i]
  posthocdf$D1D5[i]=betadf$D5[i]-betadf$D1[i]
  posthocdf$D1D6[i]=betadf$D6[i]-betadf$D1[i]
  posthocdf$D2D3[i]=betadf$D3[i]-betadf$D2[i]
  posthocdf$D2D4[i]=betadf$D4[i]-betadf$D2[i]
  posthocdf$D2D5[i]=betadf$D5[i]-betadf$D2[i]
  posthocdf$D2D6[i]=betadf$D6[i]-betadf$D2[i]
  posthocdf$D3D4[i]=betadf$D4[i]-betadf$D3[i]
  posthocdf$D3D5[i]=betadf$D5[i]-betadf$D3[i]
  posthocdf$D3D6[i]=betadf$D6[i]-betadf$D3[i]
  posthocdf$D4D5[i]=betadf$D5[i]-betadf$D4[i]
  posthocdf$D4D6[i]=betadf$D6[i]-betadf$D4[i]
  posthocdf$D5D6[i]=betadf$D6[i]-betadf$D5[i]
}

posthocdf$pseudoGroup=factor(rep(c(1:20),e=nrow(posthocdf)/20))

Suy diễn Bayes được thực hiện bằng cả 2 phương pháp : ROPE (Kruschke) và Bayes Factor.

Trước hết, ta hãy xem qua phân tích ROPE phân phối hậu định của 6 khác biệt trung bình :

Nhắc lại, ROPE là một khoảng vô nghĩa thực dụng, thí dụ ta giả địng rằng nếu khác biệt trung bình nằm trong khoảng +/- 5 xem như vô nghĩa.

# Bayes Inference

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

posthocdf=posthocdf%>%gather(D1D2:D5D6,key="Pairs",value="Difference")

posthocdf%>%
  split(.$Pairs)%>%
  map(~summaryKruschke(MCMC = .$Difference,
                       compVal=0.0,
                       rope=c(-5,5),
                       credMass=0.975)
  )%>%as.data.frame()->kruschke

colnames(kruschke)=unique(posthocdf$Pairs)

kruschke%>%t()%>%round(.,2)

##       Mean Median  Mode HDIlevel    LL    UL CompVal PcntGtCompVal ROPElow
## D1D2  0.36   0.34  0.15     0.98 -9.81 10.31       0         53.48      -5
## D1D3  3.59   3.59  3.34     0.98 -6.61 14.04       0         80.22      -5
## D1D4 15.42  15.35 14.54     0.98  5.76 24.89       0         99.85      -5
## D1D5 25.76  25.79 26.09     0.98 15.88 36.53       0        100.00      -5
## D1D6 25.15  25.14 25.01     0.98 15.31 35.93       0        100.00      -5
## D2D3  3.23   3.25  3.40     0.98 -2.85  8.86       0         90.60      -5
## D2D4 15.06  14.99 14.96     0.98  9.19 20.56       0        100.00      -5
## D2D5 25.40  25.43 25.64     0.98 20.05 31.21       0        100.00      -5
## D2D6 24.79  24.74 24.19     0.98 19.06 30.38       0        100.00      -5
## D3D4 11.82  11.78 11.67     0.98  6.37 17.62       0        100.00      -5
## D3D5 22.17  22.19 22.63     0.98 16.41 28.06       0        100.00      -5
## D3D6 21.56  21.53 21.19     0.98 15.47 26.94       0        100.00      -5
## D4D5 10.35  10.34 10.08     0.98  5.01 16.23       0         99.97      -5
## D4D6  9.74   9.74  9.42     0.98  3.88 15.28       0        100.00      -5
## D5D6 -0.61  -0.65 -1.35     0.98 -6.12  5.37       0         39.40      -5
##      ROPEhigh PcntLtROPE PcntInROPE PcntGtROPE
## D1D2        5      11.12      74.67      14.20
## D1D3        5       2.85      59.30      37.85
## D1D4        5       0.00       0.88      99.12
## D1D5        5       0.00       0.00     100.00
## D1D6        5       0.00       0.00     100.00
## D2D3        5       0.08      77.47      22.45
## D2D4        5       0.00       0.03      99.97
## D2D5        5       0.00       0.00     100.00
## D2D6        5       0.00       0.00     100.00
## D3D4        5       0.00       0.55      99.45
## D3D5        5       0.00       0.00     100.00
## D3D6        5       0.00       0.00     100.00
## D4D5        5       0.00       1.58      98.42
## D4D6        5       0.00       2.97      97.03
## D5D6        5       3.25      95.22       1.52

posthocdf%>%ggplot()+
  geom_joy(aes(y=reorder(Pairs,Difference),
               x=Difference,
               fill=reorder(Pairs,Difference),
               col=pseudoGroup),
           alpha=0.05,
           show.legend = F,
           scale=1)+
  geom_vline(xintercept=c(-5,0.0,5),col="red4",linetype=2)+
  ylab("Pairs")+
  theme_bw()+ggtitle("Bayesian Pairwise comparison")+
  scale_color_brewer(palette="PuRd")+
  scale_fill_viridis(discrete=T,option="C",direction = -1)

Kết quả: Ngoại trừ cặp D1D3, D2D3,D1D2 và D5D6, tương ứng với việc chuyển từ liều thấp nhất sang liều kế cận cho đến D4, và từ liều cao nhất D5 sang 1 liều cao hơn là D6 không cho thấy khác biệt rõ nét; thì 2 cặp D3D4 và D4D5 đều cho thấy một sự tương phản khá rõ, vì mật độ phân phối hậu định của khác biệt trung bình nằm trong ROPE là 0.55% (cặp D3D4),và 1.58% (D4D5), trong khi 98.4% và 99.5% mật độ khác biệt trung bình nằm ngoài ROPE và > 5.

Kết quả này hoàn toàn tương đồng với kết quả của Contrast analysis. Chúng tôi không muốn áp đặt ý kiến riêng nào cả về kết quả này. Các bạn có toàn quyền Tin hoặc không Tin về khác biệt có ý nghĩa cho việc chuyển liều từ D3 sang D4 và D5.

Có một cách khác nếu bạn muốn phân định Trắng/Đen: đó là dùng Bayes Factor, với ngưỡng so sánh là zero (không có khác biệt), tỉ trọng chứng cứ ủng hộ cho giả thuyết H1 (khác biệt < 0):

Do có quá nhiều cặp so sánh (15), ta có thể khảo sát trực quan độ lớn của BF cho 15 ngưỡng so sánh từ -10 đến +20 trong 1 heatmap như sau:

# Bayes Factors

library(brms)

thres=c(-10:20)

threshold=rep(NA,15*length(thres))
BayesFactor=rep(NA,15*length(thres))
Pair=rep(NA,15*length(thres))

pairlev=unique(posthocdf$Pairs)

n=0

for(i in (1:15)){
  for (j in (1:length(thres))){
    tempdf=subset(posthocdf,Pairs==pairlev[i])
    Pair[n+j]=pairlev[i]
    thr=thres[j]
    threshold[n+j]=thr
    hyp=paste("Difference>",thr,sep="")
    bf=hypothesis(tempdf,hyp,alpha=0.05)
    BayesFactor[n+j]=bf$hypothesis$Evid.Ratio
  }
  n=n+length(thres)
}

tempbfdf=cbind(Pair,threshold,BayesFactor)%>%as_tibble()
tempbfdf$threshold=as.numeric(tempbfdf$threshold)
tempbfdf$BayesFactor=as.numeric(tempbfdf$BayesFactor)

for(i in (1:nrow(tempbfdf))){
  tempbfdf$BayesFactor[i]<-ifelse(tempbfdf$BayesFactor[i]>9999,9999,
                                  tempbfdf$BayesFactor[i])
}

tempbfdf%>%ggplot(aes(x=threshold,y=reorder(Pair,BayesFactor),fill=BayesFactor,col=BayesFactor))+
  geom_tile(alpha=0.7)+
  theme_bw()+
  geom_vline(xintercept=c(0,5),col="blue3")+
  scale_color_viridis(discrete=F,option="D",direction=-1)+
  scale_fill_viridis(discrete=F,option="D",direction=-1)

Kết quả cho thấy BF rất cao cho ngưỡng so sánh =0 và các cặp D4D5 (BF khoảng 3000) và D3D4 (BF > 9000)

8 Diễn đạt văn bản khoa học

Phương pháp thống kê

Hiệu ứng của liều thuốc phenylbiguanide trên sự thay đổi của huyết áp thỏ được khảo sát bằng một mô hình Bayes cho phép đo lặp lại với liều thuốc (k=6) là hiệu ứng chính, sai biệt cơ bản ở từng cá thể (i=5) là hiệu ứng ngẫu nhiên. Mô hình này cho ra phân phối hậu nghiệm của trung bình khác biệt của huyết áp khi chuyển từ liều thấp sang liều cao hơn. Suy diễn thống kê dựa vào HDI (khoảng mật độ cao nhất), ROPE (theo Kruschke) và tỉ trọng chứng cứ (Bayes factor).

Kết quả

Mô hình Bayes cho thấy liều thuốc gây ra hiệu ứng ý nghĩa làm tăng huyết áp (F(5,20)=58.78, p_value<0.0001 ; omegasquared = 0.897). Phân phối hậu định của giá trị huyết áp trung bình và trung bình khác biệt giữa các liều thuốc được trình bày trong hình … Phân tích ROPE (bảng…) và BF factor (Hình…) cho thấy khuynh hướng như sau : Huyết áp không thay đổi ý nghĩa trong phạm vi liều D1 đến D3, nhưng tăng đột ngột khi chuyển từ D3 sang D4, sau đó tăng nhẹ khi liều tăng từ D4 sang D5 và giảm dần khi chuyển từ D5 qua D6.

9 Bàn luận và Tổng kết

Tuy thí dụ minh họa trong bài chưa hoàn hảo đối với thí nghiệm ANOVA repeated measured vì không phân biệt rõ between, within effect và yếu tố thời gian; nhưng đây là một trường hợp thú vị, cho thấy ưu thế của phương pháp Bayes so với ANOVA cổ điển; vì thực ra, trên dữ liệu này ta không thể làm ANOVA được (cỡ mẫu quá thấp: n=5 cho mỗi phân nhóm, vi phạm giả định phân phối chuẩn); nhưng khi dùng Bayes không cầnbận tâm đến chuyện này. Việc quy định likelihood và prior trong mô hình Bayes cũng là giả định, nhưng nó không bắt buộc, không cố định. Nếu muốn, ta dễ dàng dùng một likelihood theo phân phối khác cho Y (thí dụ Gamma hay logNormal…), ta cũng có thể thay đổi prior tùy thích. Lý do ta chọn phân phối normal chỉ đơn giản vì đây là giả định hợp lý nhất cho biến liên tục khi ta không có một thông tin nào khác.

Đây là lần đầu tiên các bạn làm Mixed model trong STAN, vấn đề là: Khi dùng phương pháp Bayes, các bạn không cần quan tâm đến cơ chế riêng của algorithm cho từng loại mô hình (việc đó đã có MCMC sampler lo), mà chỉ cần quan tâm đến việc giải phẫu cấu trúc mô hình và mô tả cấu trúc này trong STAN. Chỉ với 2 parameters, ta có thể thêm Random effect vào model. Trong thời gian tới, nhóm BAV sẽ giới thiệu cho các bạn những mô hình còn phức tạp hơn,như thay đổi family distribution cho likelihood thành Poisson chẳng hạn.

Các thông điệp quan trọng được chuyển tải trong bài:

1. ANOVA cho thí nghiệm lặp lại là một trường hợp tổng quát cho paired sample t test

2. Cần phân biệt giữa Repeated measure và Longitudinal study, để áp dụng phù hợp.

3. Nên tiếp cận Repeated measure ANOVA bằng Mixed model thay cho giải pháp đa biến và tính thủ công. Mixed model linh hoạt hơn nhiều.

4. Phương pháp suy diễn Bayes linh hoạt và có nhiều ưu thế hơn so với Null hypothesis testing và P_value

Nhóm BAV xin chân thành cảm ơn sự ủng hộ của các bạn. Chúng tôi sẽ gặp lại các bạn trong một bài khác.