ANOVA in R: aov(), die ANOVA-Tabelle und Tukey HSD

Warum ANOVA und kein Haufen t-Tests

Ein t-Test vergleicht zwei Mittelwerte. Bei drei oder mehr Gruppen ist der verlockende Zug, jedes Paar per t-Test zu prüfen - und genau das ist die Falle. Jeder auf dem 0,05-Niveau durchgeführte Test trägt ein Falsch-positiv-Risiko von 5 %, und die Risiken summieren sich: Bei 4 Gruppen sind das 6 paarweise Tests und rund 26 % Chance auf mindestens ein falsches "signifikantes" Ergebnis; bei 5 Gruppen (10 Tests) rund 40 %. Du würdest Entdeckungen aus reinem Rauschen fabrizieren.

Die ANOVA (ANalysis Of VAriance) behebt das, indem sie eine Frage mit einem p-Wert stellt: Sind alle Gruppenmittelwerte gleich, oder weicht mindestens einer ab? Trotz des Namens vergleicht sie Mittelwerte - sie tut das nur, indem sie Varianz analysiert: Streuen die Gruppendurchschnitte stärker, als das Rauschen innerhalb der Gruppen erklären kann, ist etwas Reales im Spiel.

Einfaktorielle ANOVA mit aov()

PlantGrowth ist dafür gemacht: Pflanzengewichte unter einer Kontrolle und zwei Behandlungsbedingungen. Passe mit aov() an und gib dann - wichtig - die Tabelle mit summary() aus:

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

Lies die Formel als "hängt weight von group ab?" Die Gruppierungsvariable muss ein Faktor sein - PlantGrowth$group ist bereits einer, aber wenn deine Gruppen als Zahlen codiert sind (Dosen, Chargen-IDs), pack sie ein: aov(y ~ factor(dose), ...). Sonst passt aov() still eine Regressionsgerade durch die Gruppencodes an, statt Gruppenmittelwerte zu vergleichen - falsches Modell, keine Fehlermeldung.

Die ANOVA-Tabelle lesen

Die Tabelle hat zwei Zeilen - den Faktor und die Residuen - und fünf Spalten. Zelle für Zelle:

            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 - Freiheitsgrade. Die Faktorzeile bekommt Gruppen − 1 (3 Gruppen → 2); die Residuenzeile bekommt Beobachtungen − Gruppen (30 − 3 = 27). Eine schnelle Plausibilitätsprüfung, dass R das beabsichtigte Design gesehen hat.
• Sum Sq - Variation, aufgeteilt in zwei Haufen. Die group-Zeile ist der Zwischen-Gruppen-Haufen: wie weit die Gruppenmittelwerte vom Gesamtmittelwert entfernt liegen. Residuals ist der Innerhalb-Gruppen-Haufen: wie stark Pflanzen um den Mittelwert ihrer eigenen Gruppe streuen. Zusammen ergeben sie die Gesamtvariation der Daten.
• Mean Sq - jede Sum Sq geteilt durch ihre Df, was Variationshaufen in vergleichbare Raten pro Freiheitsgrad verwandelt. Die Residuen-Mean-Sq (0,389) ist der Rauschpegel.
• F value - das Verhältnis: Mean Sq von group über Mean Sq von Residuals (1,8832 / 0,3886 ≈ 4,85). Wären alle Gruppenmittelwerte wirklich gleich, schwankt dieses Verhältnis um 1. Je größer F wird, desto stärker übertreffen die Unterschiede zwischen den Gruppen das, was Rauschen erklären kann.
• Pr(>F) - der p-Wert: die Wahrscheinlichkeit eines so großen F, wenn alle drei wahren Mittelwerte identisch wären. Hier 0,016 - klein genug auf dem konventionellen 0,05-Niveau, um "alle gleich" zu verwerfen. Wie immer ist das nicht die Wahrscheinlichkeit, dass die Nullhypothese stimmt, und es sagt nichts darüber, welche Gruppen sich um wie viel unterscheiden.

Dieser letzte Punkt ist die entscheidende Einschränkung: Ein signifikantes F sagt "irgendein Unterschied, irgendwo". Nicht mehr.

Welche Gruppen unterscheiden sich? TukeyHSD()

Die Anschlussfrage braucht einen Post-hoc-Test. Tukeys Honest Significant Difference prüft jedes Paar und hält dabei die familienbezogene Fehlerrate über alle Vergleiche zusammen bei 5 %:

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

Eine Zeile pro Paar, vier Zahlen pro Zeile:

• diff - die geschätzte Mittelwertdifferenz (zweitgenannte Gruppe minus erste).
• lwr, upr - das 95-%-familienbezogene Konfidenzintervall für diese Differenz.
• p adj - der p-Wert, bereits für drei Vergleiche angepasst.

Für PlantGrowth: trt2-trt1 zeigt eine Differenz von etwa 0,87 mit p adj ≈ 0,012 und einem Intervall, das die Null nicht enthält - Behandlung 2 wächst stärker als Behandlung 1. Beide Behandlung-vs-Kontrolle-Zeilen (trt1-ctrl, trt2-ctrl) haben Intervalle, die die Null überspannen, und p adj deutlich über 0,05 - keine der Behandlungen ist in dieser Stichprobe von der Kontrolle unterscheidbar. Die Faustregel spiegelt die von Konfidenzintervallen überall: Intervall schließt die Null aus ⇔ p adj unter 0,05.

Beachte, wie das übergreifende F sagte "ein Unterschied existiert", während Tukey ihn in genau einem Paar verortet - die Zweischrittstruktur ist der ganze Entwurf: ein ehrlicher Gesamttest, dann sauber korrigierte paarweise Detektivarbeit.

Zweifaktorielle ANOVA: zwei Faktoren und ihre Wechselwirkung

Bei zwei Gruppierungsvariablen prüft ein Aufruf beide - dazu, ob sie wechselwirken. ToothGrowth kreuzt den Präparattyp (supp) mit der Dosis (0,5, 1, 2 mg - numerisch, braucht also factor()):

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

supp * factor(dose) entfaltet sich zu drei Effekten, je einer Tabellenzeile:

• supp - spielt der Präparattyp, über die Dosen gemittelt, eine Rolle? (Ja: p ≈ 0,0002.)
• factor(dose) - spielt die Dosis, über die Präparate gemittelt, eine Rolle? (Nachdrücklich: p ist winzig.)
• supp:factor(dose) - die Wechselwirkung: Ändert sich der Effekt des Präparats je nach Dosis? Hier p ≈ 0,022 - er tut es. Orangensaft schlägt Ascorbinsäure bei niedrigen Dosen, aber die Lücke schließt sich bei der 2-mg-Dosis.

Eine signifikante Wechselwirkung ist ein Warnhinweis auf den Haupteffekten: "Der Präparattyp spielt eine Rolle" stimmt nur im Mittel, und der Mittelwert verbirgt eine dosisabhängige Geschichte. Ist die Wechselwirkung signifikant, beschreibe die Kombinationen (ein TukeyHSD(fit2, "supp:factor(dose)") oder eine Tabelle der Gruppenmittelwerte über gruppenweise Zusammenfassungen), statt allein die Haupteffekte zu berichten. Nutze + statt * nur, wenn du bewusst ein Modell ohne Wechselwirkung willst.

Annahmen - und der rangbasierte Rückfallplan

Die Mathematik der ANOVA stützt sich auf drei Annahmen, in absteigender Verhandelbarkeit:

• Unabhängigkeit - Beobachtungen beeinflussen einander nicht. Eine Verletzung lässt sich nachträglich durch nichts beheben; sie ist eine Eigenschaft des Studiendesigns.
• Gleiche Varianzen über die Gruppen - die Residuen-Mean-Sq ist eine gepoolte Rauschschätzung, also sollten die Gruppen etwa gleich verrauscht sein. Einzeilige Prüfung: bartlett.test(weight ~ group, data = PlantGrowth) (ein großer p-Wert bedeutet keine Evidenz für ungleiche Varianzen).
• Grob normale Residuen - zählt bei balancierten Designs und ordentlichen Gruppengrößen am wenigsten; beäug ein Histogramm oder einen Boxplot der Residuen über residuals(fit).

Sind die Daten stark schief, ordinal oder von Ausreißern durchsetzt, ist die rangbasierte Alternative zur einfaktoriellen ANOVA kruskal.test(weight ~ group, data = PlantGrowth) - das Mehrgruppen-Geschwister von wilcox.test(). Es ist eine Zeile und tauscht etwas Teststärke gegen Robustheit.

Was du mitnimmst

• Die ANOVA stellt eine Frage zu 3+ Gruppenmittelwerten mit einem p-Wert - die Lösung für die Multiplizität von t-Tests.
• fit <- aov(y ~ group, data = df); summary(fit) - und die Gruppierungsvariable muss ein Faktor sein.
• Der F-Wert ist Signal zwischen den Gruppen über Rauschen innerhalb der Gruppen; ein kleines Pr(>F) bedeutet "irgendein Unterschied irgendwo", nicht mehr.
• TukeyHSD(fit) findet, welche Paare sich unterscheiden, mit eingebauter Korrektur für multiple Vergleiche.
• a * b passt zwei Faktoren plus ihre Wechselwirkung an; eine signifikante Wechselwirkung bedeutet, dass die Haupteffekte die Geschichte nicht allein erzählen.
• Prüfe gleiche Varianzen (bartlett.test) und die Normalität der Residuen; kruskal.test() ist der rangbasierte Rückfallplan.

Als Nächstes: vom Vergleich von Gruppenmittelwerten zum Modellieren einer Beziehung - lineare Regression mit lm().