Mô hình hồi quy Negative Binomial lạm phát giá trị zeros

knitr::opts_chunk$set(echo = T)

1 Bài thực hành sử dụng các thư viện sau trong R:

tidyverse, ggridges, brms, tidybayes và ggdist

library(tidyverse)
library(ggridges)

library(rstan)
library(brms)

library(tidybayes)
library(ggdist)

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

2 Bước 1: Bài toán và dữ liệu

Bài toán này dựa trên nghiên cứu của Li Lu và cs “Which Factors Are Associated With Reproductive Outcomes of DOR Patients in ART Cycles: An Eight-Year Retrospective Study” công bố năm 2022 trên tạp chí Frontiers in Endocrinology. (doi: 10.3389/fendo.2022.796199).

Nhi trích xuất một phần dữ liệu để tạo ra một bài toán thống kê với mục tiêu khảo sát hiệu quả của chu kì ART sử dụng 4 phác đồ khác nhau, bao gồm Kích thích nhẹ bằng letrozole, đối vận GnRH, Follicular Phase Long-Acting protocol và Luteal Phase Short-Acting Long protocol, kèm theo các hiệp biến AMH, AFC và Tuổi, đối với kết cục số lượng noãn thu được sau chu kì kích thích BT.

Stim	Age	AMH	AFC	FSH	LH	Retrieved	MII	MII_rate	Clin_preg
1_GnRH	39	0.17	0	11.88	4.51	0	0	0	0
1_GnRH	47	0.04	4	13.69	4.89	1	1	1	0
1_GnRH	44	0.22	1	20.50	9.63	0	0	0	0
1_GnRH	42	0.08	2	14.45	3.90	1	1	1	0
1_GnRH	44	0.50	2	14.39	6.58	0	0	0	0

3 Bước 2: Thống kê mô tả

Bài toán: So sánh số lượng noãn thu được giữa 4 loại phác đồ, số trường hợp không thu được noãn nào cả và tỷ lệ của kết quả thất bại này:

Stim	sum	zeros	p0	mean	sd	median	p5	p95
0_Mild	58	41	0.466	0.659	0.829	1.000	0.000	2.000
1_GnRH	225	33	0.236	1.607	1.433	1.000	0.000	5.000
2_Short	366	13	0.092	2.577	2.317	2.000	0.000	7.950
3_Long	729	3	0.021	5.170	3.789	5.000	1.000	13.000

Kết quả cho thấy có một khuynh hướng tăng dần của số noãn trung bình, theo thứ tự : Letrozole (Mild), GnRH_ant, Luteal Phase Short-Acting Long và Follicular Phase Long-Acting. Đặc biệt, trong dữ liệu có một tần suất khá cao của những giá trị zeros (không thu được noãn nào cả)

Thực hiện 1 kiểm định phi tham số Kruskal-Wallis H và hậu kiểm (so sánh bắt cặp) bằng Wilcoxon với hiệu chỉnh Benjamini-Hogberg

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Retrieved by Stim
## Kruskal-Wallis chi-squared = 181.92, df = 3, p-value < 2.2e-16

## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  df_mod$Retrieved and df_mod$Stim 
## 
##         0_Mild  1_GnRH  2_Short
## 1_GnRH  4.1e-08 -       -      
## 2_Short < 2e-16 0.00023 -      
## 3_Long  < 2e-16 < 2e-16 2.5e-11
## 
## P value adjustment method: BH

So sánh trực quan phân phối của noãn thu được giữa 4 phân nhóm phác đồ

4 Bước 3: Mô hình hồi quy Negative Binomial lạm phát giá trị 0

Mô hình ước lượng số noãn thu được như 1 biến ngẫu nhiên Y theo phân phối Negative Binomial. lạm phát giá trị 0 (ZINBI). Phân phối ZINBI có 2 bộ phận, với bộ phận NBI : Mật độ phân phối của Y được xác định bằng 2 tham số \(\mu\) và \(\phi\) theo biến thể của hàm hỗn hợp Gamma-Poisson:

\[Pr(Y=y) = \frac{Γ(𝑦+𝜙)}{𝑦!Γ(𝜙)}(1 − \frac{𝜙}{(𝜙+𝜇)})^𝜙 (\frac{𝜇}{(𝜙+𝜇)})^𝑦\] Tham số \(\mu\) tương ứng với giá trị trung bình của Y, tham số kiểu hình \(\phi\) cho phép liên kết giữa trung bình và variance của Y theo công thức:

\[\phi = \frac{\mu^2}{Var(Y) + \mu}\] Sau đó tham số \(\mu\) được ước lượng bằng một mô hình tuyến tính tổng quát (GLM) với hàm liên kết logarit

\[log⁡(𝜇 ) = 𝛽_0+𝛽_1 𝑋_1+𝛽_2 𝑋_2+ … +𝛽_𝑗 𝑋_𝑗\]

Ghi chú: Ta có thể mô hình hóa cả 2 tham số, nhưng để đơn giản, ở đây \(\phi\) được xem như hằng số (chỉ có Intercept)

Bộ phận Zero inflated (ZI) mở rộng khả năng của phân phối NBI cơ bản bằng cách áp dụng hàm xác suất có điều kiện riêng cho trường hợp y = 0 và y > 0.

\[𝑦=0 :𝑓_𝑧 (𝑦)=𝑧+(1−𝑧)𝑓(0)\]

\[𝑦>0 :𝑓_𝑧 (𝑦)=(1−𝑧)𝑓(𝑦)\]

Với z là xác suất không thu được noãn nào cả (thất bại hoàn toàn) (y=0), và 1-z là xác suất thu được ít nhất 1 noãn (y>0).

Tham số \(z\) được ước lượng bằng một mô hình tuyến tính tổng quát với hàm liên kết là logit.

\[log⁡it(z ) = 𝛽_0+𝛽_1 𝑋_1+𝛽_2 𝑋_2+ … +𝛽_𝑗 𝑋_𝑗\]

Tiếp theo, ta dịch chuyển thang đo của các hiệp biến về vị trí trung bình = 0

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

prior	class	coef	group	resp	dpar	nlpar	lb	ub	source
	b								default
	b	AFC							default
	b	Age							default
	b	AMH							default
	b	Stim1_GnRH							default
	b	Stim2_Short							default
	b	Stim3_Long							default
student_t(3, 0.7, 2.5)	Intercept								default
gamma(0.01, 0.01)	shape						0		default
	b				zi				default
	b	Stim1_GnRH			zi				default
	b	Stim2_Short			zi				default
	b	Stim3_Long			zi				default
logistic(0, 1)	Intercept				zi				default

Cuối cùng, ta khớp mô hình với dữ liệu

5 Bước 4: Diễn giải kết quả mô hình

Kết quả thô của mô hình Bayes ZINNBI:

##  Family: zero_inflated_negbinomial 
##   Links: mu = log; shape = identity; zi = logit 
## Formula: Retrieved ~ Stim + Age + AMH + AFC 
##          zi ~ Stim
##    Data: df_mod (Number of observations: 511) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept         -0.18      0.15    -0.47     0.11 1.00     1304     1774
## zi_Intercept      -8.28      5.25   -22.04    -2.34 1.01      772      605
## Stim1_GnRH         0.79      0.17     0.45     1.13 1.00     1588     2139
## Stim2_Short        1.04      0.17     0.71     1.37 1.00     1394     1819
## Stim3_Long         1.51      0.18     1.16     1.86 1.00     1374     2015
## Age               -0.03      0.01    -0.04    -0.02 1.00     3309     2604
## AMH                0.90      0.16     0.57     1.21 1.00     2308     2217
## AFC                0.03      0.02    -0.02     0.08 1.00     2767     2258
## zi_Stim1_GnRH      1.22      6.57   -11.19    16.31 1.01      770      752
## zi_Stim2_Short     0.35      6.33   -11.34    15.18 1.00      787      673
## zi_Stim3_Long     -0.19      6.43   -12.00    14.53 1.01      740      713
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape     4.39      0.73     3.19     6.02 1.00     3026     2562
## 
## 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).

Từ kết quả của mô hình này, ta có thể thực hiện những suy diễn thống kê gồm:

• Incidence rate ratio (IRR) cho mỗi hiệp biến (Age, AFC, AMH), và mỗi phác đồ (GnRH, Short, Long) so với phác đồ tham chiếu là Mild stim.

• Ước lượng số noãn trung bình cho mỗi phác đồ, và sự khác biệt giữa 2 phác đồ bất kì (tương tự so sánh bắt cặp tuần tự trong posthoc test của phân tích ANOVA cổ điển)

• Ước tính tỉ lệ thất bại tuyệt đối (zi) cho mỗi phác đồ, và so sánh giữa 2 phác đồ bất kì

• Khảo sát hiệu ứng độc lập của AMH, Age, AFC…đối với số noãn thu được

Ta có thể áp dụng hàm exponential cho các hệ số hồi quy để ước tính Incidence rate ratio (IRR):

.variable	variable	mean	median	sd	25%	97.5%
AFC	.value	1.0309363	1.0312572	0.0242166	1.0144307	1.0786185
Age	.value	0.9702414	0.9701979	0.0062599	0.9661008	0.9824082
AMH	.value	2.4808685	2.4510229	0.4091415	2.1863523	3.3621153
GnRH	.value	2.2259410	2.1945718	0.3846597	1.9500285	3.0914756
Long	.value	4.6122198	4.5300683	0.8250409	4.0265463	6.4114664
Short	.value	2.8767170	2.8355408	0.4894839	2.5282423	3.9351100

Khác biệt về số lượng noãn trung bình (\(\mu_i\)) giữa các cặp phác đồ:

.variable	variable	mean	median	sd	25%	97.5%
GnRH_vs_Mild	.value	0.9900693	0.9929425	0.1963433	0.8644554	1.366586
Long_vs_GnRH	.value	1.9573479	1.9516406	0.3198211	1.7501867	2.608256
Short_vs_GnRH	.value	0.5342784	0.5363448	0.2251581	0.3828056	0.976336
Short_vs_Long	.value	1.4230695	1.4156581	0.2983059	1.2197591	2.025322

Xác suất thất bại hoàn toàn (zi) của 5 phác đồ

.variable	variable	mean	median	sd	25%	97.5%
p0_GnRH	.value	0.0105286	0.0027253	0.0173505	0.0001789	0.0621976
p0_Long	.value	0.0030241	0.0005943	0.0055914	0.0000373	0.0202003
p0_Mild	.value	0.0120355	0.0010459	0.0252059	0.0000253	0.0876489
p0_Short	.value	0.0040295	0.0009298	0.0074193	0.0000758	0.0254620

Khảo sát trực quan phân phối hậu nghiệm của IRR cho 4 phác đồ (so với Mild), AMH và Tuổi

Khảo sát trực quan phân phối hậu nghiệm của xác suất zi (không thu được noãn nào cả) của 4 phác đồ

Khảo sát trực quan phân phối hậu nghiệm của IRR của 4 phác đồ

Hiệu ứng của loại phác đồ với số lượng noãn thu được

Hiệu ứng của AMH với số lượng noãn thu được

Hiệu ứng của Tuổi với số lượng noãn thu được

df_mod %>%
  modelr::data_grid(AMH = modelr::seq_range(AMH, n = 5),
                    Age = modelr::seq_range(Age, n = 5),
                    AFC = modelr::seq_range(AFC, n = 5),
                    Stim = Stim,
  ) %>%
  add_epred_draws(mod_1, )%>%
  ggplot(aes(x = Age,
             y = .epred)) +
  stat_lineribbon(aes(y = .epred), 
                  fill = 'purple',
                  .width = c(.95, .75, .50), 
                  alpha = 1/8,
                  show.legend = F) +
  stat_dotsinterval(quantiles = 100,
                    alpha = 0.5, 
                    .width = c(0.75, 0.95),
                    point_interval = 'mean_qi',
                    show.legend = F)+
  scale_y_continuous(breaks = seq(0,20,1))+
  labs(y="Number of retrieved eggs", x = "Age") + 
  theme_bw(10)+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(axis.text = element_text(size = 12, color = "black"),
        axis.title = element_text(size = 15, color = "black"))+
  facet_wrap(~Stim, ncol = 2)

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