Ước tính cỡ mẫu bằng mô phỏng
BS. Lê Ngọc Khả Nhi
24 Tháng 3 năm 2025
1 Giới thiệu
Ước tính cỡ mẫu là một bước quan trọng trong thiết kế nghiên cứu, giúp xác định số lượng cá thể cần thiết để đạt được độ chính xác tối ưu cho một kết quả suy diễn thống kê. Trong nghiên cứu lâm sàng, mục tiêu ước tính cỡ mẫu thường nhắm đến một giả thuyết về ưu thế hoặc tính tương đương về hiệu quả của một can thiệp điều trị so với một can thiệp khác.
Trong một số hoàn cảnh, bài toán thống kê có thể phức tạp và không có sẵn công thức tính toán trực tiếp phù hợp. Lúc này, mô phỏng là một phương pháp phổ quát để ước lượng cỡ mẫu. Bài thực hành này sẽ hướng dẫn các bạn sử dụng mô phỏng để ước tính cỡ mẫu cho một bài toán kiểm định không thua kém (non-inferiority testing) dựa trên phân tích hồi quy với phân phối gamma.
Trong thí dụ minh họa này, một bác sĩ thuộc chuyên khoa hiếm muộn dự định tiến hành một nghiên cứu lâm sàng để chứng minh giả thuyết rằng phác đồ Progestin-primed ovarian stimulation (PPOS) cũng hiệu quả không thua kém so với phác đồ GnRH antagonist về khả năng kiểm soát hormone LH trong chu kỳ kích thích buồng trứng ở phụ nữ vô sinh. Theo dự kiến, thí nghiệm sẽ được tiến hành trên 2 nhóm: đối chứng (phác đồ GnRh_ant) và can thiệp (phác đồ PPOS). Kết cục chính của thí nghiệm là nồng độ LH ngày trigger (LH_TD, mIU/ml) ở mỗi nhóm. Giả thuyết cần chứng minh là giá trị LH_TD trung bình của nhóm can thiệp không cao hơn nhóm đối chứng quá một biên \(\delta\) = 1.5 mIU/ml. Bác sĩ này dự định dùng phân tích hồi quy với phân phối Gamma để ước lượng hiệu ứng điều trị và chứng minh giả thuyết của mình. Vấn đề đặt ra là cần ước lượng cỡ mẫu cần thiết để đạt được mức ý nghĩa 0.025 và power 0.8 cho kiểm định không thua kém.
Bác sĩ này cũng có trong tay một số dữ kiện từ y văn về giá trị trung bình và độ lệch chuẩn của mức LH ngày trigger cho phác đồ quy ước và phác đồ PPOS. Một cách cụ thể:
Nhóm đối chứng: \(\mu_0\) = 0.6725, \(\sigma_0\) = 0.526
Nhóm can thiệp: \(\mu_1\) = 1.8542, \(\sigma_1\) = 1.04
2 Phân tích bài toán
Bản chất của bài toán thống kê hiện thời là so sánh một biến định lượng giữa 2 phân nhóm độc lập và kiểm định giả thuyết về sự khác biệt giữa 2 trung bình này. Lẽ ra ta có thể dễ dàng áp dụng một công thức quy ước. Tuy nhiên, có một số vấn đề của kế hoạch phân tích khiến nó trở nên phức tạp hơn:
Biến kết quả (nồng độ LH) chắc chắn không thể có phân phối chuẩn, nhưng cần được ước lượng bằng một phân phối phức tạp hơn là gamma;
Suy diễn thống kê cho hiệu ứng can thiệp được thực hiện thông qua mô hình hồi quy gamma với hàm liên kết logarit;
Giả thuyết cần kiểm định có tính “không thua kém” (non-inferiority), tức là chứng minh rằng hiệu ứng can thiệp không kém hơn nhóm đối chứng quá một biên \(\delta\);
Ta thấy rằng không có sẵn công thức nào cho bài toán ước tính cỡ mẫu này.
3 Giải pháp mô phỏng
Một giải pháp khả thi cho bài toán ước tính cỡ mẫu này, đó là sử dụng dữ liệu mô phỏng. Vì bản chất của việc ước tính cỡ mẫu là xác định số lượng cá thể cần thiết để đạt được mức statistical power mong muốn, ta hãy tưởng tượng nếu lặp lại rất nhiều lần phân tích thống kê trên các mẫu có kích thước tăng dần, và tính power cho mỗi lượt, kết quả sẽ tạo ra một đồ thị có hình ảnh đường cong, gọi là power curve, diễn tả quy luật biến thiên của power tùy theo kích thước mẫu. Ta có thể dùng nó làm thước đo để ước lượng cỡ mẫu tối thiểu để đạt mức power mong muốn, ví dụ 0.8.
Mô phỏng cho phép ta tạo ra dữ liệu cho mỗi lượt phân tích thống kê trong quy trình này.
library(tidyverse)
library(fitdistrplus)
library(progressr) # Progress bar
library(furrr) # Parallel processing with purrr
library(future.apply)
library(marginaleffects)
plan(multisession, workers = parallel::detectCores() - 1)4 Tình huống thứ nhất: Có thông tin của cả hai nhóm
4.1 Cách làm thứ nhất: Mô phỏng trực tiếp và lặp lại nhiều lần
Trong tình huống này, ta có thông tin về cả hai nhóm đối chứng và can thiệp, bao gồm giá trị trung bình và độ lệch chuẩn của mức LH ngày trigger. Ta sẽ sử dụng thông tin này để mô phỏng dữ liệu cho mỗi lượt và lặp lại nhiều lần để ước lượng power curve.
Quy trình thực hiện như sau:
Trước tiên, chúng ta cần các tham số từ tài liệu, bao gồm:
Trung bình và độ lệch chuẩn của Y cho nhóm đối chứng (X=0): \(\mu_0 = 0.6725\), \(\sigma_0 = 0.526\).
Trung bình và độ lệch chuẩn của Y cho nhóm điều trị (X=1): \(\mu_1 = 1.8542\), \(\sigma_1=1.04\).
Biên độ không thua kém \(\delta - 1.5\)
Để sử dụng phân phối gamma, chúng ta cần chuyển đổi \(\mu\) và \(\sigma\) thành tham số shape \(\alpha\) và scale \(\theta\) (xem bài trước dể biết cách làm), cho mỗi nhóm
Với mỗi kích thước mẫu n (ví dụ, từ 30 đến 500, tăng dần theo bước 5), chúng ta thực hiện 100 lần mô phỏng:
Tạo dữ liệu cho nhóm đối chứng: Y0 = rgamma(n/2, shape = \(\alpha_0\), scale = \(\theta_0\))
Tạo dữ liệu cho nhóm điều trị: Y1 = rgamma(n/2, shape = \(\alpha_1\), scale = \(\theta_1\))
Kết hợp dữ liệu thành một khung dữ liệu với biến Y và X (X=0 cho nhóm đối chứng, X=1 cho nhóm điều trị), tỷ lệ 1:1.
Với mỗi lần mô phỏng, ta khớp một mô hình GLM với family=gamma và hàm liên kết logarit, với công thức: \(log(E[Y]) \sim \beta_0 + \beta_1X\)
Mô hình này cho phép ước lượng hiệu ứng can thiệp trung bình dưới dạng:
\(E[Y|X=0] = exp(\beta_0)\) là giá trị LH trung bình của nhóm đối chứng.
\(E[Y|X=1] = exp(\beta_0 + \beta_1)\) là giá trị LH trung bình của nhóm can thiệp.
Hiệu ứng can thiệp là khác biệt giữa 2 trị số này:
\(diff = exp(\beta_0 + \beta_1) - exp(\beta_0) = exp(\beta_0)(exp(\beta_1) - 1)\)
Kiểm định không thua kém (non-inferiority testing) sẽ được thực hiện, nhằm chứng minh rằng hiệu quả phác đồ điều trị không kém hơn nhóm đối chứng, tức là giá trị Y của nhóm điều trị không cao hơn nhóm đối chứng quá biên \(\delta\)
Với giả định Y thấp hơn là tốt hơn,ta sẽ so sánh giới hạn trên của hiệu ứng can thiệp với biên độ \(\delta\). Nếu giới hạn trên nhỏ hơn hoặc bằng \(\delta\), ta kết luận rằng hiệu ứng can thiệp không thua kém so với nhóm đối chứng.
Giả thuyết H0: \(diff > \delta\); Giả thuyết H1: \(diff \leq \delta\)
Đây là kiểm định một phía (left-tailed) với mức ý nghĩa 0.025, nghĩa là xác suất sai lầm loại I là 2.5%.
Tuy nhiên, để có thể sử dụng biên độ \(\delta\) cho kiểm định không thua kém dựa trên đơn vị đo nguyên bản của LH (không dùng hàm log hay hàm exponential), ta cần tính toán giới hạn trên cho hiệu ứng can thiệp trực tiếp bằng phương pháp delta.
Đạo hàm theo \(\beta_0\): \(exp(\beta_0)(exp(\beta_1)-1)\)
Đạo hàm theo \(\beta_1\): \(exp(\beta_0)exp(\beta_1)\)
Sai số chuẩn (SE) của hiệu ứng được tính từ ma trận hiệp phương sai (vcov) của các hệ số trong mô hình, với công thức:
\(SE_d^2 = (\frac {\partial diff}{\partial \beta_0})^2.Var(\beta_0) + (\frac {\partial diff}{\partial \beta_1})^2.Var(\beta_1) + 2.\frac {\partial diff}{\partial \beta_0}.\frac {\partial diff}{\partial \beta_1}.Cov(\beta_0,\beta_1)\)
Giá trị SE và \(diff\) cho phép ước tính được giới hạn trên (UL) của khoảng tin cậy 95% cho hiệu ứng can thiệp (hay phân vị thứ 97.5 của phân phối hiệu ứng):
\(UL = diff + qnorm(0.975).SE_d\), với \(qnorm(0.975) \approx 1.96\)
Kiểm định không thua kém sẽ được thực hiện bằng cách so sánh giới hạn trên này với biên độ \(\delta\). Nếu \(UL \leq \delta\), ta kết luận rằng hiệu ứng can thiệp không thua kém so với nhóm đối chứng.
Power của kiểm định sẽ được ước lượng bằng tỷ lệ số lần giả thuyết H1 được chấp nhận trong 100 lần mô phỏng.
Ta áp dụng quy trình này cho các giá trị kích thước mẫu \(n\) tăng dần, ví dụ từ 30 đến 500, với bước nhảy là 5. Kết quả thu được là một dataframe ghi lại power cho mỗi bậc giá trị n.
Đồ thị mô tả sự biến thiên của power theo kích thước mẫu (power curve) cho phép ta ước lượng cỡ mẫu cần thiết để đạt được mức power kỳ vọng, ví dụ 0.8
Sau đây là R script để triển khai quy trình này (lưu ý: code được viết theo lập trình hàm thay vì vòng lặp cơ bản, và cho phép tính toán song song nhiều CPU, sử dụng packages furrr và future.apply):
Bước 1:Định nghĩa tham số:
mu0, sigma0: Trung bình và độ lệch chuẩn của nhóm đối chứng.
mu1, sigma1: Trung bình và độ lệch chuẩn của nhóm điều trị.
delta: Biên không thua kém
n_sim: Số lần mô phỏng mỗi bậc n (100)
n_seq: Dãy kích thước mẫu
mu0 <- 0.6725 # Mean of control group (X=0)
sigma0 <- 0.526 # SD of control group (X=0)
mu1 <- 1.854264 # Mean of treatment group (X=1)
sigma1 <- 1.04 # SD of treatment group (X=1)
# alpha <- 0.025 # Significance level for one-sided test
delta <- 1.5 # Non-inferiority margin
n_sim <- 100 # Number of simulations per sample size
n_seq <- seq(30, 500, by = 5) # Sequence of sample sizes to testViết 1 hàm để chuyển đổi trung bình và độ lệch chuẩn sang tham số shape và scale của phân phối gamma:
gamma_params <- function(mu, sigma) {
shape <- mu^2 / sigma^2
scale <- sigma^2 / mu
return(list(shape = shape, scale = scale))
}
params0 <- gamma_params(mu0, sigma0)
params1 <- gamma_params(mu1, sigma1)Bước 2: Hàm mô phỏng và kiểm định cho mỗi lượt
Tạo n quan sát từ phân phối gamma cho nhóm đối chứng (Y0) và điều trị (Y1).
Tạo dataframe với Y (kết quả) và X (0 cho đối chứng, 1 cho điều trị), tỷ lệ 1:1.
Khớp mô hình GLM với family=gamma và hàm liên kết logarit
Tính hiệu ứng can thiệp trung bình và sai số chuẩn
Tính giới hạn trên của hiệu ứng can thiệp và kiểm định không thua kém
Trả về TRUE nếu giới hạn trên nhỏ hơn hoặc bằng delta, FALSE nếu ngược lại.
simulate_and_test_simple <- function(n) {
data <- tibble(
X = rep(c(0, 1), each = n / 2),
Y = c(rgamma(n / 2, shape = params0$shape, scale = params0$scale),
rgamma(n / 2, shape = params1$shape, scale = params1$scale))
)
fit <- tryCatch(
glm(Y ~ X, family = Gamma(link = "log"),
data = data),
error = function(e) NULL
)
if (is.null(fit)) return(NA)
#coefs <- coef(fit)
#vcov_mat <- vcov(fit)
#beta0 <- coefs[1]
#beta1 <- coefs[2]
#diff_est <- exp(beta0 + beta1) - exp(beta0)
#grad <- c(exp(beta0 + beta1) - exp(beta0), exp(beta0 + beta1))
#se_diff <- sqrt(t(grad) %*% vcov_mat %*% grad)
#ul <- diff_est + qnorm(1 - alpha) * se_diff
out <- avg_comparisons(fit)
ul <- out$conf.high
return(ul <= delta)
}Bước 3: Hàm chạy mô phỏng cho mỗi kích thước mẫu
run_simulations_simple <- function(n, n_sim) {
set.seed(123)
results <- future_replicate(n_sim, simulate_and_test_simple(n), future.seed = TRUE)
power <- mean(results, na.rm = TRUE)
return(tibble(n = n, power = power))
}Bước 4: Mô phỏng và tính toán power cho mỗi bậc n
with_progress({
p <- progressor(along = n_seq)
power_curve <- map_dfr(n_seq, ~{
p()
run_simulations_simple(.x, n_sim)
})
})Bước 5: Vẽ power curve
ggplot(power_curve, aes(x = n, y = power)) +
geom_smooth(fill = "gold",
linewidth = 0.8,
alpha = 0.5,
color = "black") +
geom_point(size = 1, alpha = 0.3) +
geom_hline(yintercept = 0.8,
linetype = "dashed",
color = "red3") +
geom_vline(xintercept = 295,
linetype = "dashed",
color = "black") +
scale_x_continuous(breaks = seq(0, 500, by = 50)) +
scale_y_continuous(breaks = seq(0, 1, by = 0.1)) +
labs(title = "Power Curve for Non-Inferiority Test",
x = "Sample Size (n)",
y = "Power") +
theme_minimal()Kết quả cho thấy power curve cho thấy mức power tăng dần theo kích thước mẫu, và đạt mức 0.8 ở khoảng 296 (148*2) cá thể. Điều này có nghĩa rằng cỡ mẫu tối thiểu cần thiết để đạt được mức ý nghĩa 0.025 và power 0.8 cho kiểm định không thua kém là 148 bệnh nhân cho mỗi nhóm.
5 Cách làm khác: Kết hợp mô phỏng và tái chọn mẫu
Một cách làm khác, hiệu quả hơn về mặt tính toán, là tạo ra một tập dữ liệu rất lớn (1 triệu cá thể) cố định từ phân phối gamma cho cả hai nhóm, đại diện cho quần thể, sau đó tái chọn mẫu cho mỗi lượt mô phỏng. Cách làm này giảm thiểu sự biến động ngẫu nhiên giữa các lượt mô phỏng và tăng hiệu quả tính toán.
mu0 <- 0.6725 # Mean of control group (X=0)
sigma0 <- 0.526 # SD of control group (X=0)
mu1 <- 1.854264 # Mean of treatment group (X=1)
sigma1 <- 1.04 # SD of treatment group (X=1)
delta <- 1.5 # Non-inferiority margin
# alpha <- 0.025 # Significance level for one-sided test
n_sim <- 100 # Number of simulations per sample size
n_values <- seq(30, 500, by = 5) # Sequence of sample sizes to test
theta0 <- sigma0^2 / mu0 # Scale parameter for control
k0 <- mu0^2 / sigma0^2 # Shape parameter for control
theta1 <- sigma1^2 / mu1 # Scale parameter for treatment
k1 <- mu1^2 / sigma1^2 # Shape parameter for treatment
set.seed(123) # For reproducibility
pop_control <- rgamma(1e6, shape = k0, scale = theta0) # Control population
pop_treatment <- rgamma(1e6, shape = k1, scale = theta1) # Treatment populationsimulate_power <- function(n, p) {
# Inner progress bar for simulations within each n
p_inner <- progressor(along = 1:n_sim)
# Simulate n_sim times
results <- future_map(1:n_sim, ~{
# Update inner progress bar
p_inner(sprintf("n = %d", n))
# Sample from populations
sample_control <- sample(pop_control, size = n/2, replace = T)
sample_treatment <- sample(pop_treatment, size = n/2, replace = T)
# Combine into data frame
df <- tibble(
y = c(sample_control, sample_treatment),
x = rep(c(0, 1), each = n/2)
)
fit <- tryCatch(
glm(y ~ x, family = Gamma(link = "log"),
data = df),
error = function(e) NULL
)
if (is.null(fit)) return(NA)
#coefs <- coef(fit)
#vcov_mat <- vcov(fit)
#beta0 <- coefs[1]
#beta1 <- coefs[2]
#diff_est <- exp(beta0 + beta1) - exp(beta0)
#grad <- c(exp(beta0 + beta1) - exp(beta0), exp(beta0 + beta1))
#se_diff <- sqrt(t(grad) %*% vcov_mat %*% grad)
#upper_limit <- diff_est + qnorm(1 - alpha) * se_diff
out <- avg_comparisons(fit)
ul <- out$conf.high
return(ul <= delta)
}, .options = furrr_options(seed = TRUE))
# Calculate power as proportion of successes
mean(unlist(results))
}plan(multisession, workers = parallel::detectCores() - 1) # Use all but one core
with_progress({
# Outer progress bar for sample sizes
p <- progressor(along = n_values)
power_comb <- future_map_dbl(n_values, ~{
# Update outer progress bar
p(sprintf("Sample size = %d", .x))
simulate_power(.x, p)
}, .options = furrr_options(seed = TRUE))
})
power_df <- tibble(
n = n_values,
power = power_comb
)power_df%>%ggplot(aes(x = n, y = power)) +
geom_smooth(fill = "gold",
linewidth = 0.8,
alpha = 0.5,
color = "black") +
geom_point(size = 1, alpha = 0.3) +
geom_hline(yintercept = 0.8,
linetype = "dashed",
color = "red3") +
geom_vline(xintercept = 276,
linetype = "dashed",
color = "black") +
scale_x_continuous(breaks = seq(0, 500, by = 50)) +
scale_y_continuous(breaks = seq(0, 1, by = 0.1)) +
labs(title = "Power Curve for Non-Inferiority Test",
x = "Sample Size (n per group)",
y = "Power") +
theme_minimal()Cách làm này cho ra kết quả tương tự như cách làm trước, với cỡ mẫu tối thiểu cần thiết để đạt được mức power 0.8 là 276 cá thể (138 cho mỗi nhóm).
6 Phân tích ưu và nhược điểm của hai cách làm
Với cách làm thứ nhất, ta mô phỏng trực tiếp bằng hàm rgamma cho mỗi lượt và lặp lại nhiều lần. Cách này có một số ưu điểm:
Chính xác về mặt lý thuyết: Mỗi mẫu được sinh trực tiếp từ phân phối gamma lý thuyết, đảm bảo tính độc lập giữa các mẫu và phân phối chính xác.
Không phụ thuộc bộ nhớ: Không cần lưu trữ object dữ liệu kích thước lớn, chỉ sinh dữ liệu có tính tạm thời theo nhu cầu và được dọn dẹp sau khi xong quy trình, phù hợp với máy tính có RAM hạn chế.
Linh hoạt: Dễ dàng điều chỉnh tham số của phân phối gamma (shape, rate) trong quá trình mô phỏng nếu cần thử nghiệm nhiều kịch bản khác nhau.
Tuy nhiên, cách làm này cũng có một số nhược điểm:
Hiệu suất thấp: Gọi rgamma nhiều lần tốn nhiều thời gian tính toán hơn so với resampling từ một tập dữ liệu cố định, đặc biệt khi số lượt mô phỏng lớn.
Biến động ngẫu nhiên: Mỗi lần gọi rgamma đều tạo ra dữ liệu mới, dẫn đến sự khác biệt ngẫu nhiên giữa các lượt mô phỏng. Điều này có thể làm tăng độ biến thiên trong kết quả.
Ngược lại, cách làm thứ hai ta chỉ thực hiện mô phỏng một lần để tạo một tập dữ liệu cố định và rất lớn (ví dụ 1 triệu đơn vị quan sát) - đại diện cho dữ liệu quần thể, sau đó tái chọn mẫu cho mỗi lượt. Cách làm này có một số ưu điểm:
Tính nhất quán: Vì dữ liệu quần thể được tạo một lần duy nhất nên giảm thiểu sự biến động ngẫu nhiên giữa các lần tạo mẫu, giúp kết quả ổn định hơn. Hành động chọn mẫu con từ quần thể cũng phù hợp về mặt lý thuyết với quy trình chọn mẫu trong thực tế (mẫu nghiên cứu là phần được rút ngẫu nhiên từ quần thể xác định)
Rất hiệu quả về mặt thời gian tính toán khi n nhỏ so với kích thước quần thể, vì resampling nhanh hơn nhiều so với việc sinh dữ liệu mới từ rgamma (đặc biệt khi dùng replace = TRUE).
Tái lập dễ dàng: Với một tập dữ liệu lớn cố định, kết quả mô phỏng có thể được tái lập chính xác hơn khi kết hợp với set.seed(), vì không có thêm biến động từ việc sinh dữ liệu mới.
Tuy nhiên, cách làm này cũng có một số nhược điểm:
Dữ liệu quần thể cần đảm bảo đủ lớn để đại diện đầy đủ cho phân phối gamma. Nếu kích thước không đủ lớn thì các mẫu con có thể không phản ánh chính xác toàn bộ quần thể.
Tốn kém bộ nhớ: cần lưu trữ hai object dữ liệu lớn trong RAM (trường hợp này ta có 2 object với 1 triệu phần tử, tốn khoảng 16 Mb).
Như vậy:
Cách 1 chỉ có lợi thế nhỏ về độ chính xác lý thuyết, nhưng với N lớn, cách 2 gần tương đương và đủ chính xác cho hầu hết các ứng dụng thực tế.
Cách 2 vượt trội về hiệu quả tính toán, đặc biệt khi kết hợp với tính toán song song (furrr), miễn là máy tính có đủ RAM.
Như vậy cách 2 nên dùng khi chỉ cần mô phỏng cho một bộ tham số cố định, số bậc n lớn và cần tối ưu thời gian tính toán. Cách 1 chỉ nên dùng khi cần thử nghiệm nhiều kịch bản tham số khác nhau
7 Tình huống thứ hai: Chỉ có thông tin của nhóm đối chứng
Vì không có dữ liệu từ nhóm điều trị, chúng ta cần phải giả định nhiều kịch bản khác nhau cho giá trị tham số \(\mu_1\) và \(\sigma_1\), là trung bình LH ngày trigger ở nhóm can thiệp, và thực hiện mô phỏng để dựng power curve.
Với mỗi kịch bản giả định, ta tiến hành mô phỏng 100 lần phân phối gamma với kích thước mẫu \(n\) tăng dần, tỷ lệ 1:1 giữa hai nhóm, sử dụng tham số \(\mu_0\), \(\sigma_0\) (nhóm đối chứng) từ y văn và tham số \(\mu_1\), \(\sigma_1\) từ kịch bản giả định.
Như vậy quy trình chỉ phức tạp hơn một chút với vòng lặp bổ sung cho các kịch bản giả định, nhưng vẫn giữ nguyên cấu trúc chung của quy trình mô phỏng.
Do phải mô phỏng cho nhiều kịch bản, ta không dùng cách làm thứ hai với tập dữ liệu cố định mà sẽ sử dụng cách làm thứ nhất với mỗi lượt mô phỏng tương ứng với mỗi kịch bản.
Dưới đây là R script cho quy trình này:
Bước 1: Định nghĩa tham số và kịch bản giả định
mu0 <- 0.6725 # Mean of control group (X=0)
sigma0 <- 0.526 # SD of control group (X=0)
#mu1 <- 1.854264 # Mean of treatment group (X=1)
#sigma1 <- 1.04 # SD of treatment group (X=1)
delta <- 1.5 # Non-inferiority margin
#alpha <- 0.025 # Significance level for one-sided test
n_sim <- 100 # Number of simulations per sample size
sample_sizes <- seq(40, 500, by = 10) # Sequence of sample sizes to test
scenarios <- tibble(
mu1 = c(1.,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9), # Different means for treatment group
sigma1 = 1. # Assume same SD as control for simplicity
)Bước 2: Hàm mô phỏng và kiểm định cho mỗi lượt với kịch bản giả định
# Function to calculate gamma parameters with validation
calc_gamma_params <- function(mu, sigma) {
if (!is.numeric(mu) || !is.numeric(sigma) || mu <= 0 || sigma <= 0) {
stop("Mean and SD must be positive numbers")
}
variance <- sigma^2
shape <- mu^2 / variance
rate <- mu / variance
if (shape <= 0 || rate <= 0) stop("Invalid gamma parameters")
list(shape = shape, rate = rate)
}# Function to simulate data and perform non-inferiority test
simulate_and_test <- function(n, mu0, sigma0, mu1, sigma1, delta, alpha) {
# Calculate gamma parameters for both groups
params0 <- calc_gamma_params(mu0, sigma0)
params1 <- calc_gamma_params(mu1, sigma1)
# Simulate data with suppression of warnings
y0 <- suppressWarnings(rgamma(n / 2, shape = params0$shape, rate = params0$rate))
y1 <- suppressWarnings(rgamma(n / 2, shape = params1$shape, rate = params1$rate))
# Replace any NA or Inf with a small positive value
y0 <- ifelse(is.na(y0) | is.infinite(y0), 0.001, y0)
y1 <- ifelse(is.na(y1) | is.infinite(y1), 0.001, y1)
# Ensure all values are positive
y0 <- pmax(y0, 0.001) # Avoid zero or negative values
y1 <- pmax(y1, 0.001)
data <- tibble(
y = c(y0, y1),
x = c(rep(0, n / 2), rep(1, n / 2))
)
# Fit GLM with gamma distribution and log link
fit <- tryCatch(
glm(y ~ x, family = Gamma(link = "log"), data = data, control = list(maxit = 50)),
error = function(e) NULL
)
if (is.null(fit)) return(NA) # Return NA if model fails
out <- avg_comparisons(fit)
ul <- out$conf.high
# Non-inferiority test: TRUE if upper limit <= delta
return(ul <= delta)
}Bước 3: Hàm chạy mô phỏng cho mỗi kích thước mẫu và kịch bản giả định
# Main simulation function with progress bar
run_simulation <- function(n, mu0, sigma0,
mu1, sigma1,
delta,
alpha, n_sim) {
results <- future_replicate(n_sim,
simulate_and_test(n, mu0, sigma0, mu1, sigma1, delta, alpha),
future.seed = TRUE)
mean(results, na.rm = TRUE) # Power: proportion of successes
}Bước 4: Mô phỏng và tính toán power cho mỗi bậc n và kịch bản giả định
# Run simulations across all scenarios and sample sizes
with_progress({
p <- progressor(along = 1:(nrow(scenarios) * length(sample_sizes)))
power_results <- scenarios %>%
expand_grid(n = sample_sizes) %>%
mutate(
power = future_pmap_dbl(
list(n, mu1, sigma1),
~ {
p() # Update progress bar
run_simulation(..1, mu0, sigma0,
..2, ..3,
delta,
alpha, n_sim)
}
)
)
})Đây là kết quả mô phỏng với trung bình LH nhóm can thiệp tăng dần từ 1-1.9; ta thấy giá trị trung bình này càng cao (khoảng cách với nhóm can thiệp, và với biên delta càng lớn), thì cỡ mẫu cần thiết càng cao, có thể lên đến 350 cá thể.
ggplot(power_results, aes(x = n,
y = power,
color = factor(mu1))) +
geom_smooth(se = FALSE, linewidth = 0.8) +
geom_hline(yintercept = 0.8,
linetype = "dashed",
color = "red3") +
scale_x_continuous(breaks = seq(0, 500, by = 50)) +
scale_y_continuous(breaks = seq(0, 1, by = 0.1))+
labs(
title = "Power curves for different mean in treatment group",
x = "Total sample Size (n)",
y = "Power",
color = "Mean LH_td in Treatment group"
) +
theme_minimal()8 Kết luận
Qua bài thực hành này, ta có thể thấy sức mạnh và sự linh hoạt của phương pháp tối ưu hóa power curve bằng mô phỏng để ước lượng cỡ mẫu cần thiết cho bất kì một kế hoạch phân tích phức tạp nào, ngay cả khi ta không có dữ liệu từ nghiên cứu thăm dò.
Tặng Khả Nhi một ly cà phê
https://ko-fi.com/khanhilengoc
