Phân tích hiệp phương sai theo BAYES

1 Giới thiệu

Thân chào các bạn, sau nhiều tháng gián đoạn, Nhi tiếp tục tái khởi công Project Bayes for Vietnam. Dự án này được phát động vào năm 2017 với mục tiêu 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. Đối tượng của chúng tôi là các bạn bác sĩ và sinh viên y khoa.

Qua những bài trước, Nhi và các bạn trong nhóm đã giới thiệu với các bạn cách thay thế những quy trình thống kê phổ biến như Bảng chéo, t-test, tương quan Pearson, và đặc biệt là thiết kế ANOVA (2 bài). Hôm nay, Nhi sẽ hướng dẫn các bạn mở rộng thiết kế ANOVA thành ANCOVA (phân tích hiệp phương sai) theo trường phái Bayes, thay thế cho ANCOVA cổ điển.

Do đã gián đoạn một thời gian, Nhi xin nhắc lại một số vấn đề cơ bản của phương pháp suy diễn thống kê theo trường phái Bayes như sau:

1. Nguyên tắc tổng quát cho mọi suy diễn Bayes có thể tóm tắt bằng công thức :

\[p(\theta|outcome, data)\propto p(outcome|\theta) p(\theta ,data)\]

2. Mục tiêu của chúng ta là mô tả một phân phối hậu định của một tham số Theta (vế bên trái, chính là xác suất điều kiện giá trị của theta khi có thông tin về kết quả outcome và dữ liệu data). Đại lượng này tỉ lệ với tích của hàm likelihood (xác suất có điều kiện cho phép ước tính kết quả khi có theta và dữ liệu) và một phân phối tiền định (prior= một giả thuyết về phân phối của theta trước khi ta nhìn thấy dữ liệu và kết quả).

3. Tham số theta là mục tiêu trong mọi bài toán, chúng ta chỉ cần xác định được những bộ phận còn lại (likelihood, priors) và tạo ra một mô hình phù hợp có chứa theta.

Khi đã phác thảo được mô hình trên giấy thì việc viết code không khó khăn lắm. Thí dụ chỉ nắm vững code cho mô hình GLM, ta có thể bao quát hầu hết những công cụ thống kê cổ điển, bao gồm t-test, logistic, ANOVA, ANCOVA, MANOVA, vân vân.

Do đó ta có thể thay thế tất cả những công cụ thống kê truyền thống dùng null hypothesis testing và p_value bằng phân tích Bayes.

2 Mục tiêu

Bài hôm nay, như thường lệ sẽ đi theo 3 bước như sau:

1. Trước hết, Nhi ôn lại lý thuyết về thiết kế ANCOVA, đưa ra 1 bài toán tiêu biểu và tái hiện quy trình ANCOVA cổ điển.

2. Sau đó Nhi 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.

3. Cuối cùng, Nhi 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.

3 Thí dụ minh họa

Trong bài này, Nhi sử dụng một bộ dữ liệu có thực từ một nghiên cứu pilot nhằm kiểm định một công nghệ xét nghiệm chức năng hô hấp hiện đại kết hợp vi khí dung NaCl và đo hệ số khuếch tán của 2 loại khí CO/NO qua màng phế nang/mao mạch, cho phép ước tính bề dày của lớp màng này. Như ta biết, chức năng trao đổi khí ở người phụ thuộc vào đặc tính vật lý của màng phế nang, mao mạch như tỉ lệ thuận với diện tích khả dụng, và tỉ lệ nghịch với độ dày. Một số bệnh lý hô hấp làm biến đổi những chỉ số này, thí dụ khí phế thũng (emphysema) phá hủy cấu trúc mô phế nang còn bệnh xơ phổi (fibrosis) làm tăng độ dày của màng.

Ta sẽ dùng thiết kế ANCOVA để so sánh giá trị bề dày màng phế nang mao mạch (Thickness) giữa 3 nhóm: người bình thường (Control), Khí phế thũng (Emphysema) và Xơ phổi (Fibrosis).

Nguyên nhân chúng ta phải dùng thiết kế ANCOVA chứ không phải ANOVA đơn biến, vì qua thăm dò cho thấy biến Thickness có tương quan với BMI, nên BMI được xét như một yếu tố gây nhiễu mà ta cần khắc phục. Phần tiếp theo Nhi sẽ giải thích rõ hơn về điều này.

4 Bước 1: Thăm dò dữ liệu

library(tidyverse)

dat=read.csv("https://raw.githubusercontent.com/kinokoberuji/R-Tutorials/master/Membthickness.csv",sep=";")

dat$Diagnostic%<>%as.factor()%>%
  recode_factor(.,`N`= "Control",`E`="Emphysema",`F`="Fibrosis")

data=dat

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

BMI	Diagnostic	Thickness
19.46740	Emphysema	0.1427005
31.83391	Emphysema	0.4760602
19.94450	Emphysema	0.2029516
23.43750	Emphysema	0.3968931
24.33748	Emphysema	0.2606986
27.44455	Emphysema	0.3658132

table(dat$Diagnostic)

## 
##   Control Emphysema  Fibrosis 
##        15        14         9

Như các bạn thấy, dữ liệu gồm 1 biến kết quả Thickness là số liên tục, 1 biến phân nhóm Diagnostic là biến rời rạc với 3 bậc giá trị; và BMI là biến số liên tục.

Trước hết, cỡ mẫu rất hạn chế cho từng phân nhóm, đặc biệt là nhóm Fibrosis chỉ có 9 trường hợp. Ngay cả cho một phân tích ANOVA đơn biến cổ điển, cỡ mẫu thấp là một nhược điểm rõ ràng, đây là một trong những nguyên nhân mà Nhi muốn dùng Bayes. Lưu ý rằng nếu cỡ mẫu quá thấp, suy diễn Bayes cũng có nhiều nguy cơ sai lầm không kém gì phái frequentist, vì lúc này phân phối hậu định sẽ phụ thuộc rất lớn vào prior, nếu prior sai thì ta sẽ dễ dàng suy diễn sai. Tuy nhiên, giữa một kết quả có sức mạnh thống kê chắc chắn thấp (frequentist) và một cơ hội thành công (Bayes), Nhi chọn điều thứ hai.

dat%>%group_by(Diagnostic)%>%
  summarise_at("Thickness",
               funs(Mean=mean,
               Median=median,
               Skewness=e1071::skewness(.),
               kurtosis=e1071::kurtosis(.),
               LL=quantile(.,probs=0.05),
               UL=quantile(.,probs=0.95)
               ))%>%knitr::kable()

Diagnostic	Mean	Median	Skewness	kurtosis	LL	UL
Control	0.5909805	0.5801316	0.1429082	-0.7085155	0.4638749	0.7465380
Emphysema	0.3074973	0.2893625	0.1528214	-1.5464368	0.1639136	0.4798438
Fibrosis	0.9229586	0.9901514	-0.2625314	-1.7776112	0.7236907	1.0795619

library(ggridges)

sum_df=dat%>%group_by(Diagnostic)%>%
  summarise_at("Thickness",median)

dat%>%ggplot(aes(x=Thickness,y=Diagnostic,fill=Diagnostic,col=Diagnostic))+
  geom_density_ridges(alpha=0.6,size=1,scale=1)+
  geom_point(alpha=0.3)+
  geom_point(data=sum_df,aes(x=Thickness),
             shape=23,size=5,fill="white",stroke=1.2)+
  geom_rug(alpha=0.5)+
  coord_flip()+
  scale_color_manual(values=c("#021ce5","#350091","#d8021e"))+
  scale_fill_manual(values=c("#1677ff","#a616ff","#ff1654"))+
  theme_bw()

## Picking joint bandwidth of 0.0637

Thống kê mô tả gợi ý rằng giá trị Thickness tăng cao ở nhóm Fibrosis, và giảm ở nhóm Emphysema so với nhóm Control; Tuy không thể khẳng định chắc chắn rằng Thicknes phân bố bình thường ở các phân nhóm, nhất là nhóm Fibrosis nhưng trung vị gần trung bình và skewness rất gần 0 gợi ý về tính đối xứng của phân bố.

Một phân tích ANOVA đơn biến khẳng định có sự khác biệt ý nghĩa về Thickness giữa 3 phân nhóm Diagnostic

car::Anova(lm(Thickness~Diagnostic,dat),type="III")

## Anova Table (Type III tests)
## 
## Response: Thickness
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 5.2389  1 394.995 < 2.2e-16 ***
## Diagnostic  2.0916  2  78.852 1.085e-13 ***
## Residuals   0.4642 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cor.test(dat$BMI,dat$Thickness)

## 
##  Pearson's product-moment correlation
## 
## data:  dat$BMI and dat$Thickness
## t = 3.151, df = 36, p-value = 0.003271
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1706224 0.6830964
## sample estimates:
##       cor 
## 0.4649452

Một phân tích tương quan cho thấy Thickness có quan hệ tuyến tính với BMI, yếu nhưng có ý nghĩa.

5 Ôn lại về ANOVA

Trong bài trước, chúng ta đã biết ANOVA có bản chất là một mô hình hồi quy tuyến tính có dạng Y ~ X nhằm ước lượng hiệu ứng của một yếu tố X (phân nhóm) làm thay đổi một kết quả định lượng Y.

ANOVA thường được áp dụng để so sánh, nhưng thực ra nó là một quy trình nhiều gồm nhiều bước :

1. Đầu tiên, ta bắt buộc phải kiểm tra những giả định cho thiết kế ANOVA, theo thứ tự quan trọng gồm các giả định về Phương sai đồng nhất giữa các phân nhóm, giả định về phân phối chuẩn của dữ liệu (hay của sai số phần dư của mô hình), tính độc lập cho mỗi quan sát cá thể.

2. Dựng mô hình hồi quy tuyến tính (thí dụ bằng hàm lm() của R) và kiểm định tính phù hợp dữ liệu của mô hình bằng Fisher’s F test (hàm aov() của R hoặc Anova của package car)

3. Hậu kiểm (posthoc test) hoặc phân tích mô hình với trọng số tương phản (contrast analysis), nhằm mục tiêu định vị sự khác biệt, so sánh bắt cặp tuần tự giữa các phân nhóm

4. Ước tính kích thước hiệu ứng (effect-size) của hiệu ứng (eta squared)

6 Từ ANOVA đến ANCOVA

Vì ANOVA là một mô hình hồi quy tuyến tính, ta hoàn toàn có thể đưa thêm vào mô hình những biến khác (gọi là hiệp biến – covariate) để tạo ra một mô hình hồi quy tuyến tính đa biến. Lúc này, thiết kế ANOVA được mở rộng sẽ trở thành thiết kế ANCOVA – phân tích hiệp phương sai (Analysis of covariance).

Một cách đơn giản, ANCOVA chính là một mô hình tuyến tính có cùng mục tiêu như ANOVA, nhưng chứa đồng thời biến phân nhóm X và một hiệp biến số C độc lập :

Y ~ (X + C)

Tại sao ta lại (phải) làm như thế ?

Sơ đồ sau đây sẽ giúp bạn hình dung điều gì xảy ra khi ta đưa thêm biến số C vào mô hình đơn biến của ANOVA: Như ta thấy, mô hình Y ~ X chỉ có 1 biến duy nhất, nó quá đơn giản. Sự đơn giản này có thể là vừa đủ để đáp ứng mục tiêu thực dụng: kiểm chứng hiệu ứng của thí nghiệm / so sánh giữa 3-4 phân nhóm. Nhưng đôi khi, suy nghĩ ngây thơ và chỉ nhìn thấy mục tiêu so sánh sẽ dẫn đến một kết quả F test yếu hoặc thậm chí âm tính; vì trị số F là tỉ lệ giữa phần phương sai mà mô hình cho phép giải thích (MSM) và phần còn lại mà mô hình không thể giải thích được (MSR).

Một mô hình đơn giản không bao giờ đủ để giải thích đủ thế giới hiện thực.

Khi đưa thêm 1 hiệp biến C vào mô hình, C cho phép giải thích thêm một phần phương sai và giá trị Residual sum of square (SSR) do đó sẽ thu nhỏ lại. Mô hình 2 biến mạnh hơn, phù hợp thực tế hơn.

Thiết kế ANCOVA sẽ hữu ích cho nhiều hoàn cảnh, như :

1. Khi nghiên cứu sinh muốn điều chỉnh hiệu ứng của một yếu tố gây nhiễu đã được tiên liệu trước. Yếu tố nhiễu này có thể gây sai lệch kết quả suy diễn thống kê của ANOVA một cách hệ thống.

Đa số nghiên cứu y học cơ bản đều phải dùng đến ANCOVA để hiệu chỉnh các yếu tố như tuổi, cân nặng, chiều cao cho các đại lượng sinh lý, sinh hóa vì cơ thể con người là một mạng lưới tương quan phức tạp…

Thí dụ : một nghiên cứu về chức năng tim mạch được thực hiện đồng thời ở 2 địa điểm là Cao nguyên và Đồng bằng, thì những kết quả phân tích khí máu phải được hiệu chỉnh theo độ cao, vì độ cao làm thay đổi áp suất Oxy trong khí quyển. Tương tự, ta cũng cần hiệu chỉnh các chỉ số về chức năng trao đổi khí theo giá trị Haemoglobin hoặc giới tính, vì Nam giới có nồng độ Haemoglobin trong máu cao hơn so với phụ nữ. Hầu hết những nghiên cứu định lượng trong Y học đều hiệu chỉnh kết quả theo Tuổi của bệnh nhân để loại trừ ảnh hưởng của sự lão hóa.

2. Trong các nghiên cứu thử nghiệm hiệu quả điều trị, có khi chúng ta cần chuẩn hóa giá trị cơ bản của outcome Y ở thời điểm trước khi điều trị (Baseline) cho tất cả bệnh nhân, để đảm bảo đánh giá chính xác hiệu quả điều trị của thuốc (với giả định rằng các bệnh nhân có trạng thái bệnh lý tương đương nhau khi bắt đầu điều trị).

3. Khi nghiên cứu sinh chủ quan muốn khảo sát một giả thuyết đặc biệt, nhằm xác định vai trò của nhiều yếu tố khác nhau có ảnh hưởng đến bệnh lý.

4. Nghiên cứu sinh cũng có thể chưa hài lòng với ý nghĩa thống kê và Effect-size của hiệu ứng chính trong mô hình đơn giản ban đầu và hy vọng rằng việc mở rộng mô hình sẽ giúp họ làm cho kết qủa đẹp hơn ?

Cũng như mô hình Mixed model và thiết kế thử nghiệm lâm sàng có xét random effect theo block , mục tiêu của ANCOVA cũng là thu nhỏ phần phương sai mà mô hình không giải thích được; tuy nhiên khác biệt giữa 2 phương pháp Mixed model và ANCOVA nằm ở chỗ trong Mixed model, random effect nằm bên ngoài hiệu ứng chính và phần lớn trường hợp, cơ chế của sai số ngẫu nhiên này không thể định nghĩa hay giải thích được (thí dụ cơ địa của bệnh nhân). Trong khi đó, thiết kế ANOVA chủ động sử dụng một hiệp biến số thực sự được khảo sát, được tiên liệu và có tham gia vào hiệu ứng chính trong mô hình.

7 ANCOVA theo Frequentist

Bây giờ Ta sẽ làm một phân tích ANCOVA theo trường phái frequentist trong R:

Vì ANCOVA cơ bản cũng là một thiết kế ANOVA, nó cũng yêu cầu những giả định tương tự như ANOVA. Bên cạnh đó, thiết kế ANCOVA còn có một giả định khác riêng cho nó, đó là :

1. Giả định về quan hệ tuyến tính giữa Y và hiệp biến số C : Do ta đang sử dụng hồi quy tuyến tính, lý tưởng nhất là Y tương quan tuyến tính với C. Tương quan này càng mạnh, ảnh hưởng của C đối với phương sai của Y càng lớn, thậm chí lấn át cả hiệu ứng của X. Nếu giữa Y và C không có quan hệ tuyến tính thì không có cơ sở để xét C như hiệp biến số.

2. Giả định về sự độc lập giữa X và C : Không có sự khác biệt về C giữa các phân nhóm

3. Giả định về sự đồng nhất của liên hệ tuyến tính giữa biến kết quả Y và hiệp biến số C cho từng phân nhóm X : nói cách khác, đồ thị tuyến tính Y = C gần như song song với nhau (cùng Slope) giữa các phân nhóm, hay : không có hiệu ứng do tương tác giữa C và X.

Điều kiện này không phải luôn được thỏa mãn trên thực tế, và không hoàn toàn bắt buộc – nó chỉ thể hiện cho chính khuyết điểm của thiết kế ANCOVA vì trong mô hình Y = X+C, hiển nhiên chỉ có 1 slope được ước tính cho C.

Nếu dùng mô hình GLM tổng quát, ta không còn lệ thuộc vào các thiết kế kinh điển nữa nhưng hoàn toàn có thể mở rộng mô hình chứa quan hệ tương tác Y = C+X+C :X hay Y=C*X và phân tích hiệu ứng của cả 3 thành phần X, C, C :X. Dù sao trong bài này Nhi vẫn bám sát theo thiết kế cổ điển với hy vọng là mô hình Y=C cùng slope value cho cả 3 phân nhóm.

Đầu tiên, Levene test cho phép xác nhận giả định phươn sai đồng nhất:

# Levene test
aov(Thickness ~ Diagnostic, dat)%>%car::leveneTest()

## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  2  1.1918 0.3157
##       35

Giả định Slope đồng nhất có thể kiểm tra một cách trực quan:

dat%>%ggplot(aes(x=BMI,y=Thickness,fill=Diagnostic,col=Diagnostic))+
  geom_point(alpha=0.5,size=2)+
  geom_smooth(method="lm",alpha=0.3)+
  scale_color_manual(values=c("#021ce5","#350091","#d8021e"))+
  scale_fill_manual(values=c("#1677ff","#a616ff","#ff1654"))+
  theme_bw()

Biểu đồ cho thấy Slope có vẻ giống nhau giữa Control và Fibrosis, nhưng có khác biệt ở nhóm Emphysema.

Ta cũng có thể kiểm chứng lại giả định này dựa vào mô hình có cả hiệu ứng tương tác giữa BMI và Diagnostic:

lmod0=lm(Thickness ~ Diagnostic * BMI, data)

car::Anova(lmod0)

## Anova Table (Type II tests)
## 
## Response: Thickness
##                 Sum Sq Df F value    Pr(>F)    
## Diagnostic     1.61433  2 76.2006 6.764e-13 ***
## BMI            0.07520  1  7.0989   0.01198 *  
## Diagnostic:BMI 0.05005  2  2.3625   0.11041    
## Residuals      0.33896 32                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Có thể thấy rằng không có hiệu ứng tương tác ý nghĩa giữa BMI và Diagnostic (p=0.11)

Như vậy ta có thể áp dụng ANCOVA với mô hình : Thickness ~ BMI + Diagnostic. Tuy nhiên, trước hết ta sẽ đưa trung bình BMI về 0, vì thực chất biến BMI không phải là mục tiêu mà ta quan tâm trong mô hình. Nó có mặt ở đó chỉ nhằm hỗ trợ cho F test và hiệu chỉnh cho hiệu ứng của Phân nhóm:

`

data$BMI=scale(data$BMI,scale = F)

Nhi dựng 2 mô hình lm1 và lm2 lần lượt có công thức:

Thickness ~ BMI + Diagnostic

Thickness ~ BMI + Diagnostic -1

lm1 <- lm(Thickness ~ BMI + Diagnostic + BMI, data = data)

lm2 <- lm(Thickness ~ BMI + Diagnostic -1, data = data)

summary(lm1)

## 
## Call:
## lm(formula = Thickness ~ BMI + Diagnostic + BMI, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.195395 -0.085894 -0.000105  0.071973  0.217252 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          0.625631   0.030748  20.347  < 2e-16 ***
## BMI                  0.013195   0.005147   2.564 0.014951 *  
## DiagnosticEmphysema -0.316207   0.041749  -7.574 8.49e-09 ***
## DiagnosticFibrosis   0.236581   0.058470   4.046 0.000284 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.107 on 34 degrees of freedom
## Multiple R-squared:  0.8478, Adjusted R-squared:  0.8344 
## F-statistic: 63.13 on 3 and 34 DF,  p-value: 5.555e-14

summary(lm2)

## 
## Call:
## lm(formula = Thickness ~ BMI + Diagnostic - 1, data = data)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.195395 -0.085894 -0.000105  0.071973  0.217252 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## BMI                 0.013195   0.005147   2.564    0.015 *  
## DiagnosticControl   0.625631   0.030748  20.347  < 2e-16 ***
## DiagnosticEmphysema 0.309423   0.028598  10.820 1.49e-12 ***
## DiagnosticFibrosis  0.862212   0.042811  20.140  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.107 on 34 degrees of freedom
## Multiple R-squared:  0.9735, Adjusted R-squared:  0.9704 
## F-statistic: 312.6 on 4 and 34 DF,  p-value: < 2.2e-16

Như ta thấy, ý nghĩa của 2 mô hình là khác nhau, mô hình lm1 có chứa Intercept , được hiểu như giá trị Thickness trung bình ở nhóm Control và có BMI=0, sau đó 2 dummy variable cho phân nhóm bệnh lý cho phép diễn giải: Thickness giảm 0.316 µm ở nhóm Emphysema và tăng 0.236 µm ở phân nhóm Fibrose.

Mô hình thứ hai không chứa Intercept, và hệ số hồi quy cho mỗi Dummy variable cũng chính là giá trị Thickness trung bình ở từng phân nhóm.

Trong cả 2 mô hình, hệ số hồi quy cho BMI như nhau = 0.013, cho biết Thickness sẽ tăng 0.013 µm cho mỗi đơn vị BMI.

Khi áp dụng hàm Anova, ta có kết quả của F test cho riêng BMI và Diagnostic. Kết quả cho thấy sau khi hiệu chỉnh cho hiệu ứng của BMI, Yếu tố bệnh lý có hiệu ứng độc lập làm thay đổi bề dày màng phế nang: F(3,34)=400.64; p_value < 2.10^-16

car::Anova(lm1, type = "III")

## Anova Table (Type III tests)
## 
## Response: Thickness
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 4.7368  1 413.9951 < 2.2e-16 ***
## BMI         0.0752  1   6.5722   0.01495 *  
## Diagnostic  1.6143  2  70.5467 7.935e-13 ***
## Residuals   0.3890 34                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(lm2, type = "III")

## Anova Table (Type III tests)
## 
## Response: Thickness
##             Sum Sq Df  F value  Pr(>F)    
## BMI         0.0752  1   6.5722 0.01495 *  
## Diagnostic 13.7520  3 400.6443 < 2e-16 ***
## Residuals   0.3890 34                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Marginal effect của Phân nhóm bệnh lý có thể trình bày như sau:

lmod <- lm(Thickness ~ as.vector(BMI)+Diagnostic, data = data)

newdata <- data_frame(Diagnostic = levels(data$Diagnostic), 
                      BMI = mean(data$BMI, na.rm = TRUE))

adj_df <- predict(lmod, newdata = newdata, interval = "confidence")

adj_df=cbind(newdata,adj_df)

ggplot(adj_df,aes(y = fit, x = Diagnostic))+
  geom_pointrange(aes(ymin = lwr, ymax = upr,col=Diagnostic),size=1.2) +
  geom_line(aes(group=1),linetype=2,size=1.2)+
  scale_y_continuous("Thickness")+ 
  scale_x_discrete("Diagnostic")+
  theme_bw()+
  scale_color_manual(values=c("#021ce5","#350091","#d8021e"))

Lưu ý, giá trị thể hiện trên biểu đồ chính là hiệu ứng riêng phần của từng bệnh lý, sau khi hiệu chỉnh cho BMI.

Còn đây là hiệu ứng riêng phần của BMI

newdata <- expand.grid(Diagnostic = levels(data$Diagnostic), 
                       BMI = seq(min(data$BMI),max(data$BMI),l=100))

pred_df <- predict(lmod, newdata = newdata, interval = "confidence")

pred_df <- cbind(newdata,pred_df)

obs_df <- cbind(data, obs = fitted(lmod) + resid(lmod))

ggplot(pred_df, aes(y = fit, x = BMI, group = Diagnostic))+ 
  geom_point(data = obs_df,aes(y = obs, shape = Diagnostic,col=Diagnostic),size=1.5)+ 
  geom_line(aes(col=Diagnostic),size=1)+
  geom_ribbon(aes(ymin = lwr,ymax = upr,group=Diagnostic,fill=Diagnostic),alpha = 0.2)+
  scale_y_continuous("Thickness")+
  scale_x_continuous("BMI")+
  theme_bw()+
  scale_color_manual(values=c("#021ce5","#350091","#d8021e"))+
  scale_fill_manual(values=c("#1677ff","#a616ff","#ff1654"))

Trong bài trước (Anova cho phép đo lặp lại), Nhi có giới thiệu về package afex rất hiệu quả cho phân tích phương sai (mọi thiết kế ANOVA). Ta sẽ thử áp dụng nó cho thiết kế ANCOVA như sau:

library(afex)

dat%>%mutate(id=rownames(.))%>%
  aov_ez(id="id",
         dv="Thickness", 
         data=.,
         between="Diagnostic",
         covariate="BMI",
         factorize = F,
         anova_table = list(es="pes",
                            p_adjust_method="bonferroni"))->cov1

Kết quả của afex “đẹp” hơn, vì ta có áp dụng một hiệu chỉnh Bonferroni cho p_value,ngoài F test , afex còn tính cho ta cả effect-size bộ phận cho mỗi yếu tố (partial eta squared: pes). Theo đó, Bệnh lý có kích thước hiệu ứng cao hơn nhiều so với BMI (pes=0.805 sv 0.162). Kết quả F test cho Diagnostic là F(2,34)=70,55; p< 1.5*10^-12

summary(cov1)

## Anova Table (Type 3 tests, bonferroni-adjusted)
## 
## Response: Thickness
##            num Df den Df      MSE       F     pes    Pr(>F)    
## Diagnostic      2     34 0.011442 70.5467 0.80582 1.587e-12 ***
## BMI             1     34 0.011442  6.5722 0.16199    0.0299 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ta còn có thể trình bày kết quả trung bình và Ci của Thickness sau khi hiệu chỉnh cho BMI. So với kết quả thực tế (bên dưới), kết quả sau hiệu chỉnh có một số khác biệt nhỏ, thí dụ 95%CI mở rộng hơn, trung bình cao hơn

lsm=lsmeans(cov1,specs = "Diagnostic")

lsm

##  Diagnostic    lsmean         SE df  lower.CL  upper.CL
##  Control    0.6094171 0.03887345 32 0.5302345 0.6885997
##  Emphysema  0.3114069 0.02753148 32 0.2553271 0.3674867
##  Fibrosis   0.9062413 0.04952127 32 0.8053698 1.0071128
## 
## Confidence level used: 0.95

dat%>%group_by(Diagnostic)%>%
  summarise_at("Thickness",
               funs(Mean=mean,
               LL=quantile(.,probs=0.05),
               UL=quantile(.,probs=0.95)
               ))%>%knitr::kable()

Diagnostic	Mean	LL	UL
Control	0.5909805	0.4638749	0.7465380
Emphysema	0.3074973	0.1639136	0.4798438
Fibrosis	0.9229586	0.7236907	1.0795619

Một lần nữa, ta lại có thể vẽ biểu đồ Marginal effect của Diagnostic

lsm%>%summary()%>%as_tibble()%>%
  ggplot(aes(x=Diagnostic, 
             y=lsmean,
             fill=lsmean,
             group=1))+
  geom_errorbar(aes(ymin=lsmean-SE, 
                    ymax=lsmean+SE), 
                width=0.2,size=1) +
  geom_line(size=1,col="grey",linetype=2)+
  geom_point(size=5,shape=21,col="black",show.legend = F)+
  scale_y_continuous("Thickness")+
  scale_fill_gradient(low="blue",high="red")+
  theme_bw()

package afex còn cho phép làm phân tích tương phản (contrast analysis) rất dễ dàng:

Thí dụ một giả thuyết tương phản với trọng số 0,-1/3,+1/3, hoặc tương phản theo hàm đa thức bậc 2

contrast(lsm, 
         list("Cont1"=c(0,-1,1)/3),
         adjust="bonferroni")

##  contrast  estimate         SE df t.ratio p.value
##  Cont1    0.1982781 0.01888662 32  10.498  <.0001

contrast(lsm,method="poly")

##  contrast   estimate         SE df t.ratio p.value
##  linear    0.2968242 0.06295634 32   4.715  <.0001
##  quadratic 0.8928446 0.08363869 32  10.675  <.0001

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

colnames(pcdf)=c("Step","linear","Quadratic")

pcdf%>%gather(linear:Quadratic,key="Poly",value="Coef")%>%
  ggplot(aes(x=Step,y=Poly,fill=Coef))+
  geom_tile(col="black")+
  scale_fill_gradient(low="white",high="black")+
  theme_bw()

Ta cũng có thể làm một post-hoc test với hiệu chỉnh Bonferroni :

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

##  contrast               estimate         SE df t.ratio p.value
##  Emphysema - Control  -0.2980102 0.04763536 32  -6.256  <.0001
##  Fibrosis - Control    0.2968242 0.06295634 32   4.715  0.0001
##  Fibrosis - Emphysema  0.5948344 0.05665985 32  10.498  <.0001
## 
## P value adjustment: bonferroni method for 3 tests

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

colnames(cdf)=c("Group","E/C","F/C","F/E")

cdf%>%gather(`E/C`:`F/E`,key="Pairs",value="Coef")%>%
  ggplot(aes(x=Group,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 có sự khác biệt ý nghĩa về giá trị Thickness giữa 2 bệnh lý Emphysema và Fibrosis so với nhóm chứng, cụ thể: Bệnh nhân khí phế thũng có bề dày màng phế nang giảm trung bình -0.298 µm so với người bình thường (p<0.0001) , trong khi bệnh xơ phổi làm tăng bề dày màng trung bình +0.297 µm (p=0.0001)

8 ANCOVA theo BAYES

Bây giờ, chúng ta sẽ thực hiện ANCOVA bằng phương pháp Bayes:

Phân tích Bayes cho thiết kế ANCOVA được bắt đầu với giả định: Giá trị quan sát Yik ở mỗi cá thể i trong phân nhóm Xk là kết quả của một biến ngẫu nhiên Y có phân phối chuẩn (Gaussian) được xác định bằng 2 tham số : Mu (trung bình) và Sigma (độ lệch chuẩn).

Giá trị tham số Mu được ước tính bằng một mô hình tuyến tính có nội dung:

\[\mu _{ik} = \beta _{c}C + \beta _{1}X_{1} + \beta _{2}X_{2} + \beta _{3}X_{3}\]

Giả thuyết tiền định (prior) cho beta và sigma gồm:

\[\beta \sim \mathit{N} (0,\sigma _{\beta })\]

\[\sigma \sim \mathit{Cauchy} (0,5)\]

Trước hết ta viết nội dung mô hình bằng STAN code

# stancode

model_code = "
// data block

data {
int<lower=1> n;    // sample size
int<lower=1> nX;   // factor level
vector [n] y;      // response vector Y
matrix [n,nX] X;   // design matrix of predictors
}

// para block

parameters {
vector[nX] beta;   // beta coefs = vector
real<lower=0> sigma;  // sigmaY = real number
}
transformed parameters { // MuY = vector
vector[n] mu;         

mu = X*beta;  // Estimate Mu
}

// Model block

model {

// Likelihood : y follows Gaussian dist.law, determined by sigma and mu

y~normal(mu,sigma);

// Priors for beta and sigma

beta ~ normal(0,10);
sigma~cauchy(0,5);
}

// Supp values block

generated quantities {
vector[n] log_lik;

for (i in 1:n) {
log_lik[i] = normal_lpdf(y[i] | mu[i], sigma); 
}
}
"

STAN code được lưu dưới dạng chuỗi kí tự (string),

model.text <- “…..”

Một chương trình STAN được chia theo cấu trúc từng khối (block) và mỗi block có vai trò khác nhau:

1. Block data: Tên gọi, loại và kích thước của từng phần trong dữ liệu đầu vào

• n là một số nguyên dương chỉ cỡ mẫu

• nX là 1 số nguyên dương, chỉ số predictors trong design matrix của mô hình (lưu ý: không có intercept, Diagnostic bị chia thành 3 dummy variables)

• y là biến kết quả (Thickness), 1 vector số thực có độ dài = n

• X là 1 design matrix có kích thước (n,nX) : lưu ý; tính kế thừa của STAN code, biến sau sử dụng thông tin của biến trước đã được khai báo

2. Block parameter:

Như tên gọi của nó, đây là nơi khai báo về các tham số mà ta cần trong mô hình. Như đã trình bày, ta có 2 loại tham số: beta (trung bình mỗi phân nhóm, hay tham số hồi quy trong mô hình), là một vector số thực có kích thước nX = số cột trong design matrix; còn sigma là sd của residual, là 1 số thực.

3. Block Transformed parameters:

Block này không bắt buộc phải có cho mọi mô hình vì nó tùy thuộc vào mục đích của người viết code và độ phức tạp của mô hình. Cần chú ý là mọi tham số tạo ra trong block này cũng sẽ được lấy mẫu cho chuỗi MCMC và xuất ra kết quả mô hình như các tham số trong block parameter.

Trong block này, ta khai báo tham số Mu là trung bình dự báo của Y, xác định từ beta và design matrix X.

4. Block Model: Trong block này ta sẽ khai báo về thành phần likelihood và priors (giả thuyết tiền định về phân phối của từng tham số).

Trong mô hình này, likelihood là y ~ normal (mu, sigma);

Ta có 2 tham số là beta và Sigma. Nếu ta không khai báo prior, STAN sẽ mặc định dùng prior Uniform cho tham số đó. Ở đây ta dùng prior “non informative” là half-Cauchy cho sigma, và 1 giả thuyết tiền định rằng beta có phân phối chuẩn, dao động quanh trung bình = 0 và sd=10).

5. Ngoài ra, ta còn có thể dùng thêm block thứ 5 để lấy mẫu cho các tham số phụ hay phái sinh từ các parameter có sẵn trong mô hình; tuy nhiên việc này không thực sự quan trọng, nhất là việc tính toán nhiều tham số bằng cách lấy mẫu sẽ làm quy trình tạo MCMC quá tải một cách không cần thiết.

Sau khi viết xong STAN code, ta sẽ tạo design matrix cho mô hình và data cho mô hình là 1 list

Xmat <- model.matrix(~BMI+Diagnostic-1, data)

data_list <- with(data, list(y = Thickness, 
                           X = Xmat, 
                           nX = ncol(Xmat), 
                           n = nrow(data)))

data_list

## $y
##  [1] 0.1427005 0.4760602 0.2029516 0.3968931 0.2606986 0.3658132 0.4868705
##  [8] 0.1753361 0.4272775 0.3836717 0.2935065 0.2852184 0.1908612 0.2171024
## [15] 1.0247708 0.7409633 0.8390008 0.8316315 0.7121756 1.0190631 1.0488030
## [22] 1.1000678 0.9901514 0.5403678 0.5984456 0.6568714 0.5801316 0.6411785
## [29] 0.5188145 0.7352794 0.5209054 0.6217350 0.6701464 0.4871185 0.5732800
## [36] 0.4096397 0.7728080 0.5379854
## 
## $X
##            BMI DiagnosticControl DiagnosticEmphysema DiagnosticFibrosis
## 1  -4.85173804                 0                   1                  0
## 2   7.51477071                 0                   1                  0
## 3  -4.37463707                 0                   1                  0
## 4  -0.88163933                 0                   1                  0
## 5   0.01834039                 0                   1                  0
## 6   3.12540572                 0                   1                  0
## 7   3.69108440                 0                   1                  0
## 8  -1.84793056                 0                   1                  0
## 9   0.74016770                 0                   1                  0
## 10  1.97934851                 0                   1                  0
## 11  0.78719338                 0                   1                  0
## 12 -1.26738442                 0                   1                  0
## 13 -1.13574832                 0                   1                  0
## 14 -5.54079466                 0                   1                  0
## 15 10.91631112                 0                   0                  1
## 16  5.61913228                 0                   0                  1
## 17 -0.16032136                 0                   0                  1
## 18 12.40933868                 0                   0                  1
## 19 -0.88163933                 0                   0                  1
## 20  6.47424311                 0                   0                  1
## 21  3.69127262                 0                   0                  1
## 22  1.56143131                 0                   0                  1
## 23  1.80330965                 0                   0                  1
## 24 -0.93045277                 1                   0                  0
## 25 -2.61121594                 1                   0                  0
## 26 -3.10082213                 1                   0                  0
## 27 -6.39092390                 1                   0                  0
## 28 -1.29345472                 1                   0                  0
## 29 -0.93478559                 1                   0                  0
## 30 -4.16292870                 1                   0                  0
## 31  2.84135450                 1                   0                  0
## 32 -5.41223831                 1                   0                  0
## 33 -2.71420106                 1                   0                  0
## 34 -3.66159933                 1                   0                  0
## 35 -3.11603410                 1                   0                  0
## 36 -6.31359916                 1                   0                  0
## 37  1.10087528                 1                   0                  0
## 38 -2.68949057                 1                   0                  0
## attr(,"assign")
## [1] 1 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$Diagnostic
## [1] "contr.treatment"
## 
## 
## $nX
## [1] 4
## 
## $n
## [1] 38

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 2 chuỗi song song: Trước hết sampler sẽ chạy 500 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,rút gọn 1/2 kích thước mỗi chuỗi 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 = data_list, 
          model_code = model_code , 
          chains = 2,
          iter = 2500, 
          warmup = 500, 
          thin = 2)

## In file included from C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
##                  from C:/Users/Admin/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
##                  from file35ec56f4195.cpp:8:
## C:/Users/Admin/Documents/R/win-library/3.5/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 converge rất nhanh, kết quả xuất ra là những chuỗi MCMC cho 4 tham số beta (tương ứng trung bình của 3 phân nhóm và beta cho BMI), 1 tham số sigma và 38 giá trị dự báo cho Mu của từng cá thể

Ta tóm tắt kết quả của 4 tham số beta và tham số sigma như sau:

library(broom)
tidyMCMC(fit, conf.int = TRUE, conf.method = "HPDinterval", pars = c("beta", "sigma"))

##      term   estimate   std.error    conf.low  conf.high
## 1 beta[1] 0.01324467 0.005430898 0.001801404 0.02333015
## 2 beta[2] 0.62519093 0.032571757 0.564763063 0.69212344
## 3 beta[3] 0.31030750 0.030162298 0.248693968 0.36638148
## 4 beta[4] 0.86053624 0.045440858 0.779930691 0.95755980
## 5   sigma 0.11104550 0.014152905 0.086291494 0.13969127

Ta khảo sát trực quan phân phối hậu nghiệm của beta cho 3 phân nhóm

postdf=fit%>%as.data.frame()%>%
  as_tibble()%>%.[,c(1:5)]%>%
  mutate(.,Iteration=as.numeric(rep(c(1:1000),2)),
         Chain=as.factor(rep(c(1:2),each=1000)))

colnames(postdf)=c("BMI","Control","Emphysema","Fibrosis","Sigma","Iteration","Chain")

postdf$Pseudo=factor(rep(c(1:10),e=nrow(postdf)/10))

library(ggridges)

postdf%>%gather(Control,Emphysema,Fibrosis,
                key="Predictors",value="Effect")%>%
  ggplot()+
  geom_density_ridges_gradient(aes(y=reorder(Predictors,Effect),
                                   x=Effect,fill=-..x..),
                               linetype=1,col="black",
                               alpha=0.6,scale=1,
                               show.legend = F)+
  geom_vline(xintercept=adj_df$fit[1],col="red3",linetype=2)+
  theme_bw()+
  scale_y_discrete("Factors")+
  scale_x_continuous("Estimated Thickness")+
  scale_fill_gradient(low="#f90046",
                       high="#f9cc00")+
  coord_flip()

postdf%>%gather(Control,Emphysema,Fibrosis,
                key="Predictors",value="Effect")%>%
  ggplot()+
  geom_density_ridges(aes(y=reorder(Predictors,Effect),
                          x=Effect,col=Pseudo),
                      linetype=1,fill="red",
                      alpha=0.03,scale=1,
                      show.legend = F)+
  geom_vline(xintercept=adj_df$fit[1],col="red3",linetype=2)+
  theme_bw()+
  scale_y_discrete("Factors")+
  scale_x_continuous("Estimated Thickness")+
  coord_flip()+ 
  scale_colour_brewer(palette = "Reds")

Hoặc phẩm chất của các chuỗi MCMC:

p1=postdf%>%gather(BMI:Sigma,key="Factor",value="Effect")%>%
  ggplot(aes(x=Effect,fill=Chain,col=Chain))+
  geom_density(alpha=0.3,show.legend = F)+
  facet_wrap(~Factor,ncol=1,scale="free_y")+
  theme_bw(8)

p2=postdf%>%gather(BMI:Sigma,key="Factor",value="Effect")%>%
  ggplot(aes(y=Effect,x=Iteration,col=Chain))+
  geom_path(alpha=0.5,show.legend = F)+
  facet_wrap(~Factor,ncol=1,scale="free_y")+
  theme_bw(8)

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

Từ mô hình Bayes, ta có thể trình bày Marginal effect của BMI điều kiện hóa theo phân nhóm:

library(coda)

mcmc = as.matrix(fit)

newdata=expand.grid(BMI = seq(min(data$BMI), max(data$BMI),len=100),
                    Group = levels(data$Diagnostic))

Xmat2 = model.matrix(~BMI + Group-1, newdata)
coefs = mcmc[, c("beta[1]", "beta[2]", "beta[3]", "beta[4]")]
pred = coefs %*% t(Xmat2) 

newdata = newdata %>% cbind(tidyMCMC(pred, conf.int = TRUE, conf.method = "HPDinterval"))

newdata%>%
  ggplot(aes(x=BMI,
             y=estimate,
             fill=Group,
             col=Group))+
  geom_ribbon(aes(ymin = conf.low,
                  ymax = conf.high), alpha = 0.2,col=NA)+ 
  geom_line()+
  geom_smooth(alpha=0.5)+
  scale_y_continuous("Thickness")+
  theme_bw()+
  scale_color_manual(values=c("#021ce5","#350091","#d8021e"))+
  scale_fill_manual(values=c("#1677ff","#a616ff","#ff1654"))

Mô hình này có khả năng giải thích đến 83.48 % phương sai của Thickness (Ci=0.785 đến 0.867)

mcmc <- as.matrix(fit)
Xmat = model.matrix(~BMI+Diagnostic-1,data)
wch = grep("beta", colnames(mcmc))
coefs = mcmc[, wch]
pred = coefs %*% t(Xmat)
resid = sweep(pred, 2, data$Thickness, "-")
var_f = apply(pred, 1, var)
var_e = apply(resid, 1, var)
R2 = var_f/(var_f + var_e)
tidyMCMC(as.mcmc(R2), conf.int = TRUE, conf.method = "HPDinterval")

##   term  estimate  std.error  conf.low conf.high
## 1 var1 0.8340875 0.02476514 0.7873155 0.8673647

Tiếp theo, Nhi sẽ làm 1 phân tích post-hoc theo Bayes. Nội dung của quy trình này là tạo ra các chuỗi MCMC cho khác biệt bắt cặp tuần tự giữa 3 phân nhóm

posthocdf=data.frame(
  CvsE=rep(NA,2000),
  CvsF=rep(NA,2000),
  EvsF=rep(NA,2000)
)

for (i in 1:2000) {
  posthocdf$CvsE[i]=postdf$Control[i]-postdf$Emphysema[i]
  posthocdf$CvsF[i]=postdf$Fibrosis[i]-postdf$Control[i]
  posthocdf$EvsF[i]=postdf$Fibrosis[i]-postdf$Emphysema[i]}

Sau đó ta có thể khảo sát những chuỗi MCMC này theo nhiều cách, hoặc dùng phản nghiệm (một loại Z-test), có thể xuất ra giá trị p

library(coda)
mcmcpvalue <- function(samp) {
  if (length(dim(samp)) == 0) {
    std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - mean(samp),
                     transpose = TRUE)
    sqdist <- colSums(std * std)
    sum(sqdist[-1] > sqdist[1])/length(samp)
  } else {
    std <- backsolve(chol(var(samp)), cbind(0, t(samp)) - colMeans(samp),
                     transpose = TRUE)
    sqdist <- colSums(std * std)
    sum(sqdist[-1] > sqdist[1])/nrow(samp)
  }
  
}

mcmcpvalue(posthocdf$CvsE)

## [1] 0

mcmcpvalue(posthocdf$CvsF)

## [1] 0

mcmcpvalue(posthocdf$EvsF)

## [1] 0

Kết quả cho thấy p_value cho cả 3 cặp so sánh đều rất thấp, cho thấy khác biệt có ý nghĩa thống kê

Một cách làm khác là suy diễn bằng ROPE và CompVal như Kruschke đã làm : Các bạn xem những bài trước để hiểu cơ chế của phép suy diễn này.

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

longpost=posthocdf%>%gather(CvsE,CvsF,EvsF,key="Pairs",value="Difference")

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

colnames(kruschke)=unique(longpost$Pairs)

kruschke%>%knitr::kable()

	CvsE	CvsF	EvsF
Mean	0.3148834	0.2353453	0.5502287
Median	0.3147589	0.2375636	0.5504648
Mode	0.3148886	0.2536895	0.5539946
HDIlevel	0.9750000	0.9750000	0.9750000
LL	0.2183531	0.1037730	0.4313252
UL	0.4166552	0.3791614	0.6757934
CompVal	0.0000000	0.0000000	0.0000000
PcntGtCompVal	100.0000000	100.0000000	100.0000000
ROPElow	0.0000000	0.0000000	0.0000000
ROPEhigh	0.1000000	0.1000000	0.1000000
PcntLtROPE	0.0000000	0.0000000	0.0000000
PcntInROPE	0.0000000	1.6000000	0.0000000
PcntGtROPE	100.0000000	98.4000000	100.0000000

Tóm tắt: CompVal là 1 ngưỡng vô nghĩa, thí dụ khác biệt = 0; ROPE là một khoảng vô nghĩa, gồm 2 đầu, thí dụ 0 và +0.1. Kết quả cho thấy 100% phân phối hậu nghiệm của khác biệt lớn hơn 0 cho 2 cặp Control vs Emphysema và Control vs Fibrosis, 99.95% phân bố hậu nghiệm lớn hơn 0 cho cặp Emphysema vs Fibrosis.

Tương tự: 100% phân bố hậu nghiệm nằm ngoài ROPE cho 2 cặp C vs E, Cvs F; và 98.45% nằm ngoài ROPE cho cặp E vsF.

Đơn giản hơn là mô tả trực quan phân phối hậu nghiệm của 3 cặp so sánh, người đọc có thể tự cảm nhận mức độ khả tín so với ngưỡng 0

longpost%>%
  ggplot()+
  geom_density_ridges_gradient(aes(y=reorder(Pairs,Difference),
                                   x=Difference,fill=-..x..),
                               linetype=1,col="black",
                               alpha=0.6,scale=1,
                               show.legend = F)+
  geom_vline(xintercept=0,col="red3",linetype=2)+
  theme_bw()+scale_y_discrete("Pairs")+scale_x_continuous("Difference")+
  scale_fill_gradient(low="#f90046",high="#f9cc00")+coord_flip()

Cuối cùng là suy diễn bằng Bayes Factor hay tỉ trọng chứng cứ phủ quyết cho 1 ngưỡng vô nghĩa nhất định từ 0 đến 0.7

threshold=rep(NA,33)
BayesFactor=rep(NA,33)
Pair=rep(NA,33)

thres=c(0,0.05,0.1,0.15,0.2,0.25,0.3,0.4,0.5,0.6,0.7)

pairlev=unique(longpost$Pairs)

n=0

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

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

tempbfdf%>%
  ggplot(aes(x=threshold,y=BayesFactor,fill=BayesFactor))+
  geom_path()+
  geom_point(show.legend = F,size=3,shape=21,col="black")+
  theme_bw(10)+
  geom_text(aes(label=round(tempbfdf$BayesFactor,2)),col="black",show.legend = F,
            angle = 50,nudge_y=BayesFactor+300,nudge_x=0,size=3)+
  geom_hline(yintercept =500,linetype=2,col="red")+
  facet_wrap(~Pair,ncol=1,scales = "free")+
    scale_fill_gradient(low="#f90046",high="#f9cc00")

tempbfdf%>%knitr::kable()

Pair	threshold	BayesFactor
CvsE	0.00	Inf
CvsE	0.05	Inf
CvsE	0.10	Inf
CvsE	0.15	Inf
CvsE	0.20	221.2222222
CvsE	0.25	13.4927536
CvsE	0.30	1.7322404
CvsE	0.40	0.0288066
CvsE	0.50	0.0000000
CvsE	0.60	0.0000000
CvsE	0.70	0.0000000
CvsF	0.00	Inf
CvsF	0.05	499.0000000
CvsF	0.10	61.5000000
CvsF	0.15	10.6279070
CvsF	0.20	2.6166365
CvsF	0.25	0.7241379
CvsF	0.30	0.1621150
CvsF	0.40	0.0040161
CvsF	0.50	0.0000000
CvsF	0.60	0.0000000
CvsF	0.70	0.0000000
EvsF	0.00	Inf
EvsF	0.05	Inf
EvsF	0.10	Inf
EvsF	0.15	Inf
EvsF	0.20	Inf
EvsF	0.25	Inf
EvsF	0.30	Inf
EvsF	0.40	221.2222222
EvsF	0.50	4.7971014
EvsF	0.60	0.2202563
EvsF	0.70	0.0035123

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

Phần này nhằm trình bày một văn bản khoa học mẫu cho phân tích ANCOVA theo trường phái Bayes.

Phương pháp: Một mô hình ANCOVA được thực hiện nhằm khảo sát thay đổi của bề dày màng PNMM ở các bệnh nhân khí phế thũng và xơ phổi so với nhóm chứng, sau khi hiệu chỉnh cho hiệu ứng của BMI. Các tham số trong mô hình được xác định bằng phương pháp Bayes với giả định kết quả có phân bố chuẩn,và prior cho Mu và sigma lần lượt là Normal(0,10) và Cauchy (0,5). Sự khác biệt giữa các phân nhóm được kiểm định bằng phương pháp suy diễn Bayes theo John Kruschke và Bayes Factor.

Kết quả:

Thăm dò dữ liệu cho thấy bề dày màng phế nang có quan hệ tuyến tính với BMI nhưng không có hiệu ứng tương tác giữa BMI và phân nhóm bệnh lý. Sau khi hiệu chỉnh cho BMI, kết quả mô hình ANCOVA Bayes cho thấy có sự thay đổi ý nghĩa về bề dày màng phế nang ở bệnh nhân có bệnh hô hấp so với nhóm chứng. Cụ thể, bệnh Khí phế thũng làm màng phế nang mỏng đi trung bình 0.32 µm (HDI:-0.20 đến 0.41 µm; BayesFactor>100 cho ngưỡng vô hiệu = -0.20 µm), trong khi bệnh xơ phổi làm màng phế nang dày thêm trung bình 0.24 µm (HDI:0.10 đến 0.38 µm; Bayesfactor =63 và 499 tương ứng với ngưỡng vô hiệu = +0.1 và +0.05 µm).

10 Kết luận

Do hầu hết những biến số trong nghiên cứu y học đều có mối tương quan phức tạp, thiết kế ANCOVA rất quan trọng. Trong bài thực hành này, các bạn đã làm quen với nó (nếu các bạn chưa biết) và thực hiện được một phân tích ANCOVA bằng cả 2 phương pháp Frequentist và Bayes. Bản chất của ANCOVA cũng như ANOVA, là một mô hình hồi quy tuyến tính nhưng có xét thêm hiệu ứng của hiệp biến số C mà trung bình đã được chuyển về zero. Mô hình này cho phép hiệu chỉnh tác động riêng của C đối với kết quả để đánh giá chính xác hơn về hiệu ứng chính của phân nhóm X. Sự hiệu chỉnh này có thể mang lại nhiều bất ngờ, thậm chí thay đổi kết quả suy diễn thống kê cho X trong một số trường hợp.

Chúc các bạn thực hành vui !