-
Notifications
You must be signed in to change notification settings - Fork 1
/
03-NBI.Rmd
676 lines (481 loc) · 29.3 KB
/
03-NBI.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
# Khảo sát kết cục số đếm: Hồi quy Negative Binomial
## Giới thiệu
...
## Bối cảnh của thí nghiệm
...
## Kế hoạch phân tích
...
## Công cụ cần thiết
Quy trình phân tích sẽ sử dụng các thư viện sau trong R:
+ Thư viện dplyr trong hệ sinh thái tidyverse[@tidyverse] để thao tác dữ liệu và thống kê mô tả
+ Thư viện fitdistrplus[@fitdistrplus] để đối chiếu 3 quy luật phân phối Gaussian, Poisson và NBI, nhằm chứng minh tính phù hợp hơn của phân phối NBI.
+ Thư viện ggplot2, ggrides[@ggridges], tidybayes[@tidybayes] và ggdist[@ggdist] để vẽ một số biểu đồ thống kê; thư viện patchwork để ghép các biểu đồ với nhau.
+ Thư viện gamlss[@gamlss1] để dựng mô hình GLM với phân phối NBI
+ Thư viện marginaleffects[@marginal] để ước tính avegrage marginal effect (AME) và suy diễn thống kê từ kết quả mô hình.
+ Thư viện broom[@broom] nhằm trích xuất dữ liệu kết quả ước lượng từ mô hình.
```{r}
library(tidyverse)
library(ggridges)
library(tidybayes)
library(ggdist)
library(patchwork)
library(fitdistrplus)
library(gamlss)
library(broom)
library(marginaleffects)
```
...
```{r,message = FALSE,warning=FALSE}
df = read.csv('POSEIDON_ART.csv',
sep = ';',
dec = ',',
fileEncoding = "UTF-8-BOM")%>%
dplyr::filter(POSEIDON>0)
df$POSEIDON= as.factor(df$POSEIDON)
df%>%dplyr::sample_frac(0.1)%>%
head(6)%>%
knitr::kable(digits = 3)
```
Giải thích:
Dữ liệu được đọc từ một file CSV có tên 'POSEIDON_ART.csv', với ; làm dấu phân cách và , làm dấu thập phân.
Các dòng có giá trị POSEIDON nhỏ hơn hoặc bằng 0 sẽ bị loại bỏ.
Chuyển đổi kiểu dữ liệu: Biến POSEIDON được chuyển thành kiểu factor.
Chọn mẫu ngẫu nhiên và hiển thị dữ liệu:
Lấy ra một mẫu ngẫu nhiên từ df chiếm 10% dữ liệu.
Hiển thị 6 dòng đầu tiên của mẫu này dưới dạng bảng với 3 chữ số sau dấu phẩy.
## Thống kê mô tả
...
```{r,message = T,warning=T}
tot = nrow(df)
df %>%
group_by(POSEIDON) %>%
summarize(n = n(),
prop = 100*n/tot,
mean = mean(Collected, na.rm = T),
sd = sd(Collected, na.rm = T),
median = median(Collected, na.rm = T),
p5 = quantile(Collected, 0.05, na.rm = T),
p95 = quantile(Collected, 0.95, na.rm = T)
) %>%
knitr::kable(digits = 3)
```
Giải thích:
Tính tổng số dòng trong DataFrame df: Biến tot được gán giá trị là tổng số dòng của DataFrame df.
Nhóm dữ liệu theo cột POSEIDON: Sử dụng hàm group_by để nhóm dữ liệu trong df dựa trên giá trị của cột POSEIDON.
Tính toán và tóm tắt thông tin: Các hàm sau đây được sử dụng để tóm tắt thông tin cho từng nhóm (group):
n = n(): Tính số lượng dòng trong mỗi nhóm.
prop = 100*n/tot: Tính tỷ lệ phần trăm của số lượng dòng trong mỗi nhóm so với tổng số dòng tot.
mean = mean(Collected, na.rm = T): Tính trung bình của cột Collected trong mỗi nhóm, bỏ qua các giá trị NA (Not Available).
sd = sd(Collected, na.rm = T): Tính độ lệch chuẩn của cột Collected trong mỗi nhóm, bỏ qua các giá trị NA.
median = median(Collected, na.rm = T): Tính trung vị của cột Collected trong mỗi nhóm, bỏ qua các giá trị NA.
p5 = quantile(Collected, 0.05, na.rm = T): Tính phân vị thứ 5 của cột Collected trong mỗi nhóm, bỏ qua các giá trị NA.
p95 = quantile(Collected, 0.95, na.rm = T): Tính phân vị thứ 95 của cột Collected trong mỗi nhóm, bỏ qua các giá trị NA.
Hiển thị kết quả: Kết quả được hiển thị dưới dạng bảng với 3 chữ số sau dấu phẩy sử dụng knitr::kable().
Kết quả sẽ là một bảng tóm tắt các thông số đã tính cho từng nhóm của biến POSEIDON. Bảng này sẽ có các cột như n, prop, mean, sd, median, p5, p95, và sẽ được sắp xếp theo các giá trị unique của POSEIDON.
...
```{r,message = FALSE,warning=FALSE}
pals = c("#ffc403","#ff7403","#ff0357","#9203ff")
P1 = df %>% ggplot(aes(x = Collected,
y = POSEIDON,
fill= POSEIDON)) +
geom_density_ridges(stat = "binline",
scale = 2,
bins = 50,
alpha = 0.7,
draw_baseline = FALSE,
show.legend = F) +
labs(x="Number of retrieved oocytes", y = "POSEIDON groups") +
scale_fill_manual(values = pals,
name = "POSEIDON") +
scale_x_continuous(breaks = seq(0,80,5))+
coord_flip() +
theme_bw(9) +
geom_vline(xintercept = median(df$Collected),
linetype=2,
col="blue")+
theme(axis.text.x = element_text(angle = 45,
vjust = 1,
hjust=1)) +
theme(axis.text = element_text(size = 10,
color = "black"),
axis.title = element_text(size = 10,
color = "black"))
P1
```
Giải thích:
Đoạn mã R này sử dụng ggplot2 để vẽ một loại biểu đồ có tên là "density ridge plot", hiển thị phân phối của số lượng "retrieved oocytes" (số trứng thu được) trong các nhóm POSEIDON.
Dưới đây là giải thích chi tiết:
Cấu hình
Đặt màu sắc: pals = c("#ffc403","#ff7403","#ff0357","#9203ff") định nghĩa một vector chứa các mã màu HEX. Mỗi màu sẽ được sử dụng cho một nhóm POSEIDON trong đồ thị.
Cấu trúc cơ bản của ggplot
Khởi tạo ggplot: ggplot(aes(x = Collected, y = POSEIDON, fill= POSEIDON)) khởi tạo một đồ thị ggplot với trục x là Collected và trục y là POSEIDON. Màu sắc (fill) sẽ dựa vào giá trị của POSEIDON.
Geoms và Layers
Thêm lớp density ridges: geom_density_ridges tạo các "ridge" (dải) phân phối dựa trên dữ liệu.
stat = "binline": Sử dụng binning thay vì KDE (Kernel Density Estimation) để tính density. Kết quả sẽ tạo ra histogram
scale = 2: Quyết định mức độ "dày" của các dải.
bins = 50: Số lượng bins sử dụng trong binning.
alpha = 0.7: Độ trong suốt của các dải.
draw_baseline = FALSE: Không vẽ đường cơ sở ở dưới các dải.
show.legend = F: Không hiển thị legend.
Thêm labels và scales:
labs(): Đặt nhãn cho các trục.
scale_fill_manual(): Điều chỉnh màu sắc (fill) dựa trên vector màu pals.
scale_x_continuous(): Điều chỉnh các điểm ngắt (breaks) trên trục x.
Thêm các yếu tố đồ họa và điều chỉnh:
coord_flip(): Đảo trục x và y.
theme_bw(): Sử dụng chủ đề đen trắng.
geom_vline(): Thêm một đường thẳng dọc tại vị trí trung vị của Collected.
theme(): Điều chỉnh kích thước và màu sắc của văn bản trên các trục.
Hiển thị đồ thị
Cuối cùng hiển thị biểu đồ đã tạo.
Kết quả sẽ là một biểu đồ histogram hiển thị sự phân phối của số lượng "retrieved oocytes" trong các nhóm POSEIDON.
...
```{r}
pals_2 = c("#b5e617","#ffc403","#ff7403","#ff0357","#9203ff")
df2 = read.csv('POSEIDON_ART.csv',
sep = ';',
dec = ',',
fileEncoding = "UTF-8-BOM")
df2$POSEIDON = factor(df2$POSEIDON)
df2%>%gather(c(AFC,Age,AMH), key = 'Parameter', value = "Value")%>%
ggplot()+
geom_jitter(aes(x = Value,
y = Collected,
col = POSEIDON),
alpha = 0.2,
size = 1)+
geom_smooth(aes(x = Value,
y = Collected,
fill = POSEIDON,
col = POSEIDON),
alpha = 0.3,
show.legend = T,
method = 'glm')+
scale_color_manual(values = pals_2)+
scale_fill_manual(values = pals_2)+
labs(y="Number of retrieved oocytes", x = "Value") +
theme_bw(10)+
facet_wrap(~Parameter, scales = "free", ncol = 2)
```
Giải thích:
Đoạn mã R này sử dụng ggplot2 để vẽ biểu đồ tán xạ "jitter plot" kết hợp với đồ thị hàm hồi quy tuyến tính "smoothing line", để mô tả mối quan hệ giữa các tham số (AFC, Age, AMH) và số lượng noãn thu được (Collected) trong các nhóm POSEIDON
Đặt màu sắc: pals_2 = c("#b5e617","#ffc403","#ff7403","#ff0357","#9203ff") định nghĩa một vector chứa các mã màu HEX, sẽ được sử dụng cho mỗi nhóm POSEIDON.
Đọc dữ liệu: Đọc file CSV POSEIDON_ART.csv vào DataFrame df2 với các thiết lập như dấu phân cách ;, dấu thập phân , và mã hóa UTF-8-BOM.
Chuyển đổi dữ liệu: df2$POSEIDON = factor(df2$POSEIDON) chuyển đổi cột POSEIDON sang dạng factor.
Tiền xử lý dữ liệu
Chuyển đổi từ wide format sang long format: gather(c(AFC, Age, AMH), key = 'Parameter', value = "Value") chuyển đổi dữ liệu để có thể vẽ đồ thị facet.
Cấu trúc cơ bản của ggplot
Khởi tạo ggplot: ggplot() khởi tạo một đồ thị ggplot mà chưa có dữ liệu hay esthetics nào cả.
Geoms và Layers
Thêm jitter plot: geom_jitter tạo một biểu đồ jitter plot (tán xạ).
aes(x = Value, y = Collected, col = POSEIDON): Định nghĩa aesthetics cho jitter plot: tô màu điểm theo phân nhóm POSEIDON
alpha = 0.2, size = 1: Điều chỉnh độ trong suốt và kích thước của các điểm.
Thêm smoothing line: geom_smooth tạo một đường để mô tả mối quan hệ giữa Value và Collected.
method = 'glm': Sử dụng Generalized Linear Model để tạo đồ thị hàm tuyến tính.
alpha = 0.3, show.legend = T: Điều chỉnh độ trong suốt và hiển thị legend.
Điều chỉnh màu sắc và fill: scale_color_manual và scale_fill_manualđiều chỉnh màu sắc và fill cho các điểm và đồ thị.
Thêm labels: labs(y="Number of retrieved oocytes", x = "Value") đặt nhãn cho các trục.
Thêm các yếu tố đồ họa và điều chỉnh:
theme_bw(10): Sử dụng chủ đề đen trắng với cỡ chữ 10.
facet_wrap(~Parameter, scales = "free", ncol = 2): Tạo các "facet" (phần) dựa trên giá trị của cột Parameter, cho phép các trục có thang đo tự do và sắp xếp chúng thành 2 cột.
...
## Biện luận về tính phù hợp của quy luật phân phối NBI
...
```{r}
fit_poisson = fitdist(df$Collected, 'pois', optim.method = 'BFGS')
fit_negbin = fitdist(df$Collected, 'nbinom', optim.method = 'BFGS')
fit_norm = fitdist(df$Collected, 'norm', optim.method = 'BFGS')
P1 = denscomp(list(fit_poisson, fit_negbin, fit_norm),
legendtext = c("Poisson", "Negative binomial","Gaussian"),
plotstyle = 'ggplot',
fitcol = c('#03a9fc','#fc036f','#32a852'),
dempcol = c('grey'),
fitlty = 1,
fitlwd = 1)+
theme_bw()
P2 = cdfcomp(list(fit_poisson, fit_negbin,fit_norm),
legendtext = c("Poisson", "Negative binomial","Gaussian"),
plotstyle = 'ggplot',
fitcol = c('#03a9fc','#fc036f','#32a852'),
fitlty = 1,
datacol = 'grey')
P2 / P1
```
Giải thích:
Đoạn mã R này thực hiện các công việc sau:
Khớp dữ liệu thực nghiệm với các mô hình phân phối: Sử dụng hàm fitdist từ thư viện fitdistrplus để ước lượng các thông số cho ba loại phân phối: Phân phối Poisson (pois), phân phối nhị thức âm (nbinom), và phân phối chuẩn (norm) từ dữ liệu trong cột Collected của DataFrame df.
'optim.method = 'BFGS': chỉ định sử dụng phương pháp tối ưu hóa BFGS.
Vẽ đồ thị PDF (Probability Density Function): Hàm denscomp từ fitdistrplus được sử dụng để vẽ đồ thị mật độ xác suất (PDF) của các mô hình phân phối.
legendtext = c("Poisson", "Negative binomial","Gaussian"): Chú thích cho từng dòng.
plotstyle = 'ggplot': Sử dụng ggplot2 cho việc vẽ đồ thị.
fitcol: Màu sắc của các dòng.
dempcol = c('grey'): Màu của dữ liệu thực nghiệm.
fitlty = 1, fitlwd = 1: Loại và độ rộng của đường.
Vẽ đồ thị CDF (Cumulative Distribution Function): Tương tự như trên, nhưng sử dụng hàm cdfcomp để vẽ đồ thị tích lũy phân phối (CDF).
Kết hợp các biểu đồ: Cuối cùng, P2 / P1 sử dụng thư viện patchwork để kết hợp P1 (biểu đồ PDF) và P2 (biểu đồ CDF) thành một đồ thị duy nhất.
fit_poisson, fit_negbin, và fit_norm giúp ta hiểu liệu dữ liệu có tuân theo một trong những phân phối này không, và nếu có, thì phân phối nào là phù hợp nhất.
Các biểu đồ trực quan để so sánh mô hình phân phối đã khớp với dữ liệu thực nghiệm. Biểu đồ này giúp xác định phân phối nào là phù hợp nhất cho việc mô hình hóa cột Collected trong DataFrame df.
...
## Phân tích hồi quy bằng mô hình GLM NBI
...
```{r,message = T,warning=T}
mod_1 = gamlss(data = df,
formula = Collected ~ POSEIDON,
sigma.formula =~ POSEIDON,
family = NBI())
```
Giải thích:
Đoạn code này sử dụng thư viện gamlss trong R để xây dựng một mô hình phân phối nhị thức âm (Negative Binomial, viết tắt là NBI).
Xây dựng mô hình gamlss: Hàm gamlss được sử dụng để xây dựng một mô hình GLMM
data = df: Sử dụng DataFrame df để lấy dữ liệu.
formula = Collected ~ POSEIDON: Mô hình dự đoán biến kết quả Collected từ biến độc lập POSEIDON.
sigma.formula =~ POSEIDON: Mô hình sẽ cũng cố gắng ước lượng tham số (sigma) dựa trên POSEIDON.
family = NBI(): Sử dụng phân phối nhị thức âm.
Kết quả của mô hình mod_1 như sau:
```{r}
summary(mod_1)
```
...
```{r}
df%>%group_by(POSEIDON)%>%
summarize("Mean" = mean(Collected))%>%
knitr::kable(digits = 3)
```
...
```{r,message = FALSE,warning=FALSE}
abs = avg_comparisons(mod_1,
what = 'mu',
newdata = df,
variables = list(POSEIDON = 'pairwise'))
abs$adj_p = p.adjust(abs$p.value,
method = "bonferroni")
abs%>%dplyr::select(-c(2,6))%>%
knitr::kable(digits = 3)
```
Giải thích:
Đoạn code này sử dụng mô hình gamlss đã được xây dựng (mod_1) để thực hiện phân tích so sánh giá trị trung bình giữa các nhóm POSEIDON và sau đó điều chỉnh giá trị p.
avg_comparisons: Hàm này được sử dụng để so sánh các giá trị trung bình (mu) của biến phụ thuộc (Collected) giữa các nhóm POSEIDON.
what = 'mu': Thực hiện so sánh trên giá trị trung bình (mu) của phân phối.
newdata = df: Sử dụng DataFrame df để thực hiện so sánh.
variables = list(POSEIDON = 'pairwise'): So sánh sẽ được thực hiện giữa tất cả các cặp nhóm trong POSEIDON.
Điều chỉnh giá trị p bằng p.adjust: Giá trị p được điều chỉnh bằng cách sử dụng phương pháp Bonferroni.
method = "bonferroni": Sử dụng phương pháp điều chỉnh Bonferroni.
Sử dụng dplyr::select để loại bỏ cột thứ 2 và thứ 6 từ DataFrame. Tiếp theo, knitr::kable được sử dụng để hiển thị kết quả dưới dạng bảng với 3 chữ số sau dấu phẩy.
...
```{r,results='hide'}
augment(mod_1)%>%mutate(UL= .fitted +.se.fit*1.96,
LL=.fitted-.se.fit*1.96)%>%
group_by(POSEIDON)%>%
summarize_at('.fitted', median) -> m1_sum
augment(mod_1)%>%mutate(UL= .fitted +.se.fit*1.96,
LL=.fitted-.se.fit*1.96)%>%
ggplot()+
geom_jitter(aes(y = POSEIDON,
x = Collected,
col = POSEIDON),
alpha = 0.2,
width = 0.01,
height = 0.1,
show.legend = F)+
geom_density_ridges(aes(y = POSEIDON,
x = Collected,
fill = POSEIDON),
stat = "binline",
scale = 1,
bins = 50,
alpha = 0.5,
draw_baseline = FALSE,
show.legend = F)+
geom_errorbar(aes(y=POSEIDON,
xmin=exp(LL),
xmax=exp(UL)),
width=0.2,
size=1) +
geom_path(data = m1_sum,
aes(y=POSEIDON,
x=exp(.fitted),
group = 1))+
geom_point(aes(y=POSEIDON,
x=exp(.fitted)),
size=3)+
labs(x="Number of retrieved oocytes", y = "POSEIDON groups") +
scale_fill_manual(values = pals) +
scale_color_manual(values = pals) +
scale_x_continuous(breaks = seq(0,60,5))+
coord_flip()+
theme_bw()
```
Giải thích:
Đoạn mã R này sử dụng mô hình gamlss (mod_1) để tạo ra dữ liệu về giá trị ước lượng của số lượng noãn thu được (Collected), sau đó vẽ một một biểu đồ so sánh giữa các nhóm POSEIDON
Các bước chi tiết:
Thêm cột mới cho DataFrame: Sử dụng augment(mod_1) để thêm các cột được tính toán từ mô hình vào dataframe. Cụ thể là .fitted (giá trị ước tính) và .se.fit (sai số chuẩn của giá trị ước tính).
UL= .fitted +.se.fit*1.96 và LL=.fitted-.se.fit*1.96 tạo ra các giới hạn trên (Upper Limit) và dưới (Lower Limit) của khoảng tin cậy 95%.
Tính giá trị trung bình cho mỗi nhóm: group_by(POSEIDON) và summarize_at('.fitted', median) được sử dụng để tính giá trị trung bình của .fitted cho mỗi nhóm POSEIDON. Kết quả được lưu vào biến m1_sum.
Tạo biểu đồ:
geom_jitter(): Vẽ các điểm dữ liệu với độ nhiễu nhỏ để tránh hiện tượng chồng lắp.
geom_density_ridges(): Vẽ các đường density cho mỗi nhóm.
geom_errorbar(): Vẽ dải tin cậy 95% dựa trên các giới hạn UL và LL đã tính.
geom_path() và geom_point(): Vẽ đường và điểm cho giá trị trung bình của .fitted từ m1_sum
...
```{r}
mod_2 = gamlss(data = df2,
formula = Collected ~ Age + AMH + AFC,
sigma.formula =~ Age + AMH + AFC,
family = NBI())
summary(mod_2)
```
Giải thích:
Đoạn code này sử dụng hàm gamlss từ thư viện gamlss để xây dựng một mô hình ước lượng biến kết quả là số noãn thu được (Collected) trong dataFrame df2. Mô hình dựa vào quy luật phân phối Negative Binomial Inflation (NBI).
Các thành phần của hàm gamlss:
data = df2: Cho biết rằng dữ liệu cho mô hình sẽ được lấy từ dataframe df2.
formula = Collected ~ Age + AMH + AFC: Đây là công thức của mô hình,với biến kết quả là Collected và các biến độc lập là Age, AMH, và AFC.
sigma.formula =~ Age + AMH + AFC: Điều này chỉ định mô hình ước lượng cả tham số sigma (có ý nghĩa phương sai) như một hàm từ các biến Age, AMH, và AFC.
Bằng cách mô hình hóa cả trung bình (mu) và phương sai (sigma), mô hình có thể bắt được cả tính trung bình và sự biến động của biến mục tiêu, tùy thuộc vào các biến độc lập.
family = NBI(): Chỉ định rằng phân phối của biến mục tiêu Collected sẽ được mô hình hóa dựa trên phân phối Negative Binomial Type I. Phân phối này thích hợp cho kết quả số đếm với tính chất quá phân tán.
...
```{r}
abs = avg_comparisons(mod_2,
what = "mu",
newdata = df2)
abs%>%dplyr::select(-c(1,6))%>%
knitr::kable(digits = 3)
```
Giải thích:
Đoạn mã này sử dụng hàm avg_comparisons để khảo sát marginal effect trên mô hình mod_2
mod_2: Đây là mô hình mà bạn muốn thực hiện phân tích
what = "mu": Bạn chỉ định rằng bạn muốn so sánh trung bình (mu) giữa các nhóm.
Phân tích này cho phép so sánh bắt cặp tuần tự giá trị trung bình giữa các nhóm khác nhau, hay so sánh hiệu ứng/ảnh hưởng của mỗi bậc giá trị biến độc lập đối với biến kết quả.
newdata = df2: Đây là dữ liệu mới mà bạn muốn sử dụng để dự đoán từ mô hình. Trong trường hợp này, nó cùng là dữ liệu đã được sử dụng để xây dựng mô hình (df2).
dplyr::select(-c(1,6)): Đoạn mã này sử dụng hàm select từ thư viện dplyr để loại bỏ các cột 1 và 6 từ DataFrame abs.
knitr::kable(digits = 3): Đây là hàm từ thư viện knitr dùng để tạo ra một bảng kết quả, với mỗi số được làm tròn đến 3 chữ số thập phân.
...
```{r}
preds = predictions(model = mod_2,
what = "mu",
newdata = datagrid(Age = seq(20,40,5),
grid_type = "counterfactual"))
preds%>%ggplot(aes(x = Age,
y = estimate)) +
stat_lineribbon(aes(y = estimate),
fill = "gold",
.width = c(.95, .75, .50),
alpha = 1/5,
show.legend = F) +
stat_dotsinterval(aes(fill = factor(Age)),
quantiles = 50,
alpha = 0.7,
.width = c(0.75, 0.95),
point_interval = 'median_qi',
show.legend = F)+
labs(y="Number of retrieved eggs", x = "Age") +
scale_y_continuous(limits = c(0,40), breaks = seq(0,40,5))+
scale_x_continuous(limits = c(20,45), breaks = seq(20,40,5))+
theme_bw(10)+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 10, color = "black"))
```
Giải thích:
Đoạn code này nhằm tạo một biểu đồ khảo sát mối liên hệ giữa giá trị ước tính và khoảng tin cậy về "Số lượng noãn thu thập được" ('Number of retrieved eggs') theo "Tuổi" ('Age').
predictions(model = mod_2, what = "mu", newdata = datagrid(...)): Hàm này tạo ra giá trị dự đoán từ mô hình mod_2 (mô hình đã được tạo trước đó) dựa trên dữ liệu mới được tạo ra từ datagrid, với 'Age' từ 20 đến 40, tăng dần 5 đơn vị.
ggplot() + ...: Đây là các đoạn mã để tạo biểu đồ.
stat_lineribbon: Tạo ra một dải băng (ribbon) thể hiện khoảng tin cậy của dự đoán. Các dải với độ rộng khác nhau (.95, .75, .50) cho thấy các khoảng tin cậy khác nhau.
stat_dotsinterval: Thêm các điểm vào biểu đồ để thể hiện trung bình dự đoán và khoảng tin cậy của chúng.
labs, scale_x_continuous, scale_y_continuous: Đặt tên cho các trục và điều chỉnh phạm vi của chúng.
theme_bw, theme: Điều chỉnh chủ đề và kiểu chữ của biểu đồ.
Biểu đồ này sẽ cho ta thấy hình ảnh về giá trị dự đoán về số lượng trứng thu được tùy thuộc vào độ tuổi và khoảng tin cậy cho những dự đoán này
...
```{r}
preds = predictions(model = mod_2,
what = "mu",
newdata = datagrid(AMH = seq(0,10,1),
grid_type = "counterfactual"))
preds%>%ggplot(aes(x = AMH,
y = estimate)) +
stat_lineribbon(aes(y = estimate),
fill = "gold",
.width = c(.95, .75, .50),
alpha = 1/5,
show.legend = F) +
stat_dotsinterval(aes(fill = factor(AMH)),
quantiles = 50,
alpha = 0.5,
.width = c(0.75, 0.95),
point_interval = 'median_qi',
show.legend = F)+
labs(y="Number of retrieved eggs", x = "AMH") +
scale_y_continuous(limits = c(0,35), breaks = seq(0,35,5))+
scale_x_continuous(limits = c(0,11), breaks = seq(0,10,1))+
theme_bw(10)+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 10, color = "black"))
```
Giải thích: cấu trúc và công dụng đoạn code này tương tự như trên, nhưng lần này thay vì Tuổi ta quan tâm đến biến AMH. Việc dự đoán được thực hiện cho các giá trị AMH từ 0 đến 10, cho phép ước lượng về sự thay đổi của số noãn thu được tùy theo mức AMH thấp hoặc cao.
...
```{r}
preds = predictions(model = mod_2,
what = "mu",
newdata = datagrid(AFC = seq(1,30,2.5),
grid_type = "counterfactual"))
preds%>%ggplot(aes(x = AFC,
y = estimate)) +
stat_lineribbon(aes(y = estimate),
fill = "gold",
.width = c(.95, .75, .50),
alpha = 1/5,
show.legend = F) +
stat_dotsinterval(aes(fill = factor(AFC)),
quantiles = 50,
alpha = 0.7,
.width = c(0.75, 0.95),
point_interval = 'median_qi',
show.legend = F)+
labs(y="Number of retrieved eggs", x = "AFC") +
scale_y_continuous(limits = c(0,40), breaks = seq(0,40,5))+
scale_x_continuous(limits = c(0,31), breaks = seq(0,30,5))+
theme_bw(10)+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 10, color = "black"))
```
Đoạn code vẽ biểu đồ này tương tự như trên, nhưng cho biến AFC
...
```{r}
preds = predictions(model = mod_2,
what = "mu",
newdata = datagrid(AFC = seq(1,30,2.5),
Age = c(25,30,35,40),
grid_type = "counterfactual"))
preds%>%ggplot(aes(x = AFC,
y = estimate)) +
stat_lineribbon(aes(y = estimate),
fill = "gold",
.width = c(.95, .75, .50),
alpha = 1/5,
size = 1,
show.legend = F) +
stat_dotsinterval(aes(fill = factor(AFC)),
quantiles = 50,
alpha = 0.7,
size = 1,
.width = c(0.75, 0.95),
point_interval = 'median_qi',
show.legend = F)+
labs(y="Number of retrieved eggs", x = "AFC") +
scale_y_continuous(limits = c(0,40), breaks = seq(0,40,5))+
scale_x_continuous(limits = c(0,31), breaks = seq(0,30,5))+
theme_bw(10)+
facet_wrap(~Age, scales = "free")+
theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
theme(axis.text = element_text(size = 10, color = "black"),
axis.title = element_text(size = 10, color = "black"))
```
Giải thích:
Đoạn R code này tạo ra một biểu đồ trình bày sự thay đổi của giá trị dự đoán và khoảng tin cậy của số lượng noãn thu được ('Number of retrieved eggs') dựa trên sự tổ hợp giữa yếu tố "AFC" (Antral Follicle Count) và "Age" (Tuổi). Đây là một phân tích đa biến: Vì AFC và Tuổi đều được sử dụng, biểu đồ này giúp chúng ta hiểu rõ hơn về mức độ các yếu tố này có thể ảnh hưởng đến số lượng trứng thu thập được.
predictions(model = mod_2, what = "mu", newdata = datagrid(...)): Kết quả dự đoán từ mô hình mod_2 dựa trên dữ liệu mới. Đây là một dữ liệu có tính "phản thực tế" (counterfactual) với giả định AFC có giá trị từ 1 đến 30 (tăng dần 2.5 đơn vị) tổ hợp với các độ Tuổi [25, 30, 35, 40].
facet_wrap(~Age, scales = "free"): Tạo các "facet" (biểu đồ con) dựa trên các giá trị Tuổi khác nhau. Mỗi "facet" có thể có tỷ lệ trục y riêng. Cách trình bày này có ý nghĩa phân tích hiệu ứng của sự tổ hợp giữa tuổi và AFC đối với số lượng noãn thu được.
## Diễn đạt kết quả phân tích
...
## Thông điệp rút gọn làm hành trang
...