Mô hình hồi quy Gamma

1 Giới thiệu

So sánh giá trị một đại lượng giữa nhiều phân nhóm độc lập là một vấn đề thường gặp trong nghiên cứu y học lâm sàng. Giải pháp thường quy là phân tích phương sai một yếu tố (đơn biến) (One-way ANOVA). Tuy nhiên, ANOVA cổ điển lại đặt ra giả định về tính phân phối chuẩn của dữ liệu và giả định này gần như không phù hợp với thực tiễn, bởi vì hầu hết những đại lượng sinh lý bệnh học trong tự nhiên đều có phân phối lệch dương.

Ở đây chúng tôi giới thiệu một giải pháp linh hoạt và phổ quát hơn cho bài toán so sánh giá trị trung bình của biến định lượng giữa nhiều phân nhóm độc lập, đó là mô hình tuyến tính tỗng quát (Generalized linear model, GLM) với phân phối Gamma.

2 Công cụ cần thiết:

Quy trình phân tích cần những công cụ sau đây:

• Hệ sinh thái tidyverse: với thư viện ggplot2 cho đồ họa thống kê và dplyr để thao tác dữ liệu

• Thư viện ggridges : để vẽ biểu đổ mật độ phân phối riêng cho từng phân nhóm (ggplot2 cơ bản không hỗ trợ loại biểu đồ này).

• Thư viện fitdistrplus cho phép kiểm tra tính phù hợp giữa dữ liệu thực nghiệm và một quy luật phân phối xác suất lý thuyết, như Gaussian hoặc Gamma.

• Thư viện gamlss: cho phép khớp các mô hình hồi quy cho nhiều họ phân phối khác nhau, bao gồm Gamma.

• Thư viện broom: cho phép thao tác trên kết quả mô hình hồi quy gamlss.

• Thư viện marginaleffects: dùng cho phân tích hậu kiểm (post-hoc analysis) trên kết quả mô hình hồi quy gamlss

• Dữ liệu đầu vào có cấu trúc gồm 1 biến kết quả Y thuộc loại định lượng liên tục (kiểu dữ liệu double, numeric, int) và 1 biến phân nhóm X thuộc loại định tính rời rạc (định dạng character hoặc factor). Lưu ý, phân phối Gamma có miền xác định là tập số thực dương R+, nên trong dữ liệu Y không được có giá trị < 0 hoặc = 0.

# Thao tác dữ liệu và đồ họa
library(tidyverse)
library(ggridges)

# Kiểm tra quy luật phân phối
library(fitdistrplus)

# Mô hình GLM
library(broom)
library(gamlss)

library(marginaleffects)

Thí dụ minh họa sử dụng dataset “Five_protocols.csv” từ một nghiên cứu có thực trên 1005 phụ nữ được thực hiện thụ tinh nhân tạo. Mục tiêu đặt ra là khảo sát giá trị hormone LH ở ngày trigger (biến LH) giữa 5 loại phác đồ kích thích buồng trứng (biến Protocol), gồm:

• Mild = mild stimulation sử dụng Letrozole,
• GnRH_ant = Phác đồ đối vận GnRH,
• Fol_long = Follicular long acting,
• PPOS: Progestin-Primed Ovarian Stimulation,
• Lut_short = Lutean phase short acting.

3 Chuẩn bị dữ liệu

Dữ liệu được lưu ở định dạng .csv với dấu ngăn cách là “;” và dấu thập phân là “.”

1. Ta dùng hàm read.csv để tải dữ liệu vào R và lưu trong dataframe df.

2. Dùng hàm na.omit() để loại bỏ những đơn vị quan sát bị thiếu sót dữ liệu (missing value)

3. Chuyển Protocol từ 1 biến kiểu kí tự (character) thành factor, để tương thích với quy trình phân tích.

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

df$Protocol = as.factor(df$Protocol)

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

Protocol	Age	BMI	LH	P
Lut_short	40	23.00	3.26	0.05
Lut_short	30	22.70	1.97	0.69
GnRH_ant	30	24.20	1.55	0.19
Fol_long	36	20.20	0.69	0.56
Lut_short	36	20.80	2.95	0.38
PPOS	39	20.05	1.87	0.60

Một cách mặc định, biến rời rạc loại factor sẽ được phân cấp theo thứ tự alphabet, mỗi bậc giá trị sẽ tương ứng với số thứ tự. Trong trường hợp này, thứ tự các loại protocol là:

1 = Fol_long (sẽ được xem là nhóm tham chiếu khi phân tích hồi quy)

2 = GnRH_ant

3 = Lut_short

4 = Mild

5 = PPOS

levels(df$Protocol)

## [1] "Fol_long"  "GnRH_ant"  "Lut_short" "Mild"      "PPOS"

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

Một quy trình phân tích gồm 4 bước sẽ được tiến hành nhằm giải đáp câu hỏi nghiên cứu đã đặt ra:

1. Đặc tính phân bố của giá trị LH ở 5 nhóm phác đồ sẽ được mô tả bằng trung bình, độ lệch chuẩn, trung vị và bách phân vị 5-95, đồng thời trình bày trực quan bằng biểu đồ mật độ phân bố (Kernel density plot) hoặc Boxplot, violin plot.

2. Dữ liệu LH sẽ được khớp lần lượt với 2 quy luật phân phối Gamma và Gaussian, sau đó đối chiếu đồ thị hàm PDF và CDF để xác nhận tính phù hợp tốt hơn của phân phối Gamma so với Gaussian.

3. Mô hình hồi quy tuyến tính tổng quát với phân phối Gamma sẽ được áp dụng để khảo sát giá trị trung bình LH giữa 5 phác đồ.

4. Hậu kiểm so sánh bắt cặp tuần tự nhằm kiểm tra một loạt giả định về sự khác biệt của giá trị LH giữa các nhóm phác đồ.

Suy diễn thống kê dựa trên kiểm định giả thuyết vô hiệu với ngưỡng ý nghĩa thống kê p = 0.05. Hiệu chỉnh Bonferroni được áp dụng khi kiểm tra hàng loạt giả định.

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

Mục tiêu của công đoạn này nhằm thăm dò và mô tả sơ khởi về đặc tính phân phối của LH giữa các nhóm phác đồ khác nhau. Nhận định rút ra từ phân tích mô tả này cho phép đặt ra những giả thuyết cần chứng minh và đề xuất hướng phân tích phù hợp.

1. Dựng bảng thống kê mô tả

• Sử dụng hàm group_by của thư viện dplyr để phân chia dữ liệu thành 5 phần, tương ứng 5 bậc giá trị của biến Protocol

• Áp dụng hàm summarize của thư viện dplyr để tính mean, sd, median, p5, p95, min và max của LH cho mỗi nhóm nhỏ

• Dùng hàm arrange của dplyr để Xếp thứ tự theo giá trị trung bình LH tăng dần

Kết quả như sau:

df%>% 
  group_by(Protocol)%>%
  summarize(n = n(),
            mean = sprintf("%0.2f", mean(LH)),
            sd = sprintf("%0.2f",sd(LH)),
            median = sprintf("%0.2f", median(LH)),
            p5 = sprintf("%0.2f", quantile(LH, 0.05)),
            p95 = sprintf("%0.2f", quantile(LH, 0.95)),
            min= sprintf("%0.2f", min(LH)),
            max = sprintf("%0.2f", max(LH)),
            )%>%
    arrange(mean)%>%
    knitr::kable()

Protocol	n	mean	sd	median	p5	p95	min	max
Fol_long	261	1.37	1.64	0.87	0.13	4.18	0.10	12.72
Mild	75	12.02	8.45	8.70	2.52	27.65	1.81	36.90
Lut_short	220	2.50	1.65	2.12	0.75	4.92	0.17	14.51
PPOS	204	3.65	3.90	2.41	0.65	10.68	0.09	29.15
GnRH_ant	245	7.78	6.01	6.05	1.55	19.06	0.41	45.02

Theo kết quả này, có sự khác biệt về giá trị LH trung bình giữa 5 phác đồ. Một cách cụ thể, LH trung bình có khuynh hướng tăng dần theo thứ tự:

Fol_long < Short < Mild < < PPOS < GnRH_ant < Mild

2. So sánh trực quan bằng biểu đồ

Ta sẽ sử dụng 2 loại biểu đồ cho phép so sánh đặc điểm phân bố của LH giữa 5 phác đồ, đó là Biểu đồ mật độ phân phối 1 chiều (Density plot) hoặc biểu đồ Boxplot

Đầu tiên, ta tạo dataframe med_df trình bày giá trị trung bình của LH giữa 5 phân nhóm, sử dụng kết hợp 2 hàm group_by và summarize_at. Dataframe này sẽ được dùng để vẽ biểu đồ tuyến kí và định vị trung bình LH cho mỗi nhóm.

Tiếp theo, ta lần lượt dùng các hàm:

• geom_density_ridges (data = df) để vẽ 5 biểu đồ density plot cho LH ở 5 nhóm, màu nền theo nhóm phác đồ. Lưu ý rằng hàm geom_density_ridges mặc định trục x = biến định lượng, y = phân nhóm.

• geom_path cho med_df để vẽ biểu đồ tuyến kí nối các điểm giá trị LH trung bình;

• geom_point cho med_df để chấm 5 điểm giá trị LH trung bình

Lưu ý rằng trục X luôn là LH, y luôn là Protocol cho 3 hàm trên

Sau đó, dùng hàm coord_flip() để xoay trục 90°, như vậy trục x và y sẽ hoán chuyển cho nhau, hình ảnh sẽ biểu diễn theo chiều dọc thay vì chiều ngang.

Hàm scale_x_continuous cho phép tùy chỉnh thang đo trục X, ở đây giới hạn limits từ -0.5 đến 40, mốc thang đo 5 đon vị

med_df = df %>% 
  group_by(Protocol)%>%
  summarize_at('LH', mean)

ggplot()+
  geom_density_ridges(data = df,
                      aes(y = Protocol, 
                          x = LH,
                          fill = Protocol),
                      scale = 0.8,
                      alpha = 0.5,
                      show.legend = F)+
  geom_path(data = med_df,
            aes(y = Protocol, x = LH, group = 1),
            color = 'black',
            size = 1,
            alpha = 0.5)+
  geom_point(data = med_df,
            aes(y = Protocol, x = LH),
            color = 'black',
            size = 2.5,
            alpha = 0.9)+
  coord_flip()+
  scale_x_continuous(limits = c(-0.5,40), 
                     breaks = seq(0,40,by = 5))+
  labs(x = 'LH level (mIU/mL)', y = 'Protocols')+
  theme_bw(10)

Kết quả thu được là 5 biểu đồ mật độ phân bố như trong hình. Các biểu đồ này có bản chất là đồ thị hàm probability desity (PDF) của biến LH, trục Y tương ứng với thang đo của LH). Chúng mô tả đặc tính phân phối của LH trong 5 nhóm phác đồ (trục X). Các chấm màu đen tương ứng với giá trị trung bình LH ở mỗi nhóm, được nối với nhau bằng các đoạn thẳng thể hiện khuynh hướng tăng/giảm của giá trị này.

Dựa vào hình ảnh này, ta có thể rút ra được hai thông tin: thứ nhất, có sự khác biệt rõ nét, không chỉ về giá trị trung bình nhưng còn về đặc tính phân phối của LH, thí dụ độ phân tán, độ lệch… giữa 5 nhóm. Thứ hai là ở mỗi phân nhóm, LH chắc chắn không thể phù hợp với quy luật phân phối chuẩn (Gaussian). Phân phối của LH không đối xứng, thậm chí có độ lệch rất cao ở các nhóm GnRh_ant, Mild và PPOS.

Ta có thể làm tương tự để tạo ra biểu đổ boxplot, chỉ cần thay hàm geom_density_ridges bằng hàm geom_boxplot. Lưu ý rằng với dạng biểu đồ boxplot mặc định thì trục x chỉ Protocol, còn trục y cho LH, ta vẫn dùng hàm coord_flip() để xoay trục, trình bày theo phương ngang, vì khả năng so sánh bằng thị giác sẽ tốt hơn theo chiều ngang so với chiều dọc.

ggplot()+
  geom_boxplot(data = df,
                      aes(x = Protocol, 
                          y = LH,
                          fill = Protocol),
                      alpha = 0.5,
                      show.legend = F)+
  geom_point(data = med_df,
            aes(x = Protocol, y = LH),
            color = 'black',
            shape = 18,
            size = 5,
            alpha = 0.9)+
  scale_y_continuous(limits = c(-0.5,40), 
                     breaks = seq(0,40,by = 5))+
  labs(y = 'LH level (mIU/mL)', x = 'Protocols')+
  coord_flip()+
  theme_bw(10)

Biểu đồ boxplot trình bày các trị số thống kê mô tả đại diện gồm trung vị, tứ phân vị (phân vị thứ 25, 50, 75), các chấm tròn nhỏ là những giá trị vượt ra khỏi khoảng tứ phân vị. Chấm màu đen lớn nhất là vị trí của trung bình. Tuy không mô tả hình ảnh của toàn thể mật độ phân phối, biểu đồ Boxplot vừa đủ để cho ta thấy cùng những thông tin về sự tương phản rõ nét giữa 5 phân nhóm, ta cũng cảm nhận được LH không có phân phối chuẩn (trung vị không trùng với trung bình, 2 đầu boxplot không cân xứng, những giá trị cao lệch hẳn về chiều dương của thang đo LH…).

6 Bước 2: Biện luận về tính phù hợp của phân phối Gamma

Mục tiêu của công đoạn này nhằm chứng minh rằng để ước lượng giá trị LH thì phân phối Gamma phù hợp hơn so với phân phối Gaussian.

Về mặt lý thuyết, ta hoàn toàn có thể chọn quy luật phân phối Gamma cho biến LH vì những lý do sau: Thang đo LH có giá trị >0, phân phối của LH có hình ảnh bất đối xứng và lệch dương.

Tuy nhiên, ta có thể dựa vào chứng cứ xác thực hơn bằng quy trình sau:

1. Đầu tiên, ta dùng hàm fitfist để lần lượt khớp dữ liệu LH với 2 quy luật phân phối Gaussian và Gamma,

2. sau đó đối chiếu hình ảnh đồ thị hàm PDF và CDF so với dữ liệu thực nghiệm, bằng 2 hàm denscomp và cdfcomp:

fit_gam = fitdist(df$LH, 'gamma', optim.method = 'BFGS')
fit_gauss = fitdist(df$LH, 'norm', optim.method = 'BFGS')

denscomp(list(fit_gauss, fit_gam), 
         legendtext = c("Gaussian", "Gamma"),
         plotstyle = 'ggplot',
         fitcol = c('#03a9fc','#fc036f'),
         dempcol = c('grey'),
         fitlty = 1,
         fitlwd = 1)+
  labs(x = 'LH on trigger day (mUI/mL)')+ 
  theme_bw()

Chú thích: Biểu đồ trên cho phép đối chiếu đồ thị hàm PDF (mật độ xác suất) của 2 phân phối Gaussian (màu xanh), Gamma (màu đỏ) với hình ảnh biểu đồ tần suất trên thực nghiệm (màu xám).

cdfcomp(list(fit_gauss, fit_gam), 
        legendtext = c("Gaussian", "Gamma"),
        plotstyle = 'ggplot',
        fitcol = c('#03a9fc','#fc036f'),
        fitlty = 1,
        datacol = 'grey')+
  labs(x = 'LH on trigger day (mUI/mL)')+ 
  theme_bw()

Chú thích: Biểu đồ trên cho phép đối chiếu đồ thị hàm CDF (mật độ xác suất tích lũy) của 2 phân phối Gaussian (màu xanh), Gamma (màu đỏ) với hàm CDF trên thực nghiệm (màu xám).

Cả 2 kết quả này cho thấy phân phối Gamma thích hợp hơn so với phân phối Gaussian để ước lượng giá trị LH, vì hàm PDF và CDF của phân phối Gamma phù hợp hơn với dữ liệu thực nghiệm.

Lợi ích của việc chọn quy luật phân phối đúng cho biến kết quả mô hình hồi quy, đó là nó cho phép ước lượng và suy luận thống kê chính xác hơn, phù hợp thực tế hơn.

7 Bước 3: Dựng mô hình hồi quy với phân phối Gamma

Như đã trình bày trong kế hoạch phân tích, công cụ thống kê chính mà ta sử dụng là một mô hình hồi quy tuyến tính tổng quát (GLM). Hệ thống mô hình GLM do 2 nhà thống kê John Nelder và Robert Wedderburn thiết lập vào năm 1972 cho phép ước lượng giá trị của biến kết quả theo nhiều quy luật phân phối đa dạng.

Một cách tổng quát, mô hình GLM ước lượng giá trị trung bình (tham số \(\mu\)) của biến ngẫu nhiên Y, điều kiện hóa theo một tập biến X thông qua một hàm liên kết \(g\) và một mô hình hồi tuyến tính \(X\beta\).

\[E(Y|X) = \mu_{Y} \sim g^{-1}(X\beta)\] Trong bài toán hiện thời, biến kết quả Y là giá trị của LH ngày trigger, được giả định là một biến ngẫu nhiên theo quy luật phân phối Gamma với 2 tham số \(\mu\) và \(\sigma\)

Hàm mật độ phân phối (PDF) của phân phối Gamma được xác định như sau:

\[f(y|\mu ,\sigma)=\frac{1}{\left ( \sigma ^{2} \mu \right )^{1/\sigma^{2}}}\frac{y^{\frac{1}{\sigma ^{2}}-1}e^{-y/(\sigma ^{2}\mu )}}{\Gamma \left ( 1/\sigma^{2} \right )}\]

\(\mu\) có ý nghĩa như trọng tâm (central tendency) hay giá trị trung bình của biến kết quả cần ước lượng, còn giá trị của \(\sigma\) tương ứng với mức độ phân tán và xác định hình ảnh của hàm mật độ phân phối. Một khi có mu và sigma, có thể ước tính độ lệch chuẩn (và phương sai) của LH theo công thức:

\[SD = \sigma * \mu\]

Phân tích hồi quy Gamma dùng thư viện gamlss. Đầu tiên, ta khớp một mô hình hồi quy Gamma chỉ chứa Intercept. Mô hình này tương ứng với việc ước lượng giá trị trung bình của LH cho toàn bộ mẫu khảo sát, không phân biệt nhóm phác đồ nào.

Ta dùng hàm gamlss, với công thức formula = LH ~ 1, và family = GA

# Mô hình chỉ có intercept

m0 = gamlss(formula = LH ~ 1,
             data=df,
             family = GA,
             trace=F,
             parallel="multicore",
             ncpus = nC)

Ta sẽ dùng mô hình cơ bản này để làm tham chiếu và so sánh với mô hình chính trong phân tích, nhằm chứng minh rằng yếu tố chia nhóm theo phác đồ thực sự có ảnh hưởng đến giá trị trung bình LH.

Nội dung mô hình cơ bản (chỉ chứa Intercept, ước lượng LH trung bình cho toàn quần thể) có nội dung như sau:

summary(m0)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = LH ~ 1, family = GA, data = df, trace = F,  
##     parallel = "multicore", ncpus = nC) 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.49042    0.03189   46.73   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.01101    0.01960   0.561    0.575
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1005 
## Degrees of Freedom for the fit:  2
##       Residual Deg. of Freedom:  1003 
##                       at cycle:  2 
##  
## Global Deviance:     5005.425 
##             AIC:     5009.425 
##             SBC:     5019.251 
## ******************************************************************

gamlss dựng 2 mô hình hồi quy riêng cho \(\mu\) và \(\sigma\), cả 2 mô hình đều sử dụng hàm liên kết là logarit, vì vậy thang đo của các hệ số hồi quy hiện thời là log(LH), để đưa về thang đo nguyên thủy (mUI/mL) ta phải dùng hàm exponential.

Mô hình thứ nhất cho biết giá trị trung bình của LH cho toàn quần thể \(\mu = exp(1.49042) = 4.438959\)

Ta có thể đối chiếu với kết quả ước lượng từ dữ liệu thực nghiệm, và thấy rằng Intercept của mô hình m0 cơ bản rất phù hợp với giá trị trung bình của LH

mean(df$LH)

## [1] 4.438955

2. Tiếp theo, ta dựng một mô hình khác với 1 biến độc lập là phân nhóm Protocol. Mô hình này cho phép ước lượng giá trị trung bình LH ở 5 phân nhóm Protocol khác nhau, hay nói cách khác, khảo sát liên hệ giữa biến protocol và giá trị trung bình LH (hay hiệu ứng của Protocol đối với giá trị trung bình của LH).

\[\mu_{LH} \sim \beta_0*Long + \beta_1*GnRHant + \beta_2*Short + \beta_3*Mild + \beta_4*PPOS\]

Mô hình này có công thức cho tham số chính (trung bình của phân phối Gamma) là LH ~ Protocol; và cho tham số sigma: sigma.formula = ~ Protocol.

# Mô hình ANOVA 1 biến

m1 = gamlss(formula = LH ~ Protocol,
            sigma.formula = ~ Protocol,
            data=df,
            family = GA,
            trace=F,
            parallel="multicore",
            ncpus = nC)

Ở đây ta ước lượng đồng thời cả 2 tham số của phân phối Gamma là mu (trung bình) và sigma (tham số kiểu hình) theo biến Protocol. Vì ta không chỉ quan tâm đến việc Protocol sẽ làm thay đổi giá trị trung bình LH (trọng tâm của phân phối Gamma, tham số \(\mu\)), mà còn muốn biết liệu phân loại Protocol có tác động gây thay đổi hình ảnh (độ phân tán, độ lệch) của phân phối Gamma (tham số \(\sigma\)) hay không ?

Nội dung của mô hình chính (ước lượng LH theo phác đồ) có nội dung như sau

summary(m1)

## ******************************************************************
## Family:  c("GA", "Gamma") 
## 
## Call:  gamlss(formula = LH ~ Protocol, sigma.formula = ~Protocol,  
##     family = GA, data = df, trace = F, parallel = "multicore",      ncpus = nC) 
## 
## 
## Fitting method: RS() 
## 
## ------------------------------------------------------------------
## Mu link function:  log
## Mu Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.31235    0.05756   5.426 7.23e-08 ***
## ProtocolGnRH_ant   1.73978    0.07409  23.481  < 2e-16 ***
## ProtocolLut_short  0.60523    0.06878   8.800  < 2e-16 ***
## ProtocolMild       2.17429    0.09639  22.556  < 2e-16 ***
## ProtocolPPOS       0.98292    0.08235  11.935  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## Sigma link function:  log
## Sigma Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)       -0.07264    0.03903  -1.861 0.063003 .  
## ProtocolGnRH_ant  -0.24182    0.05717  -4.230 2.56e-05 ***
## ProtocolLut_short -0.51025    0.05988  -8.521  < 2e-16 ***
## ProtocolMild      -0.32842    0.08573  -3.831 0.000136 ***
## ProtocolPPOS      -0.10027    0.05946  -1.686 0.092072 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## ------------------------------------------------------------------
## No. of observations in the fit:  1005 
## Degrees of Freedom for the fit:  10
##       Residual Deg. of Freedom:  995 
##                       at cycle:  2 
##  
## Global Deviance:     4275.111 
##             AIC:     4295.111 
##             SBC:     4344.238 
## ******************************************************************

Mô hình m1 này cho phép ước lượng giá trị trung bình của LH ở mỗi phân nhóm Protocol. Lưu ý rằng phác đồ Fol_long được chọn làm tham chiếu, vì vậy Intercept trong mô hình tương ứng với trung bình của log(LH) trong phân nhóm Protocol = Fol_long.

Ta có thể đối chiếu với kết quả thống kê mô tả ban đầu:

df%>% 
  group_by(Protocol)%>%
  summarize('mean' = mean(LH)) %>% knitr::kable()

Protocol	mean
Fol_long	1.366628
GnRH_ant	7.784449
Lut_short	2.503227
Mild	12.020800
PPOS	3.651961

exp(Intercept) tương ứng với trung bình của LH trong phác đồ Follicular long acting, đơn vị mUI/mL: \(exp(0.31235) = 1.366633\), ta thấy giá trị này tương đồng với kết quả mean LH ở nhóm Fol_Long

Giá trị trung bình LH nhóm GnRH_ant được tính bằng : exp(Intercept + ProtocolGnRH_ant): \(exp(0.31235 + 1.73978) = 7.784464\) …

Cứ thế, ta có thể ước lượng được trung bình của log(LH) hoặc LH ở cả 5 phân nhóm.

8 Bước 4: Suy diễn thống kê

Dựa vào 2 mô hình tuyến tính tổng quát này, ta có thể thực hiện một số suy diễn thống kê như sau:

Đầu tiên, khi so sánh mô hình chỉ chứa intercept và mô hình có thêm biến Protocol bằng một kiểm định Likelihood ratio, ta có thể biết liệu phân nhóm phác đồ thực sự có hiệu ứng ý nghĩa đối với LH hay không ? (câu hỏi này tương tự với kiểm định F theo Fisher trong phân tích phương sai cổ điển)

LR.test(m0,m1)

##  Likelihood Ratio Test for nested GAMLSS models. 
##  (No check whether the models are nested is performed). 
##  
##        Null model: deviance= 5005.425 with  2 deg. of freedom 
##  Altenative model: deviance= 4275.111 with  10 deg. of freedom 
##  
##  LRT = 730.3143 with 8 deg. of freedom and p-value= 0

Kết quả cho ra một giá trị p cực kì thấp, như vậy ta có thể khẳng định rằng phân nhóm phác đồ thực sự có tác động làm thay đổi LH.

Nói cách khác, có sự khác biệt ý nghĩa về giá trị LH giữa các nhóm phác đồ khác nhau. Tuy nhiên ta chưa thể biết rõ sự khác biệt này là ở nhóm nào so với nhóm nào.

Nội dung hệ số hồi quy trong mô hình cho phép ta kiểm định về sự khác biệt của LH trung bình ở phác đồ PPOS, GnRH, Luteal short acting và Mild (Letrozole) so với nhóm tham chiếu là Follicular long acting.

Dựa vào kết quả thô của mô hình, có thể thấy tất cả hệ số hồi quy đều >0 và có ý nghĩa thống kê, cho thấy giá trị trung bình LH ở 4 phác đồ còn lại đều cao hơn so với phác đồ tham chiếu (Longacting).

Tuy nhiên, ta sẽ không suy diễn trực tiếp dựa vào hệ số hồi quy, mà tiến hành một phân tích hậu kiểm với nội dung so sánh bắt cặp tuần tự giữa 5 loại phác đồ. Trong khi làm thống kê mô tả, ta có cảm nhận về khuynh hướng tăng dần của giá trị LH trung bình theo thứ tự: Fol_long < Short < PPOS < GnRH_ant < Mild

Ở đây ta sẽ kiểm tra lại kết quả này.

Để tiến hành phân tích hậu kiểm này, ta dùng hàm avg_comparisons của thư viện marginalseffect với tùy chỉnh variables = list(Protocol = ‘pairwise’), sau đó hiệu chỉnh giá trị p bằng hàm p.adjust trong R, với phương pháp hiệu chỉnh Bonferroni.

Kết quả của hàm này là một bảng như sau:

pw_comp = avg_comparisons(m1,
                          what = 'mu',
                variables = list(Protocol = 'pairwise'))

pw_comp$adj_p = p.adjust(pw_comp$p.value, method = 'bonferroni')

pw_comp %>% dplyr::select(-c(1,2,6,7)) %>%knitr::kable()

contrast	estimate	std.error	conf.low	conf.high	adj_p
GnRH_ant - Fol_long	6.417821	0.3715848	5.6895277	7.146114	0.0000000
Lut_short - Fol_long	1.136599	0.1227487	0.8960158	1.377182	0.0000000
Lut_short - GnRH_ant	-5.281222	0.3751863	-6.0165734	-4.545870	0.0000000
Mild - Fol_long	10.654172	0.9328177	8.8258825	12.482461	0.0000000
Mild - GnRH_ant	4.236351	0.9979211	2.2804617	6.192240	0.0002184
Mild - Lut_short	9.517573	0.9342582	7.6864603	11.348685	0.0000000
PPOS - Fol_long	2.285332	0.2290351	1.8364319	2.734233	0.0000000
PPOS - GnRH_ant	-4.132488	0.4220836	-4.9597569	-3.305220	0.0000000
PPOS - Lut_short	1.148734	0.2348331	0.6884692	1.608998	0.0000100
PPOS - Mild	-8.368839	0.9540588	-10.2387601	-6.498918	0.0000000

Cột contrast trình bày các cặp so sánh mà ta cần kiểm nghiệm, thí dụ GnRH_ant - Fol_long;

Cột estimate trình bày giá trị khác biệt trung bình về LH giữa 2 nhóm phác đồ, thí dụ GnRHant cao hơn Long acting 6.4178 mUI/mL.

2 cột conf.low conf.high trình bày ngưỡng dưới và trên khoảng tin cậy 95% của khác biệt, KTC95% có ý nghĩa thống kê khi không chứa giá trị 0.

Cột adj_p trình bày giá trị p sau hiệu chỉnh Bonferroni,

Kết quả cho thấy toàn bộ 10 giả thuyết đều có ý nghĩa thống kê. Một cách tổng quát, giá trị LH trung bình tăng dần theo thứ tự 5 phác đồ:

Fol_long < Short < PPOS < GnRH_ant < Mild

Như vậy phác đồ Mild stimulation với Letrozole có LH trung bình ngày trigger cao nhất, phác đồ Follicular long acting có LH trung bình thấp nhất, phác đồ PPOS ở vị trí trung gian, với LH cao hơn so với 2 phác đồ Long và Short.

Cuối cùng, kết quả của mô hình hồi quy có thể trình bày trực quan như sau:

m1%>%augment()%>%mutate(UL= .fitted +.se.fit*1.96,
                     LL=.fitted-.se.fit*1.96)%>%
  group_by(Protocol)%>%
  summarize_at('.fitted', median) -> m1_sum

augment(m1)%>%mutate(UL= .fitted +.se.fit*1.96,
                     LL=.fitted-.se.fit*1.96)%>%
  ggplot()+
  geom_jitter(aes(y = Protocol, 
                  x = LH,
                  col = Protocol),
              scale = 0.8,
              alpha = 0.3,
              width = 0.01,
              show.legend = F)+
  geom_density_ridges(aes(y = Protocol, 
                          x = LH,
                          fill = Protocol),
                      scale = 0.8,
                      alpha = 0.3,
                      show.legend = F)+
  geom_errorbar(aes(y=Protocol, 
                    xmin=exp(LL), 
                    xmax=exp(UL)),
                width=0.1,
                size=1) + 
  geom_path(data = m1_sum,
            aes(y=Protocol, 
                x=exp(.fitted),
                group = 1))+
  geom_point(aes(y=Protocol, 
                 x=exp(.fitted)), 
             size=3)+
  scale_x_continuous(breaks = seq(0,40,by = 5))+
  labs(x = 'LH level on trigger day (mIU/mL)', y = 'Protocols')+
  coord_flip()+
  theme_bw()  

9 Kết luận

Ta thấy mô hình GLM với phân phối Gamma là một giải pháp phổ quát và hiệu quả cho hầu hết những tình huống nghiên cứu y học lâm sàng. Mô hình GLM Gamma có thể thay thế hoàn toàn phân tích phương sai (ANOVA) cổ điển và cả mô hình hồi quy tuyến tính thông thường.