Mô phỏng phân phối Gamma từ y văn

1 Đặt vấn đề

Trong y văn, đặc tính phân phối của các đại lượng được nghiên cứu thường chỉ được mô tả một cách giản lược dưới hình thức trung bình (độ lệch chuẩn) hoặc trung vị (khoảng tứ phân vị). Những thông tin này chỉ cung cấp khái niệm về khuynh hướng trung tâm và mưc độ phân tán của dữ liệu, nhưng không đủ để tái hiện lại hình ảnh chính xác của phân phối thực tế.

Trên thực tế, các đại lượng y sinh học thường không khớp với phân phối chuẩn, nhưng có những đặc điểm như thang đo chỉ bao gồm giá trị >0, lệch phải, và đuôi dài. Ngoài ta, cơ chế tạo nên giá trị của các thông số này thường tương đồng với tiến trình tích lũy theo quy luật phân phối Gamma; do đó quy luật phân phối xác suất Gamma có thể phù hợp hơn để ước lượng chúng.

Trong một số hoàn cảnh, nhà khoa học có nhu cầu mô phỏng lại chính xác dữ liệu của một đại lượng dựa theo thông tin mô tả từ y văn, với tiêu chí càng gần với thực tế càng tốt. Việc mô phỏng này có rất nhiều ứng dụng thực tiễn, bao gồm ước tính cỡ mẫu, tối ưu hóa hiệu suất của một suy luận thống kê, thiết kế analysis plan, phân tích độ nhạy, phân tích tổng hợp, ước tính chi phí/hiệu quả của can thiệp điều trị, tối ưu hóa chính sách y tế, giảng dạy và học tập v.v.

Tuy nhiên, mô phỏng chính xác phân phối Gamma khi chỉ có trong tay thông tin về trung vị và khoảng tứ phân vị hay trung bình và độ lệch chuẩn là một bài toán không hề đơn giản, vì phân phối Gamma là một phân phối xác suất phức tạp được xác định bởi hệ tham số hoàn toàn khác so với phân phối chuẩn, và các hàm PDF, CDF và quantile của nó không có dạng đóng.

Bài thực hành sau sẽ cung cấp 2 giải pháp hiệu quả cho bài toán này.

2 Thí dụ minh họa

library(tidyverse)
library(fitdistrplus)
library(nleqslv)

Trong một nghiên cứu, người ta định lượng hormone LH (luteinizing hormone) vào ngày trigger trong một chu kì kích thích buồng trứng. Phân phối của dữ liệu thực nghiệm có hình ảnh như sau:

# Load data
df <- readxl::read_xlsx("Data_Qianwen_Xi.xlsx")%>%
  filter(Protocol == "PPOS")%>%
  dplyr::select(LH_TD)%>%
  na.omit()

# KDE plot
df %>% ggplot()+
  geom_density(aes(x = LH_TD), 
               fill ="red",
               alpha = 0.6)+
  theme_bw()

Có nhiều cơ sở cho thấy quy luật phân phối Gamma có thể phù hợp với dữ liệu này, bao gồm:

Lượng hormone là kết quả của quá trình tăng trưởng của các tế bào chế tiết, và được tích lũy theo thời gian tương ứng với tiến trình tăng trưởng theo quy luật hàm mũ;

Đại lượng sinh lý này không thể âm và có khuynh hướng lệch phải chứ không đối xứng;

Giả định này có thể được kiểm nghiệm một cách trực quan bằng cách khớp một phân phối Gamma với dữ liệu thực nghiệm bằng package fitdistrplus. Kết quả cho thấy quả thực phân phối Gamma khớp tốt với dữ liệu thực tế, với các tham số ước lượng là \(k = 3.010749\), \(\theta = 1.623632\) và \(\beta = 1/\theta = 0.61590\).

# Fit gamma
fit_gamma <- fitdist(df$LH_TD, "gamma")

# value of params
fit_gamma$estimate

##    shape     rate 
## 3.010749 1.623632

# Comparative plots with legend

df%>%ggplot(aes(x=LH_TD))+
  geom_density(aes(y = after_stat(density)), 
               fill = "gold",
               color = NA,
               linetype = 1,
               alpha = 0.6) +
  stat_function(fun = dgamma, 
                args = list(shape = fit_gamma$estimate[1], 
                            rate = fit_gamma$estimate[2]), 
                linetype = 2,
                color = "black") +
  theme_bw()

Tuy nhiên, trong bài báo, người ta chỉ báo cáo mean(sd) hoặc median(IQR) của giá trị LH như sau:

df %>% 
  summarise(median = median(LH_TD),
            iqr = IQR(LH_TD),
            mean = mean(LH_TD),
            sd = sd(LH_TD))%>%
  knitr::kable()

median	iqr	mean	sd
1.72	1.46	1.854264	1.042271

3 Phân tích bài toán

Phân phối Gamma là một phân phối xác suất liên tục, thường được tham số hóa bởi hai tham số: \(k\) (shape, hình dạng của phân phối, ảnh hưởng đến độ lệch và độ nhọn), \(\theta\) (scale) và \(\beta = 1/\theta\) (rate) (ảnh hưởng đến độ phân tán của phân phối), với hàm mật độ xác suất:

\[f(x; k, \theta) = \frac{1}{\Gamma(k) \theta^k} x^{k-1} e^{-x/\theta}, \quad x > 0, , k > 0, , \theta > 0\] Trong đó \(\Gamma(k)\) là hàm Gamma, một phiên bản tổng quát của giai thừa cho các số thực.

Mục tiêu của bài toán hiện thời là xác định các tham số \(k\) và \(\theta\) (hoặc \(\beta\)), sao cho phân phối Gamma tạo ra có trung vị và khoảng tứ phân vị, hoặc trung bình và độ lệch chuẩn khớp với các giá trị được báo cáo từ y văn.

4 Trường hợp biết trung bình và độ lệch chuẩn

Trong trường hợp biết trung bình và độ lệch chuẩn, giải pháp tương đối đơn giản vì các tham số của phân phối Gamma có thể được ước lượng dựa vào công thức tính giá trị kỳ vọng và phương sai của phân phối Gamma:

\[E(X) = k\theta, \quad Var(X) = k\theta^2\] Từ đó:

\[k = \frac{\mu^2}{\sigma^2}, \quad \theta = \frac{\sigma^2}{\mu}, \quad \beta = 1/\theta\] Ta viết hai function để mô phỏng dữ liệu từ trung bình và độ lệch chuẩn

param_from_mean_sd <- function(mu, sigma) {
  
  shape <- (mu^2)/(sigma^2)
  rate <- mu/sigma^2
  scale <- 1/rate
  
  return(list(shape = shape, 
              rate = rate, 
              scale = scale))
}

# Simulate data from mean and sd
sim_gamma_from_mean_sd <- function(mu, sigma, n) {
  
  params <- param_from_mean_sd(mu, sigma)
  
  sim_data <- rgamma(n, 
                     shape = params[['shape']], 
                     rate = params[['rate']])
  
  return(sim_data)
}

Theo kết quả mô tả ở trên, trung bình và độ lệch chuẩn của giá trị LH là 1.854264 và 1.042271, ta có thể mô phỏng dữ liệu từ thông tin này.

Kết quả cho thấy dữ liệu mô phỏng (màu xanh) khá gần với dữ liệu thực tế (màu đỏ)

set.seed(123)

sim_x <- sim_gamma_from_mean_sd(1.854264, 1.042271,1000)

df %>% ggplot()+
  geom_density(aes(x = LH_TD, y = after_stat(density)), 
               fill = "red",
               color = NA,
               linetype = 1,
               alpha = 0.3) +
   geom_density(data = tibble(sim_x),
               aes(x = sim_x), 
               fill ="blue",
               color = NA,
               alpha = 0.3)+
  theme_bw()

5 Trường hợp biết trung vị và khoảng tứ phân vị

Trường hợp này thực ra là khó hơn nhiều, vì 2 thông số này không cho phép ước lượng trực tiếp các tham số của phân phối Gamma.

5.1 Giải pháp 1: khớp phân vị

Bài toán này có thể giải bằng phương pháp “khớp phân vị” (quantiles matching). Đầu tiên, ta tính trung vị, và IQR (bách phân vị thứ 25 và 75) từ một giá trị \(k\) và \(\theta\) bất kì, cách làm là sử dụng nghịch đảo của hàm CDF, hay hàm quantile của phân phối Gamma.

Như vậy, ta tạo ra một hệ phương trình phi tuyến 2 ẩn:

\[\begin{cases} qgamma(q=0.5, shape = k, scale = \theta) - median_{k,\theta} = 0 \\ qgamma(q=0.75, shape = k, scale = \theta) - qgamma(q=0.25, shape = k, scale = \theta) - IQR_{k,\theta} = 0\end{cases}\]

Vì qgamma là một hàm phi tuyến phức tạp (không có nghiệm dạng đóng), ta phải sử dụng thuật toán số học để giải hệ phương trình \(F(k,θ)=0\)

Package nleqslv trong R cung cấp phương pháp giải hệ phương trình phi tuyến dựa trên các thuật toán như Newton-Raphson hoặc Broyden. Các bước cơ bản của thuật toán:

• Khởi tạo giá trị ban đầu: Chọn một điểm khởi đầu \(x_0 = [k_0,\theta_0]\) gần với nghiệm thực tế để đảm bảo hội tụ.

• Sử dụng phương pháp lặp để cập nhật giá trị \(x_i\) từ \(x_{i-1}\), với công thức: \(x_i = x_{i-1} - J^{-1}F(x_{i-1})\), trong đó \(J\) là ma trận Jacobian của hàm \(F\).

• Lặp lại quá trình trên cho đến khi hội tụ, tức là \(F(x_i)\) đủ nhỏ.

Sau khi giải hệ phương trình, ta thu được giá trị \(k\) và \(\theta\) tối ưu, từ đó ta có thể tạo ra một mẫu dữ liệu tuân theo phân phối Gamma với các tham số này.

Ta viết một function để triển khai phương pháp này:

quantile_matching_nleqslv <- function(m, iqr) {
  
  sys_fun <- function(x) {
    k <- x[1]
    theta <- x[2]
    
    eq1 <- qgamma(0.5, shape = k, scale = theta) - m
    
    eq2 <- qgamma(0.75, shape = k,scale = theta) - qgamma(0.25, shape = k, scale = theta) - iqr
    
    return(c(eq1, eq2))
  }
  
  init_k <- 1
  init_theta <- m / qgamma(0.5, shape = init_k, scale = 1)
  
  x_init <- c(init_k, init_theta)
  
  sol <- nleqslv(x_init, sys_fun)
  if (sol$termcd != 1) warning("nleqslv did not converge.")
  
  return(list(shape = sol$x[1], 
              rate = 1/sol$x[2], 
              scale = sol$x[2]))
}

Áp dụng function này cho trung vị và IQR của dữ liệu LH, kết quả là ta cũng thu được kết quả khá gần với thực tế:

set.seed(123)

params <- quantile_matching_nleqslv(1.72, 1.46)

sim_x <- rgamma(1000, 
                shape = params[['shape']], 
                rate = params[['rate']])

df %>% ggplot()+
  geom_density(aes(x = LH_TD, y = after_stat(density)), 
               fill = "red",
               color = NA,
               linetype = 1,
               alpha = 0.3) +
    geom_density(data = tibble(sim_x),
               aes(x = sim_x), 
               fill ="blue",
               color = NA,
               alpha = 0.3)+
  theme_bw()

5.2 Giải pháp khác: khớp tỷ số

Một phương pháp khác giản lược hơn để giải bài toán này là “khớp tỷ số” (ratio matching).

Phương pháp này tận dụng đặc tính của phân phối Gamma: tỷ số giữa median và IQR phụ thuộc duy nhất vào tham số \(k\) (shape) khi \(\theta\) được chuẩn hóa (scale = 1). Khi so sánh tỷ số này với tỷ số được cung cấp từ y văn, ta tạo ra một phương trình phi tuyến với một ẩn số là \(k\):

\[f(k) = \frac{median_k}{IQR_k} - \frac{median_{true}}{IQR_{true}}\] Như vậy ta đã giản lược bài toán từ hệ 2 phương trình xuống còn 1 phương trình và một ẩn, tuy nhiên phương trình này cũng không có công thức nghiệm đóng, nó chỉ có thể được giải bằng phương pháp số học.

Ví dụ, ta có thể dùng thuật toán uniroot trong R với cơ chế như sau:

• Cung cấp một khoảng ước lượng chứa nghiệm \([l,u]\), với \(l\) và \(u\) là hai giá trị của k, và ta biết rằng nghiệm k nằm trong khoảng này.

• Chia đôi tuần tự khoảng này thành nhiều khoảng nhỏ hơn, rồi tính giá trị của hàm \(f(m)\) tại giữa mỗi khoảng, với \(m = (l+u)/2\). Nếu \(f(m)\) có dấu trái với \(f(l)\), ta chuyển \(u\) thành \(m\), ngược lại, ta chuyển \(l\) thành \(m\).

• Lặp lại quá trình trên cho đến khi khoảng \([l,u]\) đủ nhỏ, và ta coi \(m\) là nghiệm của phương trình.

Sau khi giải được \(k\), ta có thể dùng \(k\) và tỷ số \(\frac{median_k}{IQR_k}\) để giải hệ phương trình ban đầu, từ đó tìm ra \(\theta\) hoặc \(\beta\).

Ta viết R function để triển khai giải pháp này. Kết quả của 3 tham số hoàn toàn giống với cách làm thứ nhất, và dữ liệu mô phỏng ra cũng gần với thực tế:

# Quantile matching using uniroot

quantile_match_uniroot <- function(m,iqr){
  
  f <- function(k){
    
    med_k <- qgamma(0.5, shape = k, scale = 1)
    iqr_k <- qgamma(0.75, shape = k, scale = 1) - qgamma(0.25, shape = k, scale = 1)
    
    ratio_k <- med_k/iqr_k
    target_ratio <- m/iqr
    
    return(ratio_k - target_ratio)
  }
  
  alpha_initial <- (1.35 * m / iqr)^2
  lower_k <- max(0.1, alpha_initial - 5)
  upper_k <- alpha_initial + 5
  
  k_sol <- uniroot(f, 
                   lower = lower_k, 
                   upper = upper_k)$root
  
  theta_sol <- m / qgamma(0.5, shape = k_sol, scale = 1)
  
  return(list(shape = k_sol, 
              rate = 1/theta_sol, 
              scale = theta_sol))
}

solved_d <- quantile_match_uniroot(1.72, 1.46)

# Simulate data with the solved k and theta

sim_x <- rgamma(1000, 
                   shape = solved_d[['shape']], 
                   rate = solved_d[['rate']])

# compared 2 KDE

df %>% ggplot()+
  geom_density(aes(x = LH_TD, y = after_stat(density)), 
               fill = "red",
               color = NA,
               linetype = 1,
               alpha = 0.3) +
    geom_density(data = tibble(sim_x),
               aes(x = sim_x), 
               fill ="blue",
               color = NA,
               alpha = 0.3)+
  theme_bw()

6 Kết luận

Qua bài thực hành này, chúng ta đã có thể mô phỏng phân phối Gamma từ thông tin mô tả trong y văn, thông qua 2 phương pháp: khớp trung bình và độ lệch chuẩn, và khớp trung vị và khoảng tứ phân vị. Cả hai phương pháp đều cho kết quả khá gần với dữ liệu thực tế, và có thể áp dụng trong nhiều tình huống thực tiễn khác nhau.