Ghi chú: Bài thực hành sử dụng các thư viện sau trong R: tidyverse, brms, tidybayes, ggdist

1 Giới thiệu

Trong bài này, Nhi áp dụng mô hình tuyến tính tổng quát với hiệu ứng hỗn hợp (Generalized Linear mixed model, GLMM) để phân tích một dữ liệu từ nghiên cứu của Qianwen Xi và cs (Thượng Hải, TQ) năm 2020 (Comparison Between PPOS and GnRHa-LongProtocol in Clinical Outcome with the First IVF/ ICSI Cycle: A Randomized Clinical Trial. Clinical Epidemiology 2020:12 261–272).

Bài toán thống kê hiện thời có mục tiêu khảo sát hiệu ứng của 2 loại phác đồ: gonadotrophin-releasing hormone agonist long protocol (GnRHa) và progestin-primed ovarian stimulation (PPOS) đối với diễn tiến của hormone LH qua 4 thời điểm: ngày 1, 6, 10 và ngày trigger trong chu kì COS.

Đầu tiên ta tải dữ liệu vào R

ID Timepoint Method LH
128 D1 PPOS 2.81
129 D1 PPOS 5.10
130 D1 PPOS 4.34
131 D1 PPOS 3.78
132 D1 PPOS 2.51

Đây là một dữ liệu có cấu trúc hỗn hợp (panel data) rất điển hình cho thí nghiệm lặp lại, trong đó một đại lượng (LH) được khảo sát tại các thời điểm/điều kiện khác nhau cho cùng cá thể. Mỗi phác đồ PPOS và GnRHa được áp dụng riêng biệt cho 128 và 127 cá thể, LH là biến kết quả chính trong bài toán, một biến số liên tục, được định lượng 4 lần tại N1, N6, N10 và ngày trigger cho mỗi cá thể.

2 Thống kê mô tả

Đây là một bài toán khá phức tạp vì ta sẽ phải khảo sát sự tương tác giữa 2 yếu tố: nhóm phác đồ (GnRHa so với PPOS) và Thời gian (4 thời điểm D1, D6, D10 và TD).

Ta thử làm một phân tích mô tả để so sánh đặc tính phân bố (trung vị, trung bình, phân vị 5 và 95) của LH tại 4 thời điểm và giữa 2 phác đồ:

Method Timepoint median p5 p95
GnRHa D1 0.380 0.190 0.889
GnRHa D6 0.525 0.222 1.207
GnRHa DX 0.640 0.248 1.317
GnRHa TD 0.525 0.240 1.768
PPOS D1 3.355 2.034 6.326
PPOS D6 2.475 0.895 5.483
PPOS DX 1.630 0.684 3.363
PPOS TD 1.720 0.554 3.944

Còn đây là hình ảnh trực quan của diễn tiến và phân phối LH:

3 Biện luận về mô hình GLMM

Trong thống kê cổ điển, bài toán này tương ứng với một phân tích phương sai hỗn hợp (Mixed ANOVA), hay phân tích phương sai 1 yếu tố cho phép đo lặp lại (Factorial repeated measure ANOVA). Tuy nhiên lần này ta sẽ áp dụng trực tiếp một công cụ khác linh hoạt và phổ quát hơn, đó là mô hình GLMM.

Trước hết, mô hình GLMM cho phép chúng ta vượt ra khỏi giới hạn của mô hình tuyến tính thông thường và ANOVA cổ điển khi ước lượng biến kết quả LH bằng một họ phân phối khác phù hợp hơn ngoài phân phối chuẫn (Gaussian).

Vì thang đo LH toàn giá trị dương (>0) và có hình ảnh phân phối lệch, quy luật phân phối thích hợp cho nó có thể là Gamma. Vì vậy ta chọn phân phối Gamma cho biến kết quả của mô hình.

Giá trị LH được ước lượng bởi 1 phân phối Gamma với 2 tham số : trung bình \(\mu\) và kiểu hình \(\alpha\)

\[Y \sim Gamma(\mu,\alpha)\] Tiếp theo, giá trị ước lượng của LH (tham số \(\mu\)) của mỗi cá thể chịu ảnh hưởng từ một số yếu tố có tính chất thứ bậc và phân nhóm như minh họa trong sơ đồ sau:

Cấu trúc mô hình

Bộ phận thứ nhất, gọi là hiệu ứng chính (fixed effect), hay hiệu ứng ở cấp độ quần thể (population effect)

Bậc hiệu ứng thứ nhất là g, hay nhóm phác đồ (GnRHa hay PPOS), trong đó một phần Intercept (\(\beta_0\)) tương ứng với nhóm GnRH, và \(\beta_{g=ppos}\) tương ứng với phần phương sai LH do PPOS độc lập gây ra.

Bậc hiệu ứng thứ 2, tương ứng với thời điểm (j), j có 4 bậc giá trị : D1, D6, DX và TD; trong đó j=D1 tham gia 1 phần vào giá trị Intercept, còn lại ta có 3 hệ số hồi quy \(\beta_{D6], \beta_{D10}, \beta{TD}\) tương ứng với hiệu ứng độc lập của GnRha tại 3 thời điểm này.

Vì ta xét sự tương tác giữa phác đồ g và thời điễm j, mô hình sẽ có thêm 3 hệ số hồi quy \(\beta_{gj}\) nữa là \(\beta{ppos:D6}, \beta_{ppos:D10}, \beta_{ppos::td}\), tương ứng với hiệu ứng độc lập của phác đồ PPOS tại 3 thời điểm D6,D10 và TD, so với D1.

Ngoài ra, mô hình của chúng ta còn đi xa hơn khi xét bộ phận thứ hai, hiệu ứng ngẫu nhiên ở cấp độ cá thể (random effects, Z). Ta sẽ dùng 1 mô hình random slope, cho phép ước lượng sự biến thiên ngẫu nhiên của cả Intercept và hiệu ứng của 4 thời điểm (biến Timepoint) cho mỗi cá thể i.

Như vậy ta có công thức mô hình ước lượng tham số \(\mu\) qua hàm liên kết logarit như sau;

\[log(\mu) \sim \beta_0 + \beta_{ppos} + \beta_{D6} + \beta_{D10} + \beta_{DT} + \beta_{ppos:D6} + \beta_{ppos:D10} + \beta_{ppos:DT} + (Time|ID)\] # Khớp mô hình với dữ liệu

Thiết lập công thức mô hình Bayes và xem nội dung giả thuyết tiền định (priors) mặc đinh.

# formula for brms model
form <- bf(LH ~ Method*Timepoint + (Timepoint|ID),
           family = Gamma(link='log'))

mod_priors = get_prior(formula = form,
                       data = df_mod)

mod_priors%>%knitr::kable()
prior class coef group resp dpar nlpar lb ub source
b default
b MethodPPOS default
b MethodPPOS:TimepointD6 default
b MethodPPOS:TimepointDX default
b MethodPPOS:TimepointTD default
b TimepointD6 default
b TimepointDX default
b TimepointTD default
lkj(1) cor default
cor ID default
student_t(3, -0.1, 2.5) Intercept default
student_t(3, 0, 2.5) sd 0 default
sd ID default
sd Intercept ID default
sd TimepointD6 ID default
sd TimepointDX ID default
sd TimepointTD ID default
gamma(0.01, 0.01) shape 0 default

ta khớp mô hình với dữ liệu

4 Diễn giải kết quả mô hình

Kết quả thô của mô hình như sau:

##  Family: gamma 
##   Links: mu = log; shape = identity 
## Formula: LH ~ Method * Timepoint + (Timepoint | ID) 
##    Data: df_mod (Number of observations: 878) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Group-Level Effects: 
## ~ID (Number of levels: 257) 
##                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                    0.35      0.03     0.29     0.42 1.03       64
## sd(TimepointD6)                  0.32      0.07     0.19     0.44 1.05       42
## sd(TimepointDX)                  0.33      0.09     0.13     0.47 1.05       43
## sd(TimepointTD)                  0.43      0.06     0.30     0.54 1.04       45
## cor(Intercept,TimepointD6)       0.38      0.28    -0.07     0.89 1.05       42
## cor(Intercept,TimepointDX)      -0.10      0.22    -0.40     0.47 1.03       58
## cor(TimepointD6,TimepointDX)     0.37      0.20    -0.12     0.66 1.02      144
## cor(Intercept,TimepointTD)       0.01      0.16    -0.26     0.35 1.03       56
## cor(TimepointD6,TimepointTD)     0.28      0.15    -0.08     0.50 1.02       96
## cor(TimepointDX,TimepointTD)     0.51      0.14     0.19     0.77 1.01      316
##                              Tail_ESS
## sd(Intercept)                     308
## sd(TimepointD6)                   201
## sd(TimepointDX)                   107
## sd(TimepointTD)                   215
## cor(Intercept,TimepointD6)        146
## cor(Intercept,TimepointDX)        143
## cor(TimepointD6,TimepointDX)      294
## cor(Intercept,TimepointTD)        273
## cor(TimepointD6,TimepointTD)      292
## cor(TimepointDX,TimepointTD)      461
## 
## Population-Level Effects: 
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                 -0.91      0.04    -0.99    -0.83 1.00      409
## MethodPPOS                 2.15      0.06     2.05     2.26 1.00     1135
## TimepointD6                0.30      0.04     0.21     0.38 1.00     1595
## TimepointDX                0.49      0.05     0.39     0.59 1.01     1524
## TimepointTD                0.36      0.05     0.26     0.46 1.00     1648
## MethodPPOS:TimepointD6    -0.64      0.06    -0.76    -0.51 1.00     1375
## MethodPPOS:TimepointDX    -1.15      0.08    -1.31    -0.99 1.00     1755
## MethodPPOS:TimepointTD    -1.12      0.07    -1.26    -0.98 1.00     1639
##                        Tail_ESS
## Intercept                  1617
## MethodPPOS                 1742
## TimepointD6                2252
## TimepointDX                2477
## TimepointTD                2358
## MethodPPOS:TimepointD6     2407
## MethodPPOS:TimepointDX     2654
## MethodPPOS:TimepointTD     2347
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape    16.01      6.01     9.60    32.07 1.06       37      122
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Suy luận thống kê của chúng ta sẽ chủ yếu dựa vào hiệu ứng chính (fixed effect hay population-level effect). Vì mô hình ước lượng \(\mu\) qua hàm liên kết logarit, nên các hệ số hồi quy trong bảng kết quả sẽ cho ra giá trị log(LH) (thang đo logarit), nếu muốn hoán chuyển về thang đo nguyên thủy, ta phải dùng hàm exponential.

Intercept tương ứng với giá trị LH ở nhóm GnRHa, thời điểm D1;

Intercept + MethodPPOS tương ứng với giá trị LH ở nhóm PPOS, thời điểm D1;

Intercept + TimepointD6 tương ứng LH ở D6 nhóm GnRHa;

Intercept + PPOS + MethodPPOS:TimepointD6 tương ứng LH ở D6 nhóm PPOS

…

Phân phối hậu nghiệm cho khác biệt giữa 2 phác đồ tại mỗi thời điểm:

.variable variable mean median sd 25% 97.5%
d1_g_vs_p .value 3.0681870 3.0648391 0.1395647 2.9739576 3.3575918
d10_g_vs_p .value 0.4443567 0.4403261 0.1082938 0.3712867 0.6643250
d6_g_vs_p .value 1.2934150 1.2888549 0.1441178 1.1924534 1.5868607
dt_g_vs_p .value 0.5763047 0.5754319 0.0326376 0.5535544 0.6428905

Phân phối hậu nghiệm cho tương phản giữa các thời điểm, cho mỗi phác đồ

Phân phối hậu nghiệm cho LH tại mỗi thời điểm, cho 2 phác đồ GnRHa và PPOS

Cuối cùng, ta trình bày trực quan hiệu ứng của phác đồ đối với diễn tiến của LH qua 4 thời điểm ước lượng bởi mô hình:

Kết luận: Khác với phác đồ GnRHa, phác đồ PPOS không đặt mục tiêu kiểm soát mức độ LH thấp liên tục kéo dài, tuy nhiên PPOS tác động làm giảm LH trong một giai đoạn ngắn từ ngày 1 đến ngày 10 trong chu kì kích thích buồng trứng, với mức thấp nhất trước ngày Trigger.

Kết thúc bài thực hành.