Ước tính cỡ mẫu cho ANOVA đơn biến

1 Đặt vấn đề

Trong bài thứ nhất của project R thực dụng, Nhi muốn bàn về quy trình ước tính cỡ mẫu sử dụng R. Cỡ mẫu là một chủ đề gây khó khăn cho nhiều anh chị em đồng nghiệp: một mặt đây là việc cần thiết khi xét đến thời gian và công sức phải bỏ ra để thu thập dữ liệu, nhưng mặt khác không phải ai cũng hiểu và làm đúng. Từ đó sinh ra việc sao chép một cách vô nghĩa những công thức từ thế hệ này sang thế hệ khác. Đáng tiếc rằng đến năm 2018 hiện nay việc sao chép này vẫn được làm một cách khá tùy tiện và bất hợp lý trong nhiều luận văn.

Để minh họa, Nhi trình bày một thí dụ có thật, trích từ luận án Tiến sỹ Y học có tựa đề “Nghiên cứu nồng độ leptin, adiponectin huyết tương và tỷ lệ leptin/adiponectin trên đối tượng thừa cân-béo phì” của BS.Võ Minh Phương (ĐH Y Dược Huế, năm 2018).

Bạn có thể đọc toàn văn luận án từ link sau đây:

http://hueuni.edu.vn/sdh/attachments/article/1241/NOIDUNGLA.pdf

Trước hết, Nhi muốn nhấn mạnh là bài viết này không nhằm mục đích phê bình cá nhân TS.BS Phương, hơn ai hết Nhi hiểu rõ bạn không có lỗi gì cả. Nhi chỉ muốn đề xuất một cách làm hợp lý hơn, để những đàn em sau này có thể tham khảo khi thực hiện luận văn cho mình và không vấp ngã như người đi trước.

2 Bình luận về cách làm cũ

Trong luận án, phần xác định cỡ mẫu được trình bày ngay đầu chương Đối tượng và phương pháp nghiên cứu (Mục 2.2, trang 52), điều này có nghĩa là tác giả đã thực hiện việc tính cỡ mẫu như một thủ tục hành chính, thậm chí trước khi biết mình muốn nghiên cứu về điều gì và như thế nào. Thật vậy, chưa có mục tiêu nghiên cứu nào được đặt ra sau phần Tổng quan Y văn trong chương 1, nên Cỡ mẫu không gắn kết với một giả thuyết nào cả. Trong khi đó thiết kế nghiên cứu, danh sách những biến số và phương pháp xử lý số liệu đều được trình bày sau phần tính cỡ mẫu. Đây là điểm bất hợp lý thứ nhất.

Tiếp theo, tác giả đã vận dụng công thức chưa phù hợp. Thực ra, ý nghĩa và mục tiêu của công thức quan trọng hơn nội dung của nó. Tác giả chưa thuyết phục được người đọc về câu hỏi: Mục tiêu của việc tính cỡ mẫu là gì ?

Lưu ý rằng xác định cỡ mẫu thực chất là một bài toán ngược, mục tiêu của cỡ mẫu thường là tối ưu hóa một tiêu chí nào đó liên quan đến phẩm chất của suy diễn thống kê, thí dụ kích thước hiệu ứng (Cohen’s d effect-size), sức mạnh thống kê (1-beta), giảm thiểu nguy cơ sai lầm Type I, type II… Trước hết, ta phải xác định mục tiêu, từ đó ta sẽ có phương trình để giải. Cỡ mẫu N chính là nghiệm của phương trình này. Nội dung phương trình này thường có nguồn gốc từ một trị số thống kê chuẩn hóa (Z-score) hoặc Effect-size như Cohen’s , hoặc phức tạp hơn khi sử dụng các hàm pnorm và qnorm trực tiếp cho Power, alpha…

Nếu bỏ qua bước giải thích này, cỡ mẫu và công thức là vô nghĩa, thậm chí nếu bạn dùng đúng công thức. Nội dung được trình bày, bao gồm kết quả cỡ mẫu được ước tính không đảm bảo điều gì cả cho những kiểm định thống kê và kết quả.

Trở lại với công thức, dù nó rất đơn giản nhưng người đọc không hiểu gì cả. Tác giả phát biểu “Công thức ước lượng một giá trị trung bình µ”,đây là cách giải thích rất tối nghĩa, vì như ta thấy, không có tham số µ nào xuất hiện trong nội dung công thức.

Trong công thức xuất hiện Z, đây là một trị số chuẩn hóa - nhưng có rất nhiều ý nghĩa, thí dụ Z có thể chỉ khác biệt trung bình giữa 2 phân nhóm sau chuẩn hóa, khác biệt giữa 2 tỉ lệ, hoặc khác biệt giữa trung bình 1 nhóm và 1 ngưỡng vô hiệu nào đó ?

Giả định một mục tiêu xác định N để mẫu có phân phối chuẩn, và Z là 1 giá trị đặc biệt nào đó, từ Z-score, µ và sigma thì công thức khả dĩ để tính cỡ mẫu là :

\[Z = \frac{x-\mu}{\sigma /\sqrt{n}}\]

Tuy nhiên trong tài liệu, tác giả đã sử dụng một công thức hoàn toàn khác và các tham số như Z và c được giải thích chưa rõ ràng, nhất là c được cho là mức độ chính xác của nghiên cứu, và sau đó được giả định = 0.05; Nhi không thể hiểu những nội dung này.

Sau đó, tác giả ước tính n dựa vào 1 tham số duy nhất là Sigma, với giá trị lấy từ y văn. Chỉ có 1 kết quả cỡ mẫu duy nhất.

3 Tính cỡ mẫu cho post-hoc test ANOVA đơn biến

Trong bài này, Nhi gợi ý một phương pháp ước tính cỡ mẫu khác hợp lý hơn.

Đầu tiên, mục tiêu của việc ước tính cỡ mẫu được nhắm vào 1 quy trình kiểm định / suy diễn mang tính trọng tâm của nghiên cứu. Sau khi đọc qua luận văn của TS.Bs. Phương, nhất là phần bàn luận thì Nhi nhận thấy một trong những mục tiêu chính trong nghiên cứu này là so sánh nồng độ leptin huyết tương giữa 3 phân nhóm độc lập: bệnh nhân béo phì, người thừa cân và những người bình thường (nhóm chứng).

Mục tiêu này dẫn chúng ta đến phương pháp suy diễn thống kê là post-hoc test trong thiết kế ANOVA đơn biến hay so sánh bắt cặp hàng loạt nhiều phân nhóm.

Trường hợp tổng quát, ta có k phân nhóm, tức là k(k-1)/2 cặp so sánh có thể, mỗi cặp so sánh như vậy tương ứng với một giả thuyết vô hiệu có nội dung: trung bình phân nhóm A = trung bình phân nhóm B

\[H_0:\mu_A=\mu_B\]

Hoặc giả thuyết cần chứng minh là:

\[H_1:\mu_A\lt\mu_B\]

Với µA và µB tương ứng với giá trị trung bình của 2 phân nhóm trong cặp so sánh

Trong trường hợp thiết kế ANOVA không cân xứng, cỡ mẫu các phân nhóm có thể không đồng nhất, do đó ta đặt thêm 1 tham số nữa là kappa chỉ tỉ lệ giữa cỡ mẫu 2 phân nhóm A và B:

\[n_B=\kappa n_A\]

Áp dụng công thức tính cỡ mẫu cho mục tiêu so sánh 2 giá trị trung bình có hiệu chỉnh cho việc bắt cặp nhiều lần, p_value 1 bên từ tài liệu : Chow S, Shao J, Wang H. 2008. Sample Size Calculations in Clinical Research. 2nd Ed. Chapman & Hall/CRC Biostatistics Series, trang 59, ta có

\[n_A=\left(\sigma_A^2+\sigma_B^2/\kappa\right)\left(\frac{z_{1-\alpha/\tau}+z_{1-\beta}}{\mu_A-\mu_B}\right)^2\]

Trong công thức này:

1. Sigma A,B là độ lệch chuẩn cho từng phân nhóm

2. MuA và MuB lần lượt là trung bình Leptin ở nhóm Chứng và nhóm Bệnh lý (MuA được phán đoán là nhỏ hơn muB)

3. sau này khi viết hàm ta sẽ thấy Z được tính bằng hàm qnorm (xác định phân vị của phân bố chuẩn)

4. Alpha là sai lầm Type I (thường được ấn định = 0.05)

5. Beta là sai lầm Type II, lưu ý: Power = 1- beta,

6. Tau là số cặp so sánh cần thực hiện, một hình thức hiệu chỉnh cho giá trị alpha

Như vậy để ước tính cỡ mẫu nA = nB, ta sẽ giải phương trình với nA là nghiệm số, cùng các tham số sigma, kappa, alpha, beta và z-score là giá trị chuẩn hóa khác biệt giữa muA và muB .

Chúng ta viết 1 hàm trong R có công dụng tính cỡ mẫu cho nA và nB

sample_estime=function(muA=NULL,sdA=NULL,
                       muB=NULL,sdB=NULL,
                       kappa=2,
                       tau=3,
                       alpha=0.05,
                       beta=0.2){
  nA=(sdA^2+sdB^2/kappa)*((qnorm(1-alpha/tau)+qnorm(1-beta))/(muA-muB))^2
  na=ceiling(nA)
  nb=na/kappa
  return(rbind(na,nb))
}

Nhi giả định giá trị của các tham số MuA và MuB, sdA và sdB tương ứng với trung bình và sd của Leptin ở phân nhóm người bình thường và bệnh nhân béo phì , ta có thể lấy thông tin này từ y văn, gọi là prior hay giả thuyết tiền định.

Giả sử một em sinh viên vào năm 2019 làm một nghiên cứu tương tự, em phải sử dụng thông tin từ luận án của TS BS Phương làm giả thuyết tiền định khi tính cỡ mẫu.

Trong luận văn này, Bs. Phương đã báo cáo kết quả như sau:

1. Nhóm bình thường: Leptin = 6.75 ± 5.17 µL/L

2. Nhóm thừa cân: Leptin = 9.74 ± 5.76

3. Nhóm béo phì: Leptin = 10.74 ± 5.61

Ta thi hành hàm tính cỡ mẫu cho cặp so sánh; Bình thường vs Béo phì, giả định cỡ mẫu đồng nhất: kappa=1

sample_estime(muA=6.75,
              muB=10.74,
              kappa=1,
              sdA=5.17,
              sdB=5.61,
              tau=3,
              alpha=0.05,
              beta=0.20)

##    [,1]
## na   33
## nb   33

Kết quả cho thấy mỗi nhóm cần có n=33 cho cặp so sánh này

Ta lại thử trên cặp so sánh thứ hai: Bình thường vs Thừa cân với tỉ lệ kappa=2

sample_estime(muA=6.75,
              muB=9.74,
              kappa=2,
              sdA=5.17,
              sdB=5.76,
              tau=3,
              alpha=0.05,
              beta=0.20)

##    [,1]
## na 43.0
## nb 21.5

Như vậy, ta có thể chọn n=43 nếu mục tiêu là so sánh nhóm chứng với cả 2 nhóm còn lại.

4 Mô phỏng nhiều giả thuyết

Quy trình tính cỡ mẫu không dừng ở đó. Do thông tin về tham số ta đưa vào công thức chỉ là giả định, nên ta không thể tự thỏa mãn quá dễ dàng với chỉ 1 giả định duy nhất. Quan trọng nhất là phải tính đường thoát thân sau này, khi trong quá trình thực hiện nghiên cứu vì lý do khách quan ta không thể đảm bảo đủ 60 bệnh nhân nhưng chỉ có 40 chẳng hạn, lúc này ta phải quay ngược lại biện luận để chứng minh 40 là đủ.

Chi bằng ngay từ ban đầu ta chuẩn bị cho tất cả mọi kịch bản có thể xảy ra : chúng ta cần mô phỏng hàng loạt giả định, thí dụ:

1. Ta có thể giả định: Giá trị Leptin trung bình ở nhóm chứng giữ nguyên, nhưng ở nhóm Béo phì có thể dao động từ 8 đến 20 đơn vị, lúc đó cỡ mẫu sẽ thay đổi ra sao ?

library(tidyverse)

sim_df=data_frame(MuB=seq(8,20,0.5),nA=rep(NA,25))

for (i in (1:length(sim_df$MuB))){
  sample_estime(muA=6.75,
                muB=sim_df$MuB[i],
                kappa=1,
                sdA=5,
                sdB=5,
                tau=3,
                alpha=0.05,
                beta=0.20)->out
  sim_df$nA[i]=out[1]
}

ggplot(sim_df,aes(x=MuB,y=nA))+geom_line()+geom_point(size=2,col="red")+theme_bw()

Theo kết quả này, giá trị trung bình của Leptin trong nhóm Béo phì càng thấp thì cỡ mẫu cần thiết càng cao, thậm chí có thể hơn 100, thậm chí 200 nếu MuB < 10. Dĩ nhiên đây chỉ là giả định, nhưng trong y văn đúng là có những báo cáo cho thấy Leptin ở người béo phì dao động trong khoảng này.

2. Trong một kịch bản khác: Giả sử nồng độ Leptin ở người bình thường dao động từ 4-7 đơn vị, trong khi đó giá trị trung bình này dao động từ 8-20 đơn vị ở người béo phì, cỡ mẫu sẽ thay đổi ra sao ?

MuA_sim=c(4.5,5,5.5,6)

sim_df=data_frame(MuB=rep(seq(8,20,0.5),4),nA=rep(NA,100),MuA=c(rep(4,25),rep(5,25),rep(6,25),rep(7,25)))

for (i in (1:length(sim_df$MuB))){
  sample_estime(muA=sim_df$MuA[i],
                muB=sim_df$MuB[i],
                kappa=1,
                sdA=5,
                sdB=5,
                tau=2,
                alpha=0.05,
                beta=0.20)->out
  sim_df$nA[i]=out[1]
}

sim_df$MuA=factor(sim_df$MuA)

ggplot(sim_df,aes(x=MuB,y=nA))+
  geom_line(aes(col=MuA))+
  geom_point(size=2,aes(col=MuA))+
  theme_bw()

Kết quả mô phỏng cho thấy nếu giá trị Leptin ở nhóm chứng càng cao thì cỡ mẫu càng cao, nhưng nếu giá trị Leptin ở người béo phì > 10 đơn vị thì cỡ mẫu không thể vượt quá 50, bất kể giá trị nào của nhóm chứng trong khoảng 4-7

3. Cuối cùng, ta có thể mô phỏng cả ảnh hưởng của sự thay đổi alpha và beta (hay 1-beta) lên cỡ mẫu. Thí dụ ta tăng ngưỡng alpha lên cao dần, từ 0.05 đến 0.01, 0.005 hay thậm chí 0.0001

Việc mô phỏng này có thể quan trọng, nếu vì sự hạn chế về tính khả thi mà ta buộc phải hy sinh kì vọng ban đầu của mình, thí dụ Power không thể tới 0.8 nhưng chỉ là 0.75 ? Hoặc ta kì vọng rất cao (Power =0.9, alpha=0.01)và muốn kiểm tra điều này có khả thi hay không ?

sim_df=data_frame(alpha=rep(c(0.001,0.005,0.01,0.05),5),
                  beta=rep(c(0.05,0.1,0.2,0.25,0.3),4),
                  nA=rep(NA,20))

for (i in (1:length(sim_df$alpha))){
  sample_estime(muA=6.75,
                muB=10.74,
                kappa=1,
                sdA=5.17,
                sdB=5.61,
                tau=3,
                alpha=sim_df$alpha[i],
                beta=sim_df$beta[i])->out
  sim_df$nA[i]=out[1]
}

ggplot(sim_df,aes(x=1-beta,y=nA))+
  geom_line(aes(col=factor(alpha)))+
  geom_point(size=2,aes(col=factor(alpha)))+
  theme_bw()

Kết quả mô phỏng cho thấy alpha càng thấp và/hoặc 1-beta càng cao thì cỡ mẫu yêu cầu càng lớn. Tuy nhiên nếu ta hài lòng với Power = 0.8 thì thậm chí tăng ngưỡng alpha lên 0.001 cũng không đưa đến một cỡ mẫu bất khả thi.

5 Diễn đạt

Quy trình tính toán, mô phỏng nêu trên có thể được diễn đạt trong văn bản như sau:

1. Cơ sở : Do mục tiêu chính của nghiên cứu là khảo sát sự khác biệt của Leptin giữa 3 phân nhóm: Bình thường/Thừa cân và Béo phì, việc ước tính cỡ mẫu dựa trên cơ sở đảm bảo tối ưu sức mạnh thống kê của phép so sánh lặp lại trên 3 cặp so sánh với giả thuyết vô hiệu là trung bình Leptin như nhau.

2. Cách thức: Chúng tôi sử dụng công thức ước tính cỡ mẫu:

\[n_A=\left(\sigma_A^2+\sigma_B^2/\kappa\right)\left(\frac{z_{1-\alpha/\tau}+z_{1-\beta}}{\mu_A-\mu_B}\right)^2\]

Trong đó nA là cỡ mẫu của nhóm chứng trong cặp so sánh A/B, µA,sigmaA và µB ,sigmaB lần lượt là trung bình và độ lệch chuẩn của nồng độ Leptin, Z là trị số thống kê chuẩn hóa được tính theo công thức:

(qnorm(1-alpha/tau)+qnorm(1-beta)

Tham số alpha tương ứng với sai lầm Type I, Beta là sai lầm Type II và sức mạnh thống kê Power =1-beta.

3. Dựa theo thông tin từ nghiên cứu của Võ Minh Phương và cs., chúng tôi thực hiện một số mô phỏng, sau đó đối chiếu với quỹ thời gian và nguồn lực của dự án, và đi đến kết luận:

Nhóm chứng cần 50 cá thể, 2 phân nhóm Thừa cân và Béo phì cần ít nhất 30 bệnh nhân để đảm bảo phát hiện được khác biệt giữa bệnh nhân và nhóm chứng với sức mạnh thống kê 80% và nguy cơ sai lầm Type I< 0.05.

6 Nhận xét

Quy trình xác định cỡ mẫu trong đề cương/luận văn cần được trình bày một cách hợp lý hơn với những nội dung như sau:

1. Nên dựa trên cơ sở một câu hỏi nghiên cứu cụ thể, từ đó phát biểu giả thuyết vô hiệu (so sánh hơn/kém/tương đương), và dẫn đến một công thức tính toán rõ nghĩa.

2. Giá trị của các tham số cần được biện luận sử dụng thông tin từ y văn, giả thuyết tiền định, nghiên cứu pilot hoặc kết quả mô phỏng.

3. Không chỉ trình bày một giá trị duy nhất, nhưng nên có một hay nhiều kết quả mô phỏng cho các kịch bản khác nhau, bao gồm khả năng dao động ngẫu nhiên của các tham số, tỉ lệ mất mát dữ liệu.

4. Cuối cùng, kết quả mô phỏng và giả định sẽ được kiểm tra về tính khả thi, như thời gian, kinh phí, nhân lực và phương tiện… từ đó quyết định một con số cỡ mẫu tối ưu cân bằng giữa kì vọng và tính khả thi.