Hiệu chỉnh sai số gộp khi so sánh cặp bằng gói emmeans trong r

Trong nghiên cứu thực nghiệm, việc thực hiện nhiều phép so sánh cặp cùng một lúc là vô cùng phổ biến. Tuy nhiên, nếu chúng ta chỉ sử dụng các kiểm định t thông thường mà không hiệu chỉnh, tỷ lệ sai số loại một trên toàn bộ họ kiểm định sẽ tăng lên rất nhanh. Khi đó, việc sử dụng các giá trị p hiệu chỉnh là bắt buộc để đảm bảo tính tin cậy của các kết luận thống kê. Bài viết này sẽ hướng dẫn cách xử lý vấn đề đa so sánh bằng gói lệnh emmeans trong ngôn ngữ R, đồng thời giải thích bản chất thống kê từ phân phối đơn biến đến đa biến.

Thực Nghiệm Và Mô Hình Phân Tích Phương Sai Một Chiều

Chúng ta sẽ bắt đầu với một bộ dữ liệu thực tế về khả năng diệt cỏ của ba hỗn hợp hóa chất so với nhóm đối chứng không xử lý trên cây cỏ dại thuộc họ cà trong ruộng cà chua. Biến phản hồi là khối lượng của cây cỏ dại trong mỗi chậu, và biến giải thích là các công thức xử lý cỏ.

Để đơn giản, chúng ta giả định các điều kiện cơ bản của mô hình phân tích phương sai đều được thỏa mãn. Hãy cùng xây dựng mô hình tuyến tính và thực hiện phân tích phương sai với đoạn mã dưới đây.

1library(statforbiology)
2library(emmeans)
3library(multcomp)
4dataset <- getAgroData("mixture")
5dataset$Treat <- factor(dataset$Treat)
6model <- lm(Weight ~ Treat, data = dataset)
7anova(model)
8# Ket qua phan tich phuong sai:
9# Response: Weight
10#           Df  Sum Sq Mean Sq F value    Pr(>F)    
11# Treat      3 1089.53  363.18  23.663 2.509e-05 ***
12# Residuals 12  184.18   15.35                      

Kết quả phân tích phương sai cho thấy yếu tố xử lý có ý nghĩa thống kê rất cao. Do đó, chúng ta có thể tiếp tục tiến hành so sánh trung bình giữa các công thức theo từng cặp. Nếu chúng ta so sánh mà không thực hiện bất kỳ hiệu chỉnh nào cho sai số gộp, các giá trị p thu được sẽ tương đương với việc sử dụng hàm phân phối Student đơn biến thông thường trong R.

1groupMeans <- emmeans(model, ~Treat)
2tab <- contrast(groupMeans, method = "pairwise", adjust = "none")
3tab
4# Ket qua so sanh cap khong hieu chinh:
5#  contrast                         estimate   SE df t.ratio p.value
6#  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.1697
7#  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0168
8#  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  <.0001
9#  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0012
10#  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
11#  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0038

Chúng ta có thể dễ dàng kiểm chứng lại các giá trị p này bằng cách tính toán thủ công từ trị tuyệt đối của tỷ số t thông qua hàm tích lũy phân phối Student đơn biến hai đuôi với mười hai độ tự do.

1abst <- abs(as.data.frame(tab)$t.ratio)
22 * pt(abst, 12, lower.tail = FALSE)
3# Ket qua kiem chung:
4# [1] 1.696785e-01 1.683167e-02 3.651239e-05 1.157189e-03 4.782986e-06
5# [6] 3.794451e-03

Chuyển Đổi Từ Phân Phối Đơn Biến Sang Phân Phối Đa Biến

Khi thực hiện đồng thời sáu phép so sánh cặp, xác suất xảy ra ít nhất một kết luận bác bỏ sai giả thuyết không tăng lên đáng kể. Để kiểm soát sai số gộp này, chúng ta cần chuyển đổi mô hình xác suất từ phân phối Student đơn biến sang phân phối Student đa biến.

Hãy xem xét tỷ số t đầu tiên thu được từ kết quả trên là 1.461. Câu hỏi đặt ra là: xác suất để nhận được một tỷ số t cực đoan bằng hoặc hơn mức 1.461 từ một phân phối Student đa biến sáu chiều là bao nhiêu? Trong tính toán này, chúng ta phải lưu ý rằng sáu phép kiểm định này có mối tương quan với nhau vì chúng cùng chung một mẫu số là sai số chuẩn của mô hình. Trong trường hợp dữ liệu cân bằng và đồng nhất phương sai, hệ số tương quan giữa các cặp so sánh này bằng đúng 0.5.

Phương Pháp Hiệu Chỉnh Tukey Cho Dữ Liệu Cân Bằng

Trước đây, khi năng lực tính toán của máy tính còn hạn chế, việc tính trực tiếp xác suất tích lũy từ phân phối đa biến là cực kỳ khó khăn. Để giải quyết, các nhà thống kê đã sử dụng phân phối phạm vi Student hóa, hay còn gọi là phương pháp Tukey. Đây cũng là lựa chọn mặc định của hàm contrast trong gói emmeans.

1tab <- contrast(groupMeans, method = "pairwise")
2tab
3# Ket qua so sanh dung phuong phap Tukey:
4#  contrast                         estimate   SE df t.ratio p.value
5#  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.4885
6#  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0698
7#  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
8#  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0055
9#  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
10#  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0173

Chúng ta cũng có thể tái lập chính xác các giá trị p này bằng cách sử dụng hàm ptukey trong R:

1abst <- abs(as.data.frame(tab)$t.ratio)
2ptukey(sqrt(2) * abst, 4, 12, lower.tail = FALSE)
3# Ket qua tinh tu phan phoi pham vi Student hoa:
4# [1] 4.884620e-01 6.981178e-02 1.853807e-04 5.501451e-03 2.473776e-05
5# [6] 1.725725e-02

Phương pháp Tukey mang lại kết quả chính xác tuyệt đối cho dữ liệu cân bằng và hoạt động rất tốt ngay cả khi dữ liệu có sự mất cân bằng nhẹ. Đây chính là thuật toán nền tảng của các kiểm định Tukey HSD và Tukey-Kramer kinh điển.

Tính Toán Trực Tiếp Từ Phân Phối Student Đa Biến

Với sự phát triển của phần cứng và thuật toán hiện đại, ngày nay chúng ta có thể tính toán trực tiếp xác suất từ phân phối Student đa biến. Phương pháp tổng quát này được hiện thực hóa trong gói lệnh mvtnorm thông qua hàm pmvt.

Để tính toán, chúng ta cần xác định khoảng tích lũy cho từng chiều (trong trường hợp này là khoảng từ âm 1.461081 đến dương 1.461081), số độ tự do và ma trận tương quan của các tổ hợp tuyến tính được trích xuất trực tiếp từ đối tượng đồ thị của emmeans. Đoạn mã dưới đây minh họa quy trình này:

1library(mvtnorm)
2t1 <- abs(as.data.frame(tab)$t.ratio)[1]
3ncontr <- 6
4corMat <- cov2cor(vcov(tab))
5plev <- pmvt(lower = rep(-t1, ncontr), upper=rep(t1, ncontr), df = 12, corr = corMat)[1]
61 - plev
7# Ket qua tinh toan truc tiep:
8# [1] 0.4883843

Trong thực tế sử dụng gói emmeans, chúng ta chỉ cần khai báo tham số hiệu chỉnh mvt để tự động kích hoạt phương pháp này.

1tab <- contrast(groupMeans, method = "pairwise", adjust = "mvt")
2tab
3# Ket qua khi su dung tuy chon mvt:
4#  contrast                         estimate   SE df t.ratio p.value
5#  Metribuzin__348 - Mixture_378        4.05 2.77 12   1.461  0.4885
6#  Metribuzin__348 - Rimsulfuron_30    -7.68 2.77 12  -2.774  0.0698
7#  Metribuzin__348 - Unweeded         -17.60 2.77 12  -6.352  0.0002
8#  Mixture_378 - Rimsulfuron_30       -11.73 2.77 12  -4.235  0.0054
9#  Mixture_378 - Unweeded             -21.64 2.77 12  -7.813  <.0001
10#  Rimsulfuron_30 - Unweeded           -9.91 2.77 12  -3.578  0.0172

Vì hàm số này sử dụng phương pháp tích phân số học dựa trên mô phỏng ngẫu nhiên, kết quả giữa các lần chạy có thể có sai số cực nhỏ và không lặp lại hoàn toàn giống nhau. Tuy nhiên, về mặt tiệm cận, kết quả thu được hoàn toàn tương đương với phương pháp Tukey. Do tính phức tạp trong tính toán, chúng ta không nhất thiết phải dùng hiệu chỉnh mvt cho các thí nghiệm cân bằng thông thường, nhưng đây sẽ là công cụ cực kỳ mạnh mẽ khi đối mặt với các bộ dữ liệu mất cân bằng nghiêm trọng.

✨ Việc áp dụng hiệu chỉnh đa so sánh không chỉ là một yêu cầu kỹ thuật để tránh các phát hiện dương tính giả, mà còn thể hiện tư duy phân tích dữ liệu nghiêm túc. Lựa chọn phương pháp hiệu chỉnh phù hợp với cấu trúc dữ liệu cân bằng hay không cân bằng sẽ giúp tối ưu hóa sức mạnh thống kê mà vẫn kiểm soát chặt chẽ sai số hệ thống của toàn bộ nghiên cứu.

Khi làm việc với các bộ dữ liệu thực tế bị khuyết mẫu nghiêm trọng dẫn đến mất cân bằng lớn giữa các nhóm, việc sử dụng phương pháp hiệu chỉnh Tukey truyền thống có thể dẫn đến những sai lệch nào so với phương pháp phân phối đa biến mvt? Hãy thử nghiệm trên một bộ dữ liệu không cân bằng của riêng bạn và so sánh sự khác biệt giữa hai kết quả này.