ANOVA trong R: aov(), bảng ANOVA và Tukey HSD

Vì sao dùng ANOVA chứ không phải một loạt t-test

Một t-test so sánh hai trung bình. Với ba nhóm trở lên, cách làm hấp dẫn là chạy t-test cho mọi cặp, và đó chính là cái bẫy. Mỗi phép kiểm định ở mức 0,05 mang rủi ro dương tính giả 5%, và các rủi ro cộng dồn lại: với 4 nhóm là 6 phép kiểm định theo cặp và khoảng 26% khả năng có ít nhất một kết quả "có ý nghĩa" giả; với 5 nhóm (10 phép kiểm định), khoảng 40%. Bạn sẽ tạo ra những phát hiện từ nhiễu thuần túy.

ANOVA (ANalysis Of VAriance, phân tích phương sai) khắc phục điều này bằng cách đặt một câu hỏi với một p-value: mọi trung bình nhóm đều bằng nhau, hay có ít nhất một nhóm khác? Dù tên gọi như vậy, nó so sánh các trung bình, chỉ là bằng cách phân tích phương sai: nếu trung bình các nhóm phân tán nhiều hơn mức mà nhiễu bên trong nhóm có thể giải thích, thì có điều gì đó thật sự đang diễn ra.

ANOVA một yếu tố với aov()

PlantGrowth sinh ra là để dùng cho việc này: khối lượng cây trong một điều kiện đối chứng và hai điều kiện xử lý. Dựng mô hình bằng aov(), rồi (điều quan trọng) in bảng bằng summary():

fit <- aov(weight ~ group, data = PlantGrowth)
summary(fit)

Đọc công thức là "weight có phụ thuộc vào group không?". Biến phân nhóm phải là một factor: PlantGrowth$group đã là factor sẵn, nhưng nếu các nhóm của bạn được mã hóa bằng số (liều lượng, mã lô), hãy bọc chúng lại: aov(y ~ factor(dose), ...). Nếu không, aov() sẽ lặng lẽ khớp một đường hồi quy qua các mã nhóm thay vì so sánh trung bình nhóm: sai mô hình mà không có thông báo lỗi nào.

Đọc bảng ANOVA

Bảng gồm hai hàng (yếu tố và phần dư) và năm cột. Từng ô một:

            Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159
Residuals   27 10.492  0.3886

• Df: bậc tự do. Hàng yếu tố có số nhóm − 1 (3 nhóm → 2); hàng phần dư có số quan sát − số nhóm (30 − 3 = 27). Một cách kiểm tra nhanh rằng R đã hiểu đúng thiết kế bạn định dùng.
• Sum Sq: biến thiên, chia thành hai phần. Hàng group là phần giữa các nhóm: trung bình các nhóm cách trung bình chung bao xa. Residuals là phần trong nội bộ nhóm: các cây dao động quanh trung bình của chính nhóm mình bao nhiêu. Cộng lại, chúng bằng tổng biến thiên trong dữ liệu.
• Mean Sq: mỗi Sum Sq chia cho Df tương ứng, biến các phần biến thiên thành tỷ lệ trên mỗi bậc tự do để so sánh được. Mean Sq của phần dư (0,389) là mức nhiễu.
• F value: tỷ số giữa Mean Sq của group và Mean Sq của Residuals (1,8832 / 0,3886 ≈ 4,85). Nếu mọi trung bình nhóm thực sự bằng nhau, tỷ số này dao động quanh 1. F càng lớn, chênh lệch giữa các nhóm càng vượt xa mức mà nhiễu có thể giải thích.
• Pr(>F): p-value, xác suất nhận được giá trị F lớn như vậy nếu cả ba trung bình thực đều bằng nhau. Ở đây là 0,016, đủ nhỏ ở mức thông dụng 0,05 để bác bỏ giả thuyết "tất cả bằng nhau". Như mọi khi, đây không phải xác suất giả thuyết không là đúng, và nó không nói gì về việc nhóm nào khác nhau hay khác bao nhiêu.

Điểm cuối cùng đó là giới hạn then chốt: một giá trị F có ý nghĩa chỉ nói "có khác biệt, ở đâu đó". Không hơn.

Nhóm nào khác nhau? TukeyHSD()

Câu hỏi tiếp theo cần một kiểm định hậu nghiệm (post-hoc). Phép kiểm định Honest Significant Difference của Tukey kiểm định mọi cặp, đồng thời giữ tỷ lệ sai theo cả họ so sánh ở mức 5% cho tất cả các phép so sánh cùng lúc:

fit <- aov(weight ~ group, data = PlantGrowth)
TukeyHSD(fit)

Mỗi cặp một hàng, mỗi hàng bốn con số:

• diff: chênh lệch trung bình ước lượng (nhóm được nêu thứ hai trừ nhóm thứ nhất).
• lwr, upr: khoảng tin cậy 95% theo cả họ so sánh cho chênh lệch đó.
• p adj: p-value, đã được hiệu chỉnh cho việc thực hiện ba phép so sánh.

Với PlantGrowth: trt2-trt1 cho chênh lệch khoảng 0,87 với p adj ≈ 0,012 và khoảng tin cậy nằm hẳn ngoài số 0: cây ở điều kiện xử lý 2 lớn hơn cây ở điều kiện xử lý 1. Cả hai hàng so sánh xử lý với đối chứng (trt1-ctrl, trt2-ctrl) đều có khoảng tin cậy vắt qua số 0 và p adj cao hơn hẳn 0,05: trên mẫu này, không điều kiện xử lý nào phân biệt được với nhóm đối chứng. Quy tắc kinh nghiệm giống hệt với khoảng tin cậy ở mọi nơi khác: khoảng tin cậy không chứa số 0 ⇔ p adj dưới 0,05.

Hãy để ý cách giá trị F tổng thể nói "có khác biệt" trong khi Tukey chỉ ra nó nằm ở đúng một cặp: cấu trúc hai bước chính là toàn bộ thiết kế. Một phép kiểm định tổng thể trung thực, rồi đến phần điều tra theo cặp đã được hiệu chỉnh đúng cách.

ANOVA hai yếu tố: hai yếu tố và tương tác giữa chúng

Với hai biến phân nhóm, một lời gọi kiểm định cả hai, cùng với việc chúng có tương tác với nhau không. ToothGrowth kết hợp loại chất bổ sung (supp) với liều lượng (0,5, 1, 2 mg; là số nên cần factor()):

fit2 <- aov(len ~ supp * factor(dose), data = ToothGrowth)
summary(fit2)

supp * factor(dose) mở rộng thành ba hiệu ứng, mỗi hiệu ứng một hàng trong bảng:

• supp: tính trung bình qua các liều, loại chất bổ sung có ảnh hưởng không? (Có: p ≈ 0,0002.)
• factor(dose): tính trung bình qua các chất bổ sung, liều lượng có ảnh hưởng không? (Rất rõ: p cực nhỏ.)
• supp:factor(dose): tương tác: tác động của chất bổ sung có thay đổi tùy theo liều không? Ở đây p ≈ 0,022, tức là có. Nước cam vượt trội axit ascorbic ở liều thấp, nhưng khoảng cách thu hẹp ở liều 2 mg.

Một tương tác có ý nghĩa là nhãn cảnh báo dán lên các hiệu ứng chính: "loại chất bổ sung có ảnh hưởng" chỉ đúng tính trung bình, và con số trung bình che giấu một câu chuyện phụ thuộc vào liều. Khi tương tác có ý nghĩa, hãy mô tả các tổ hợp (bằng TukeyHSD(fit2, "supp:factor(dose)") hoặc một bảng trung bình nhóm qua thống kê theo nhóm) thay vì chỉ báo cáo các hiệu ứng chính. Chỉ dùng + thay cho * khi bạn chủ ý muốn một mô hình không có tương tác.

Các giả định và phương án dựa trên thứ hạng

Phần toán của ANOVA dựa trên ba giả định, xếp từ khó thương lượng nhất đến dễ nhất:

• Tính độc lập: các quan sát không ảnh hưởng lẫn nhau. Không gì sửa được vi phạm sau khi đã xảy ra; đây là đặc tính của thiết kế nghiên cứu.
• Phương sai bằng nhau giữa các nhóm: Mean Sq của phần dư là một ước lượng nhiễu gộp chung, nên các nhóm cần có mức nhiễu xấp xỉ nhau. Kiểm tra trong một dòng: bartlett.test(weight ~ group, data = PlantGrowth) (p-value lớn nghĩa là không có bằng chứng về phương sai không bằng nhau).
• Phần dư xấp xỉ phân phối chuẩn: ít quan trọng nhất khi thiết kế cân bằng và kích thước nhóm vừa đủ; hãy nhìn qua một histogram hoặc boxplot của phần dư qua residuals(fit).

Khi dữ liệu lệch mạnh, là dữ liệu thứ bậc hoặc đầy giá trị ngoại lai, phương án thay thế dựa trên thứ hạng cho ANOVA một yếu tố là kruskal.test(weight ~ group, data = PlantGrowth), người anh em nhiều nhóm của wilcox.test(). Chỉ một dòng, đổi một chút độ mạnh kiểm định lấy sự vững chắc.

Điều cần nhớ

• ANOVA đặt một câu hỏi về trung bình của 3 nhóm trở lên với một p-value: cách khắc phục vấn đề nhiều t-test.
• fit <- aov(y ~ group, data = df); summary(fit), và biến phân nhóm phải là factor.
• Giá trị F là tín hiệu giữa các nhóm chia cho nhiễu trong nội bộ nhóm; Pr(>F) nhỏ nghĩa là "có khác biệt ở đâu đó", không hơn.
• TukeyHSD(fit) tìm ra những cặp nào khác nhau, đã có sẵn hiệu chỉnh cho nhiều phép so sánh.
• a * b dựng mô hình hai yếu tố cùng tương tác giữa chúng; tương tác có ý nghĩa nghĩa là riêng các hiệu ứng chính không kể hết câu chuyện.
• Kiểm tra phương sai bằng nhau (bartlett.test) và tính chuẩn của phần dư; kruskal.test() là phương án dựa trên thứ hạng.

Tiếp theo: từ so sánh trung bình nhóm sang mô hình hóa một mối quan hệ, hồi quy tuyến tính với lm().