Hồi quy Binomial

1 Giới thiệu

Những nghiên cứu về quy trình thụ tinh ống nghiệm đặt ra yêu cầu khảo sát nhiều loại kết cục lâm sàng với bản chất khác nhau. Mặt khác, có tính chất kế thừa trong chuỗi kết quả của quy trình. Thành phẩm thu được của công đoạn sau chính là hệ quả từ việc tiến hành một loạt phép thử (thí dụ thụ tinh, chuyển phôi…) trên tập hợp kết quả công đoạn trước. Thí dụ, số lượng phôi là hệ quả từ tỷ lệ thụ tinh thành công trên tập hợp nhiều noãn. Do đó, cùng một thực thể kết cục nhưng có thể được diễn đạt theo nhiều cách; số lượng phôi, tỷ lệ thụ tinh trên tổng số noãn, xác suất thụ tinh thành công cho mỗi noãn…

Hầu hết những nghiên cứu hiện nay chỉ mới dừng lại ở sự so sánh và kiểm định bằng những phép kiểm thống kê thô sơ như Mann-Whitney, Student t test, cho mọi kết quả là con số, không phân biệt bản chất của chúng la số đếm, số liên tục hay tỷ lệ.

Trong chương này, chúng tôi giới thiệu một phương pháp phân tích hợp lý và chính xác hơn cho kết cục là tỷ lệ, áp dụng lý thuyết xác suất cho biến rời rạc và mô hình hồi quy tuyến tính tổng quát với phân phối nhị thức (Binomial).

Tình huống minh họa như sau: Một bác sĩ muốn khảo sát hiệu quả tạo phôi (blastocytes) của kỹ thuật kích thích trưởng thành noãn kép (Dual trigger), so với 2 kỹ thuật tham chiếu dùng agonist và hCG đơn giản.

2 Kế hoạch phân tích

Khả năng tạo blastocyte khi dùng kỹ thuật Dual trigger so với 2 phương pháp tham chiếu (đồng vận hoặc hCG đơn) sẽ được ước lượng bằng marginal Odds-ratio (OR) và risk-ratio (RR), có hiệu chỉnh cho Tuổi và AFC, dựa vào một mô hình hồi quy tuyến tính tổng quát (GLM) với phân phối nhị thức (Binomial).

Suy diễn thống kê dựa vào phủ nhận giả thuyết vô hiệu (OR và RR = 1) ở ngưỡng ý nghĩa p = 0.05.

3 Công cụ cần thiết cho quy trình

• Thư viện dplyr trong hệ sinh thái tidyverse để thao tác dữ liệu và thống kê mô tả

• Thư viện ggplot2, ggrides, tidybayes và ggdist để vẽ một số biểu đồ thống kê;

• Thư viện gamlss để dựng mô hình GLM với phân phối Binomial

• Thư viện marginaleffects để ước tính OR, RR và suy diễn thống kê từ kết quả mô hình.

library(tidyverse)
library(ggridges)

library(gamlss)
library(marginaleffects)

library(tidybayes)
library(ggdist)

Đầu tiên, ta tải dữ liệu từ file Dual_trigger.csv vào dataframe df bằng hàm read.csv, sau đó áp dụng hàm factor để chuyển dạng biến Trigger thành số thứ tự (1 = a (agonist), 2 = d (dual), 3 = h (hCG)).

Ta tạo ra thêm biến p1 là tỷ lệ tạo blastocytes từ số noãn ban đầu, bằng cách chia số lượng blastocytes (biến blast) cho số noãn (biến oocytes).

df = read.csv('Dual_trigger.csv', sep = ';', 
              dec = ',', 
              fileEncoding = 'UTF-8-BOM')%>%na.omit()

df$Trigger = factor(df$Trigger)
df$p1 = df$blast/df$oocytes

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

Trigger	Age	AFC	oocytes	blast	p1
a	41.49863	2	1	1	1.0000000
h	33.70411	3	3	2	0.6666667
d	34.57534	3	3	1	0.3333333
a	36.05205	3	2	1	0.5000000
h	25.75616	3	3	1	0.3333333
a	28.35616	4	4	1	0.2500000

4 Bước 1: Phân tích mô tả

Ta kết hợp hàm group_by và summarize của thư viện dplyr để dựng một bảng thống kê mô tả đơn giản, trình bày trung vị, SD, phân vị thứ 5 và 95 của tỷ lệ tạo blastocytes (biến p1) trong 3 phân nhóm kỹ thuật trigger.

df%>%group_by(Trigger)%>%
  summarize(n = n(),
            mean = sprintf("%0.3f", mean(p1)),
            SD = sprintf("%0.3f", sd(p1)),
            p5 = sprintf("%0.3f", quantile(p1, 0.05)),
            p95 = sprintf("%0.3f", quantile(p1, 0.95)),
            )%>%knitr::kable()

Trigger	n	mean	SD	p5	p95
a	1263	0.278	0.181	0.062	0.600
d	273	0.320	0.214	0.077	0.706
h	259	0.290	0.196	0.077	0.667

Theo kết quả này, mặc dù giá trị trung bình của tỷ lệ tạo blastocyte có vẻ cao hơn ở nhóm d (dual trigger), nhưng giới hạn phân bố của tỷ lệ này lại không cho thấy sự khác biệt giữa 3 nhóm.

Sau đây, ta sẽ so sánh trực quan đặc tính phân phối của số lượng blastocytes (bằng biểu đồ tần suất) và tỷ lệ blastocyte (bằng biểu đồ boxplot).

Hình 1: Phân bố của số lượng blastocytes giữa 3 nhóm kỹ thuật trigger

pals = c("d" = "#fc035e", 
         "h" = "#b814ff", 
         "a" = "#fcb103")


df %>% ggplot(aes(x = blast, 
                      y = Trigger, 
                      fill= Trigger)) + 
  geom_density_ridges(stat = "binline", 
                      scale = 1, 
                      bins = 50,
                      alpha = 0.7,
                      draw_baseline = FALSE) +
  labs(x="Number of Blastocytes", y = "Trigger type") + 
  scale_fill_manual(values = pals, name = "Trigger") +
  scale_x_continuous(breaks = seq(0,30,5))+
  coord_flip() +
  theme_bw(10) +
  geom_vline(xintercept = median(df$blast),
             linetype=2,
             col="blue")+
  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"))

Chú thích: 3 biểu đồ tần suất trình bày phân bố của số lượng blastocytes (trục Y) ở 3 nhóm kỹ thuật trigger (trục X).

Hình 2: Phân bố của tỷ lệ tạo blastocytes từ số noãn ban đầu giữa 3 nhóm

df %>% ggplot(aes(y = p1, 
                      x = Trigger, 
                      fill= Trigger)) + 
  geom_boxplot(alpha = 0.7) +
  labs(y="Blastocyte formation rate", x = "Trigger type") + 
  scale_fill_manual(values = pals, name = "Trigger") +
  scale_y_continuous(breaks = seq(0,1,0.2))+
  coord_flip() +
  theme_bw(10) +
  geom_vline(xintercept = median(df$p1),
             linetype=2,
             col="blue")+
  theme(axis.text = element_text(size = 12, color = "black"),
        axis.title = element_text(size = 15, color = "black"))

Chú thích: 3 biểu đồ boxplot trình bày trung vị, phân vị thứ 25,75 của tỷ lệ tạo blastocytes (trục X).Những điểm tròn tương ứng với các giá trị vượt xa khỏi phân vị thứ 95.

5 Bước 2: Dựng mô hình hồi quy Binomial

Hiện thời, kết cục mà ta quan tâm là số lượng \(k\) blastocytes thu được sau quá trình thụ tinh từ tập hợp \(n\) noãn ban đầu, và tỷ lệ \(p = k/n\). Ta thấy bài toán này tương đồng với định nghĩa về quy luật phân phối xác suất Binomial.

Thực vậy, phân phối Binomial (nhị thức) cho phép ước lượng giá trị của một biến ngẫu nhiên \(Y\), với ý nghĩa \(Y\) là số lần thu được kết quả thành công từ \(n\) lượt thử nghiệm độc lập, biết rằng mỗi thử nghiệm có xác suất thành công là \(p\).

\[Y \sim B(n,p)\] Xác suất \(Y\) nhận một giá trị = k được ước tính bởi hàm:

\[Pr(k;n,p) = \binom{n}{k}p^{k}(1-p)^{n-k}\]

Với các tham số có ý nghĩa như sau:

\(k\) là số kết quả thành công (ở đây là số lượng blastocyte tạo thành)

\(n\) chỉ số lượt thử nghiệm (ở đây là số noãn được mang đi thụ tinh)

\(n \in \left\{0,1,2...\right\}\)

\(p\) là xác suất thành công của một thử nghiệm (xác suất tạo được blastocyte cho mỗi noãn)

\(p \in \left [ 0,1\right ]\)

\((1-p)\) là xác suất của kết quả đối nghịch (thất bại, không tạo được blastocyte)

Thay vì dùng tham số n và p, ta có thể đặt ra tham số \(\mu = \frac{k}{n}\) như tỷ lệ thu được blastocyte trên tổng số noãn ban đầu cho mỗi đối tượng bệnh nhân. Ta có thể ước lượng giá trị của \(\mu\) bằng một mô hình tuyến tính tổng quát với hàm liên kết logit:

\[log⁡it(\mu)=log(\frac{\mu}{1-\mu}) = \beta_0+𝛽_1 𝑋_1+𝛽_2 𝑋_2+ … +𝛽_𝑗 𝑋_𝑗\] Ta sẽ dùng mô hình này để suy diễn thống kê về hiệu ứng của kỹ thuật dual trigger trên tỷ lệ tạo blastocytes.

Để dựng mô hình binomial, ta dùng hàm gamlss của thư viện gamlss, với một số lưu ý như sau:

• Công thức của mô hình hồi quy binomial khác với các mô hình thông thường, vì ở vị trí kết quả y, ta cần kết hợp 2 biến: số lượng thành công (blast) và số lượng thất bại (oocytes - blast), bằng hàm cbind

• Tùy chỉnh family = BI(mu.link = “logit”)

Trong mô hình hiện thời, biến độc lập là Trigger, ta cho nó tương tác với 2 hiệp biến khác là Age và AFC:

“cbind(blast, oocytes-blast) ~ Trigger * (Age + AFC)”

Nội dung kết quả thô của mô hình như sau

bi_mod = gamlss(cbind(blast, oocytes-blast) ~ 
               Trigger * (Age + AFC),
             data = df,
             family = BI(mu.link = "logit"))

## GAMLSS-RS iteration 1: Global Deviance = 7646.78 
## GAMLSS-RS iteration 2: Global Deviance = 7646.78

summary(bi_mod)

## ******************************************************************
## Family:  c("BI", "Binomial") 
## 
## Call:  gamlss(formula = cbind(blast, oocytes - blast) ~ Trigger *  
##     (Age + AFC), family = BI(mu.link = "logit"), data = df) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  logit
## Mu Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.0593442  0.1152770  -9.190  < 2e-16 ***
## Triggerd      0.3033475  0.3965697   0.765  0.44442    
## Triggerh      0.5956986  0.3634420   1.639  0.10138    
## Age          -0.0031340  0.0036684  -0.854  0.39304    
## AFC           0.0017986  0.0009032   1.991  0.04661 *  
## Triggerd:Age -0.0008936  0.0113729  -0.079  0.93738    
## Triggerh:Age -0.0091438  0.0106586  -0.858  0.39107    
## Triggerd:AFC -0.0165921  0.0076311  -2.174  0.02981 *  
## Triggerh:AFC -0.0218670  0.0076811  -2.847  0.00447 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1795 
## Degrees of Freedom for the fit:  9
##       Residual Deg. of Freedom:  1786 
##                       at cycle:  2 
##  
## Global Deviance:     7646.78 
##             AIC:     7664.78 
##             SBC:     7714.215 
## ******************************************************************

Lưu ý rằng mô hình này ước lượng giá trị của \(\mu\) là 1 xác suất (hay tỷ lệ) và sử dụng hàm liên kết là logit, vì vậy các hệ số hồi quy mà ta thấy đang ở thang đo \(log(Odds)\), với odds là tỉ lệ giữa 2 xác suất \(\frac{\mu}{1-\mu}\), ta có thể chuyển về thang đo odds bằng hàm exponential, hoặc chuyển hẳn về xác xuất \(\mu\) bằng hàm plogis.

hệ số (Intercept) tương ứng với \(log(\mu/(1-\mu))\) khi kỹ thuật trigger là Agonist (nhóm tham chiếu), và Age, AFC ở vị trí trung bình của quần thể.

Áp dụng hàm plogis cho intercept, ta sẽ có giá trị trung bình của tỷ lệ blastocyte ở phân nhóm Agonist = 0.257 ở điều kiện Age = 30.55 và AFC = 16.24

plogis(-1.0593442)

## [1] 0.2574348

Triggerd và Triggerh tương ứng với giá trị log(odds) tăng thêm ở 2 nhóm Dual trigger và hCG so với nhóm tham chiếu Agonist (dĩ nhiên trong điều kiện AFC và Age không đổi).

Ở đây ta quan tâm đến nhóm Dual trigger, hệ số hồi quy của nó = 0.3033475, giá trị này > 0 nên ta có thể ước lượng là nhóm Dual trigger có khả năng tạo blastocytes cao hơn so với nhóm tham chiếu Agonist;

Áp dụng hàm plogis ta có thể ước tính tỷ lệ tạo blastocyte trung bình ở nhóm Dual trigger = 0.32

plogis(-1.0593442 + 0.3033475)

## [1] 0.3195161

Tuy nhiên, để suy diễn thống kê một cách dễ dàng hơn, ta sẽ tiến hành quy trình ước tính marginal effects như sau đây.

6 Bước 3: Suy diễn thống kê

Đầu tiên, áp dụng hàm avg_comparisons, ta có thể ước tính được khác biệt về xác suất tạo blastocytes trung bình giữa 2 kỹ thuật: A so với D, A so với H và D so với H.

dif_p = avg_comparisons(
  bi_mod,
  what = "mu",
  variables = list(Trigger = 'pairwise'))%>%
  dplyr::select(-c(1,6))

dif_p %>% knitr::kable()

term	contrast	estimate	std.error	p.value	conf.low	conf.high
Trigger	d - a	0.0029597	0.0103186	0.7742380	-0.0172643	0.0231838
Trigger	h - a	-0.0040992	0.0090471	0.6504782	-0.0218313	0.0136328
Trigger	h - d	-0.0070590	0.0128231	0.5819844	-0.0321919	0.0180739

Theo kết quả này, cả 3 cặp so sánh đều không có ý nghĩa thống kê (p value > 0.05 và KTC95% chứa giá trị 0).

Ta cũng có thể ước tính hiệu ứng của D so với A và H thông qua trị số Odds-ratio bằng cách áp dụng 2 hàm transform_pre = “lnoravg” và post = exp.

or = avg_comparisons(
  bi_mod,
  what = "mu",
  variables = list(Trigger = 'pairwise'),
  transform_pre = "lnoravg",
  transform_post = "exp")%>%
  dplyr::select(-c(1))

or %>% knitr::kable()

term	contrast	estimate	p.value	conf.low	conf.high
Trigger	ln(odds(d) / odds(a))	1.0160647	0.7734235	0.9115953	1.132506
Trigger	ln(odds(h) / odds(a))	0.9779573	0.6512904	0.8878552	1.077203
Trigger	ln(odds(h) / odds(d))	0.9624951	0.5810925	0.8402912	1.102471

Theo kết quả này, nếu chọn nhóm Agonist làm tham chiếu, OR của Dual trigger = 1.016, tuy nhiên KTC95% chứa giá trị 1 và p value > 0.05, nên không có ý nghĩa thống kê.

Tương tự, OR của nhóm hCG so với Dual trigger = 0.96 < 1, tuy nhiên cũng không có ý nghĩa thống kê.

Ta còn có thể ước tính hiệu ứng bằng Risk-ratio (RR), chỉ cần thay hàm transform_pre = “lnratioavg”.

rr = avg_comparisons(
  bi_mod,
  what = "mu",
  variables = list(Trigger = 'pairwise'),
  transform_pre = "lnratioavg",
  transform_post = exp)%>%
  dplyr::select(-c(1))

rr %>% knitr::kable()

term	contrast	estimate	p.value	conf.low	conf.high
Trigger	ln(mean(d) / mean(a))	1.0120817	0.7730269	0.9327663	1.098141
Trigger	ln(mean(h) / mean(a))	0.9832669	0.6516765	0.9138047	1.058009
Trigger	ln(mean(h) / mean(d))	0.9715292	0.5806721	0.8768966	1.076374

Ta diễn giải kết quả RR tương tự như trên, không có sự khác biệt giữa Dual trigger và nhóm tham chiếu Agonist hoặc hCG.

Ngoài ra, ta có thể đi xa hơn khi ước lượng hiệu ứng (OR hoặc RR) của Dual trigger trong những điều kiện khác nhau của hiệp biến, thí dụ AFC = 2,5,10 và 20 như sau:

or = avg_comparisons(
  bi_mod,
  what = "mu",
  by = "AFC",
  variables = list(Trigger = 'pairwise'),
  newdata = datagrid(AFC = c(2,5,10,20),
                     grid_type = 'counterfactual'),
  transform_pre = "lnoravg",
  transform_post = "exp")

or %>% knitr::kable()

type	term	contrast	AFC	estimate	p.value	conf.low	conf.high
response	Trigger	ln(odds(d) / odds(a))	2	1.2749243	0.0115669	1.0558637	1.539434
response	Trigger	ln(odds(d) / odds(a))	5	1.2130216	0.0133671	1.0409297	1.413565
response	Trigger	ln(odds(d) / odds(a))	10	1.1164571	0.0499268	1.0000352	1.246433
response	Trigger	ln(odds(d) / odds(a))	20	0.9457769	0.4721554	0.8124331	1.101006
response	Trigger	ln(odds(h) / odds(a))	2	1.3143693	0.0056288	1.0831184	1.594993
response	Trigger	ln(odds(h) / odds(a))	5	1.2309708	0.0089146	1.0534511	1.438405
response	Trigger	ln(odds(h) / odds(a))	10	1.1035630	0.0691433	0.9923044	1.227296
response	Trigger	ln(odds(h) / odds(a))	20	0.8869339	0.0901621	0.7720058	1.018971
response	Trigger	ln(odds(h) / odds(d))	2	1.0309391	0.8183369	0.7948992	1.337069
response	Trigger	ln(odds(h) / odds(d))	5	1.0147971	0.8896481	0.8246449	1.248796
response	Trigger	ln(odds(h) / odds(d))	10	0.9884509	0.8716356	0.8585411	1.138018
response	Trigger	ln(odds(h) / odds(d))	20	0.9377834	0.5286725	0.7679211	1.145219

Kết quả này khá thú vị, vì nó cho thấy rằng khi AFC <= 5, OR của Dual trigger so với Agonist lần lượt là 1.27 và 1.21 và có ý nghĩa thống kê.

Ta có thể trình bày trực quan về hiệu ứng của kỹ thuật trigger đối với tỷ lệ tạo blastocyte, dựa vào mô hình Binomial như sau:

preds = predictions(bi_mod, 
                    what = "mu", 
                    newdata = datagrid(Trigger = c('d','a','h'),
                                       grid_type = "counterfactual"),
                    new_data = df)

preds%>%ggplot()+
  stat_halfeye(alpha = 0.4, 
               aes(x = Trigger, 
                   y = estimate,
                   fill = Trigger),
               .width = c(0.75, 0.95),
               point_interval = 'median_qi',
               show.legend = T)+
  stat_lineribbon(aes(y = estimate, x = Trigger), 
                  fill = 'gold',
                  .width = c(.95, .75, .50), 
                  alpha = 1/8,
                  show.legend = F) +
  labs(y="Rate of blastocytes", x = "Trigger") + 
  scale_fill_manual(values = pals, name = "Trigger")+
  scale_y_continuous(breaks = c(seq(0,0.4,0.05)))+
  theme_bw()

Ta có thể trình bày cùng kết quả trên nhưng cho riêng 6 điều kiện khác nhau về giá trị AFC = 1,5,10,15,20 và 30

preds = predictions(model = bi_mod,
                    what = "mu",
                    newdata = datagrid(AFC = c(1,5,10,15,20,30),
                                       grid_type = "counterfactual"))

preds%>%ggplot()+
  stat_halfeye(alpha = 0.4, 
               aes(x = Trigger, 
                   y = estimate,
                   fill = Trigger),
               .width = c(0.75, 0.95),
               point_interval = 'median_qi',
               show.legend = T)+
  stat_lineribbon(aes(y = estimate, x = Trigger), 
                  fill = 'gold',
                  .width = c(.95, .75, .50), 
                  alpha = 1/8,
                  show.legend = F) +
  labs(y="Rate of blastocytes", x = "Trigger") + 
  scale_fill_manual(values = pals, name = "Trigger")+
  facet_wrap(~AFC, scales = "free")+
  theme_bw()

7 Diễn giải kết quả của phân tích

Kết quả phân tích hồi quy cho thấy ở cấp độ quần thể, kỹ thuật dual trigger có hiệu quả tương đương với 2 ký thuật Agonist và hCG về tỷ lệ tạo blastocytes. Tuy nhiên, ở những đối tượng có AFC <5, hiệu quả tạo blastocyte của Dual trigger cao hơn một cách có ý nghĩa so với Agonist (OR = 1.27, KTC95%: 1.05 - 1.54, p=0.01).