Mô hình hồi quy Negative Binomial cho nghiên cứu IVF
Lê Ngọc Khả Nhi
28 Tháng 12 năm 2022
1 Chuẩn bị thư viện cho phiên thực hành
library(tidyverse)
library(ggridges)
library(rstan)
library(brms)
library(tidybayes)
library(ggdist)
library(fitdistrplus)
library(bayestestR)
rstan::rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())2 Bước 1: Chuẩn bị dữ liệu
| Col_egg | MII_egg | Age | BMI | AFC | AMH | Infert_dt | DNG | |
|---|---|---|---|---|---|---|---|---|
| 1 | 8 | 5 | 35 | 27 | 14 | 4.81 | 39 | 0 |
| 2 | 6 | 3 | 41 | 21 | 6 | 3.56 | 22 | 0 |
| 4 | 9 | 9 | 38 | 20 | 7 | 2.42 | 0 | 0 |
| 5 | 9 | 9 | 34 | 20 | 9 | 2.28 | 48 | 0 |
| 6 | 6 | 4 | 38 | 18 | 2 | 0.34 | 61 | 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 2 phân nhóm phác đồ DNG = 1 (nhóm can thiệp) và DNG = 0 (nhóm đối chứng)
Bảng thống kê mô tả cho Col_egg
| DNG | sum | mean | sd | median | p5 | p95 |
|---|---|---|---|---|---|---|
| 0 | 676 | 10.242 | 5.836 | 9.500 | 3.000 | 20.750 |
| 1 | 433 | 6.275 | 3.609 | 6.000 | 1.400 | 12.000 |
Kiểm định Mann-Whitney
##
## Wilcoxon rank sum test with continuity correction
##
## data: Col_egg by DNG
## W = 3196, p-value = 5.044e-05
## alternative hypothesis: true location shift is not equal to 0So sánh trực quan phân phối của Col_egg giữa 2 phân nhóm
4 Bước 3: Biện luận về quy luật phân phối cho kết cục
Chứng minh phân phối Negative Binomial là phù hợp và chính xác hơn so với phân phối Poisson
## Chi-squared statistic: 465.4549 7.041656
## Degree of freedom of the Chi-squared distribution: 8 7
## Chi-squared p-value: 1.802023e-95 0.4245556
## the p-value may be wrong with some theoretical counts < 5
## Chi-squared table:
## obscounts theo Poisson theo Negative binomial
## <= 2 13 1.5693252 14.107487
## <= 3 11 3.3754144 9.736436
## <= 4 15 6.9321011 11.158186
## <= 5 11 11.3891854 11.783773
## <= 7 16 33.8928292 22.926799
## <= 8 14 18.7908628 10.306558
## <= 10 20 31.2411313 17.328319
## <= 13 14 22.2770145 17.758589
## <= 17 11 5.2473909 12.180251
## > 17 10 0.2847452 7.713602
##
## Goodness-of-fit criteria
## Poisson Negative binomial
## Akaike's Information Criterion 940.4809 802.4802
## Bayesian Information Criterion 943.3862 808.2907
# Bước 4: Thiết lập mô hình hồi quy Bayes
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. 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 (phân phối Polya):
\[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)
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 | BMI | default | |||||||
| b | DNG | default | |||||||
| b | Infert_dt | default | |||||||
| student_t(3, 2.1, 2.5) | Intercept | default | |||||||
| gamma(0.01, 0.01) | shape | 0 | default |
Ta có thể ấn định priors mới tùy ý thích:
| prior | class | coef | group | resp | dpar | nlpar | lb | ub | source |
|---|---|---|---|---|---|---|---|---|---|
| normal(0, 1) | b | DNG | NA | NA | user | ||||
| normal(0, 1) | b | Age | NA | NA | user | ||||
| normal(0, 1) | b | BMI | NA | NA | user | ||||
| normal(0, 1) | b | AFC | NA | NA | user | ||||
| normal(0, 1) | b | AMH | NA | NA | user | ||||
| normal(0, 1) | b | Infert_dt | NA | NA | user | ||||
| student_t(3, 2.1, 2.5) | Intercept | NA | NA | user | |||||
| gamma(0.01, 0.01) | shape | 0 | NA | user |
Cuối cùng, ta khớp mô hình với dữ liệu
5 Bước 5: Diễn giải kết quả mô hình
Kết quả thô của mô hình Bayes NBI:
## Family: negbinomial
## Links: mu = log; shape = identity
## Formula: Col_egg ~ DNG + Age + BMI + AFC + Infert_dt + AMH
## Data: df_mod (Number of observations: 135)
## 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 2.14 0.06 2.02 2.27 1.00 4706 2942
## DNG -0.24 0.09 -0.42 -0.06 1.00 4076 2816
## Age 0.00 0.01 -0.02 0.03 1.00 5008 3141
## BMI 0.01 0.01 -0.02 0.03 1.00 4306 2532
## AFC 0.04 0.01 0.02 0.06 1.00 3429 3160
## Infert_dt -0.00 0.00 -0.01 -0.00 1.00 5412 3395
## AMH 0.08 0.03 0.04 0.13 1.00 3264 2920
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 9.23 2.78 5.30 15.89 1.00 3931 2632
##
## 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).Ta có thể ước tính Risk-ratio cho DNG = 1 so với DNG = 0 theo cách sau:
## [1] 0.7878309 0.6453248 0.8184282Ta 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 | mad | q5 | q95 | rhat | ess_bulk | ess_tail |
|---|---|---|---|---|---|---|---|---|---|---|
| b_AFC | .value | 1.0389782 | 1.0388541 | 0.0096603 | 0.0096320 | 1.0234860 | 1.0550412 | 1.000996 | 3429.332 | 3159.878 |
| b_Age | .value | 1.0022460 | 1.0024306 | 0.0119940 | 0.0120813 | 0.9827544 | 1.0221143 | 1.001374 | 5008.034 | 3141.426 |
| b_AMH | .value | 1.0872511 | 1.0862690 | 0.0273919 | 0.0264052 | 1.0439161 | 1.1345659 | 1.000296 | 3264.167 | 2919.724 |
| b_BMI | .value | 1.0060407 | 1.0060363 | 0.0137596 | 0.0134960 | 0.9834777 | 1.0286950 | 1.000359 | 4305.999 | 2532.480 |
| b_DNG | .value | 0.7911191 | 0.7874978 | 0.0737006 | 0.0709649 | 0.6761684 | 0.9222059 | 1.001396 | 4075.577 | 2816.076 |
| b_Infert_dt | .value | 0.9963125 | 0.9963192 | 0.0016254 | 0.0016324 | 0.9936220 | 0.9989378 | 1.002325 | 5411.951 | 3394.628 |
Suy luận thống kê Bayes cho các hệ số hồi quy trong mô hình
## Summary of Posterior Distribution
##
## Parameter | Median | Mean | MAP | 95% CI | pd | ps | Rhat | ESS
## --------------------------------------------------------------------------------------------------
## (Intercept) | 2.14 | 2.14 | 2.14 | [ 2.02, 2.27] | 100% | 1.00 | 1.000 | 4721.00
## DNG | -0.24 | -0.24 | -0.24 | [-0.42, -0.06] | 99.58% | 0.93 | 1.000 | 4039.00
## Age | 2.43e-03 | 2.17e-03 | 4.52e-03 | [-0.02, 0.03] | 57.27% | 0.00 | 1.000 | 4981.00
## BMI | 6.02e-03 | 5.93e-03 | 1.21e-03 | [-0.02, 0.03] | 66.75% | 0.00 | 1.000 | 4238.00
## AFC | 0.04 | 0.04 | 0.04 | [ 0.02, 0.06] | 100% | 0.00 | 1.001 | 3418.00
## Infert_dt | -3.69e-03 | -3.70e-03 | -4.01e-03 | [-0.01, 0.00] | 99.00% | 0.00 | 1.000 | 5403.00
## AMH | 0.08 | 0.08 | 0.08 | [ 0.04, 0.13] | 100% | 0.24 | 1.000 | 3221.00
##
## # Fixed effects distributional
##
## Parameter | Median | Mean | MAP | 95% CI | pd | ps | Rhat | ESS
## ---------------------------------------------------------------------------------
## shape | 8.78 | 9.23 | 7.98 | [ 5.30, 15.89] | 100% | 1.00 | 0.999 | 3224.00Lưu ý: 1-pd sẽ cho phép ước tính giá trị p quy ước, theo kết quả này hiệu ứng của DNG là có ý nghĩa thống kê với p = 1-0.9958 = 0.0042
Vẽ biểu đồ trình bày phân phối hậu nghiệm của IRR
Trình bày trực quan marginal effect của phác đồ DNG so với phác đố
DYG
Kết thúc bài thực hành.