ANOVA w R: aov(), tabela ANOVA i test Tukeya HSD

Dlaczego ANOVA, a nie stos testów t

Test t porównuje dwie średnie. Przy trzech lub więcej grupach kusi, żeby wykonać test t dla każdej pary, i to jest właśnie pułapka. Każdy test na poziomie 0.05 niesie 5% ryzyka wyniku fałszywie dodatniego, a ryzyka się kumulują: przy 4 grupach to 6 testów dla par i około 26% szans na co najmniej jeden pozorny "istotny" wynik; przy 5 grupach (10 testów) około 40%. W ten sposób produkujesz odkrycia z czystego szumu.

ANOVA (ANalysis Of VAriance, analiza wariancji) rozwiązuje ten problem, zadając jedno pytanie z jednym p-value: czy wszystkie średnie grup są takie same, czy co najmniej jedna się różni? Wbrew nazwie ANOVA porównuje średnie, tylko robi to przez analizę wariancji: jeśli średnie grup są bardziej rozrzucone, niż może to wyjaśnić szum wewnątrz grup, dzieje się coś prawdziwego.

Jednoczynnikowa ANOVA z aov()

Zbiór PlantGrowth jest do tego stworzony: masy roślin w grupie kontrolnej i w dwóch wariantach eksperymentalnych. Dopasuj model przez aov(), a potem, co ważne, wypisz tabelę przez summary():

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

Czytaj formułę jako "czy weight zależy od group?". Zmienna grupująca musi być czynnikiem. PlantGrowth$group już nim jest, ale jeśli grupy są zakodowane liczbami (dawki, numery partii), owiń je: aov(y ~ factor(dose), ...). W przeciwnym razie aov() po cichu dopasuje prostą regresji przez kody grup zamiast porównać średnie grup. Zły model i żadnego komunikatu o błędzie.

Jak czytać tabelę ANOVA

Tabela ma dwa wiersze (czynnik i reszty) oraz pięć kolumn. Komórka po komórce:

            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: stopnie swobody. Wiersz czynnika dostaje liczba grup − 1 (3 grupy → 2), a wiersz reszt liczba obserwacji − liczba grup (30 − 3 = 27). To szybki test, czy R widzi taki układ badania, jaki zakładasz.
• Sum Sq: zmienność podzielona na dwie części. Wiersz group to część międzygrupowa: jak daleko średnie grup leżą od średniej ogólnej. Residuals to część wewnątrzgrupowa: jak bardzo rośliny różnią się od średniej własnej grupy. Razem dają całkowitą zmienność danych.
• Mean Sq: każde Sum Sq podzielone przez swoje Df, co zamienia części zmienności w porównywalne wartości na jeden stopień swobody. Mean Sq reszt (0.389) to poziom szumu.
• F value: stosunek Mean Sq dla group do Mean Sq dla Residuals (1.8832 / 0.3886 ≈ 4.85). Gdyby wszystkie średnie grup były naprawdę równe, ten stosunek oscylowałby wokół 1. Im większe F, tym bardziej różnice między grupami przewyższają to, co może wyjaśnić szum.
• Pr(>F): p-value, czyli prawdopodobieństwo uzyskania tak dużego F, gdyby wszystkie trzy prawdziwe średnie były identyczne. Tutaj 0.016, a więc przy standardowym poziomie 0.05 wystarczająco mało, by odrzucić hipotezę "wszystkie równe". Jak zawsze nie jest to prawdopodobieństwo, że hipoteza zerowa jest prawdziwa, i nie mówi nic o tym, które grupy się różnią ani o ile.

Ostatni punkt to kluczowe ograniczenie: istotne F mówi tylko "jest jakaś różnica, gdzieś". Nic więcej.

Które grupy się różnią? TukeyHSD()

Pytanie uzupełniające wymaga testu post hoc. Test HSD Tukeya (Honest Significant Difference) sprawdza każdą parę i utrzymuje poziom błędu dla całej rodziny porównań na 5% łącznie:

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

Jeden wiersz na parę, cztery liczby w każdym wierszu:

• diff: oszacowana różnica średnich (druga grupa w nazwie minus pierwsza).
• lwr, upr: 95% przedział ufności dla tej różnicy, liczony dla całej rodziny porównań.
• p adj: p-value już skorygowane o to, że wykonujesz trzy porównania.

Dla PlantGrowth wiersz trt2-trt1 pokazuje różnicę około 0.87 przy p adj ≈ 0.012 i przedziale, który nie obejmuje zera: rośliny z wariantu 2 rosną lepiej niż z wariantu 1. Oba wiersze porównujące wariant z kontrolą (trt1-ctrl, trt2-ctrl) mają przedziały obejmujące zero i p adj wyraźnie powyżej 0.05, więc na tej próbie żadnego wariantu nie da się odróżnić od kontroli. Reguła jest taka sama jak przy każdych przedziałach ufności: przedział nie obejmuje zera ⇔ p adj poniżej 0.05.

Zwróć uwagę, że ogólne F mówi "różnica istnieje", a test Tukeya lokalizuje ją w dokładnie jednej parze. Na tym polega cały pomysł dwóch kroków: jeden uczciwy test ogólny, a potem prawidłowo skorygowane porównania par.

Dwuczynnikowa ANOVA: dwa czynniki i ich interakcja

Przy dwóch zmiennych grupujących jedno wywołanie testuje oba czynniki, a do tego to, czy ze sobą oddziałują. ToothGrowth krzyżuje rodzaj suplementu (supp) z dawką (0.5, 1, 2 mg; to liczby, więc potrzebne jest factor()):

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

supp * factor(dose) rozwija się do trzech efektów, po jednym wierszu tabeli na każdy:

• supp: czy rodzaj suplementu ma znaczenie, uśredniając po dawkach? (Tak: p ≈ 0.0002.)
• factor(dose): czy dawka ma znaczenie, uśredniając po suplementach? (Zdecydowanie: p jest znikome.)
• supp:factor(dose): interakcja, czyli czy efekt suplementu zmienia się w zależności od dawki? Tutaj p ≈ 0.022, więc tak. Sok pomarańczowy wygrywa z kwasem askorbinowym przy niskich dawkach, ale przy dawce 2 mg różnica znika.

Istotna interakcja to ostrzeżenie przy efektach głównych: "rodzaj suplementu ma znaczenie" jest prawdą tylko średnio, a średnia ukrywa historię zależną od dawki. Gdy interakcja jest istotna, opisuj kombinacje (przez TukeyHSD(fit2, "supp:factor(dose)") albo tabelę średnich grup z podsumowań według grup), a nie same efekty główne. Używaj + zamiast * tylko wtedy, gdy celowo chcesz model bez interakcji.

Założenia i zamiennik oparty na rangach

Matematyka ANOVA opiera się na trzech założeniach, od najmniej do najbardziej elastycznego:

• Niezależność: obserwacje nie wpływają na siebie nawzajem. Jej naruszenia nie da się naprawić po fakcie, bo to cecha projektu badania.
• Równe wariancje w grupach: Mean Sq reszt to jedno wspólne oszacowanie szumu, więc grupy powinny być mniej więcej tak samo zaszumione. Sprawdzisz to jedną linijką: bartlett.test(weight ~ group, data = PlantGrowth) (duże p-value oznacza brak dowodów na nierówne wariancje).
• Mniej więcej normalne reszty: to założenie liczy się najmniej przy zrównoważonym układzie i rozsądnych licznościach grup. Obejrzyj histogram albo wykres pudełkowy reszt przez residuals(fit).

Gdy dane są mocno skośne, porządkowe albo pełne wartości odstających, oparta na rangach alternatywa dla jednoczynnikowej ANOVA to kruskal.test(weight ~ group, data = PlantGrowth), czyli wielogrupowy odpowiednik wilcox.test(). To jedna linijka, a w zamian za odporność traci się trochę mocy testu.

Co warto zapamiętać

• ANOVA zadaje jedno pytanie o średnie 3+ grup z jednym p-value, co rozwiązuje problem wielokrotnych testów t.
• fit <- aov(y ~ group, data = df); summary(fit), a zmienna grupująca musi być czynnikiem.
• Wartość F to sygnał międzygrupowy podzielony przez szum wewnątrzgrupowy; małe Pr(>F) oznacza "jakaś różnica gdzieś" i nic więcej.
• TukeyHSD(fit) znajduje, które pary się różnią, z wbudowaną korektą na wielokrotne porównania.
• a * b dopasowuje dwa czynniki i ich interakcję; istotna interakcja oznacza, że same efekty główne nie opowiadają całej historii.
• Sprawdź równość wariancji (bartlett.test) i normalność reszt; kruskal.test() to zamiennik oparty na rangach.

Dalej: od porównywania średnich grup do modelowania zależności, czyli regresja liniowa z lm().