1 Ziel der ANOVA mit Messwiederholung
Die ANOVA mit Msswiederholung (auch: einfaktorielle Varianzanalyse mit Messwiederholung) testet drei oder mehr Messungen von abhängigen oder denselben Testsubjekten auf unterschiedliche Mittelwerte.
Die Nullhypothese lautet, dass keine Mittelwertunterschiede (hinsichtlich der Testvariable) über die Zeit existieren. Demzufolge lautet die Alternativhypothese, dass über die Zeit Unterschiede existieren. Es ist das Ziel, die Nullhypothese zu verwerfen und die Alternativhypothese anzunehmen.
Für die ANOVA mit Messwiederholung gibt es auch ein Tutorial in SPSS.
2 Voraussetzungen der einfaktoriellen Varianzanalyse (ANOVA)
Die wichtigsten Voraussetzungen der ANOVA sind:
- drei oder mehr voneinander abhängige Messungen – Messwiederholungen bei denselben Testsubjekten
- metrisch skalierte y-Variable
- Normalverteilung der Testvariable je Messzeitpunkt
- Sphärizität – Homogene (nahezu gleiche) Varianzen der paarweisen Differenzen zwischen jeweils 2 Zeitpunkten
- Idealerweise liegen die Daten im long-Format in R vor.
3 Vorbetrachtungen zur ANOVA mit Messwiederholung in R
3.1 Beispiel und Datensatz
Der Datensatz findet sich am Ende des Artikels. Eine Faktorisierung der Zeitvariable wird empfohlen.Im fiktiven Beispiel prüfe ich die Wirksamkeit eines 3-wöchigen Trainingsplans. Dazu wird vor dem Beginn des Trainingsplans, nach Abschluss des Trainingsplans und weitere 3 Wochen nach Abschluss des Trainingsplans ein Leistungswert ermittelt. Ich vermute dahingehend Unterschiede, dass die Leistungswerte ansteigen, sich aber vermutlich nach Beendingung des Trainings nicht weiterentwickeln, eher im Gegenteil.
3.2 Deskriptive Voranalyse
Nach dem Einlesen der Daten kann direkt ein deskriptiver Vergleich gestartet werden, der im Rahmen der ANOVA mit Messwiederholung aber nicht zwingend notwendig ist. Beim Verstehen der Daten und Schreiben der Ergebnisse hilft es allerdings.
Hierzu nutze ich das Paket “dplyr“, was ich mit “install.packages” installiere und mit library(dplyr) lade. Zusätzliche verwende ich Pipe-Operatoren (%>%).
Ich berechne mit der group_by()-Funktion je Zeitpunkt (“time”) den Mittelwert (mean) und die Standardabweichung (sd) meiner Testvariable (“value”).
install.packages("dplyr")
library(dplyr)
df %>%
group_by(time) %>%
summarize(M = mean(value),
SD = sd(value)) %>%
as.data.frame()
Hier erhält man folgenden Output:
time M SD
1 T0 24.300 4.366
2 T1 27.267 3.603
3 T2 28.417 3.880
Hier ist schon erkennbar, dass die Mittelwerte von T0 zu T1 ansteigen und noch ein wenig von T1 zu T2.
Die Streuung ist in etwa vergleichbar und nimmt über die Zeit nur leicht ab, was auch mittels Boxplots veranschaulicht werden kann.
boxplot(df$value~df$time)
3.3 Prüfung auf Normalverteilung
Die Prüfung auf Normalverteilung bei der ANOVA mit Messwiederholung bezieht sich auf die Testvariable je Messzeitpunkt. Demzufolge sind mehrere Prüfungen vorzunehmen. Typisch für Prüfung auf Normalverteilung ist die Durchführung des Shapiro-Wilk-Tests. Allerdings ist hiervon abzuraten, da dessen Power bei kleinen Stichproben zu gering ist (Fehler 2. Art). Bei großen Stichproben hingegen werden aufgrund hoher Power vernachlässigbare Abweichungen von der Normalverteilung als gravierender eingestuft als sie letztlich sind (Fehler 1. Art).
Hierzu verwende ich die ggqqplot()-Funktion des ggpubr-Pakets. Besonders wichtig ist das Argument facet.by, das die Aufteilung je Zeitpunkt vornimmt.
.
library(ggpubr)
# Format: ggqqplot (Datenquelle, "Testvariable", facet.by = "Zeitvariable")
ggqqplot(df, "value", facet.by = "time")
Zur Lesart:
- Die Linie ist die perfekte Normalverteilung,
- die Punkte sind die Beobachtungen der Testvariable im Datensatz.
- Bei einer perfekten Normalverteilung der Testvariable im Datensatz würden sämtliche Punkte auf der Geraden liegen.
- In der Praxis ist dies nicht zu erwarten und Abweichungen sind die Regel.
- Speziell im jeweils oberen rechten und unteren linken Bereich sind Abweichungen typisch, sollten aber nicht zu ausgeprägt sein.
- Das graue Konfidenzband kann noch als zusätzliche Orientierung dienen, da eine ANOVA auch bei leichten Abweichungen von der Normalverteilung als robsut gilt vgl. Blanca et al. (2017).
3.4 Prüfung auf Sphärizität
Sphärizität heißt, dass die Varianzen der paarweisen Differenzen ähnlich sind. Es werden also zwischen allen Messzeitpunkten Differenzen errechnet. Die Varianzen aller dieser Differenzen sollten ähnlich sein, damit man von Sphärizität sprechen kann. Dies wird erst beim Rechnen der ANOVA mit Messwiederholung ersichtlich, sodass die Prüfung direkt im nächsten Kapitel (mit) erfolgt. Bei Verletzung wird automatisch eine Korrektur vorgenommen.
4 ANOVA mit Messwiederholung in R rechnen und interpretieren
4.1 ANOVA mit Messwiederholung
Hierzu wird die anova_test()-Funktion des rstatix-Pakets verwendet:
- data ist der Data Frame, in dem die Daten liegen
- dv (dependent variable) ist die Testvariable, hier value
- wid (within ID) ist die ID, die zur Zuordnung der Messungen zu Testsubjekten verwendet wird, hier ID
- within ist die Zeitvariable, hier time
- Ich lasse mir die Ergebnisse in den Vektor “rep_anova” speichern und den Mauchly’s Test für Sphärizität ausgeben.
library(rstatix)
rep_anova <- anova_test(data = df, dv = value, wid = ID, within = time)
rep_anova$`Mauchly's Test for Sphericity`
Einschub: Sphärizität prüfen
Effect W p p<.05
1 time 0.967 0.375
Die Nullhypothese von Mauchly's Test geht von Sphärizität aus. Der p-Wert ist hier hinreichend groß. Somit kann von Sphärizität ausgegangen werden. Sollte dies mal nicht der Fall sein, rechnet anova_test() automatisch eine Korrektur ein. Die Ergebnisse können aber analog zu den folgenden Ausführungen interpretiert werden.
Die Ergebnisse der ANOVA mit Messwiederholung fordere ich mit dem Befehl get_anova_table() an:
get_anova_table(rep_anova)
Die Ergebnisse sind übersichtlich gehalten:
ANOVA Table (type III tests)
Effect DFn DFd F p p<.05 ges
1 time 2 118 18.458 1.06e-07 * 0.163
- Das Hauptaugenmerk liegt auf dem p-Wert.
- Ist der p-Wert hinreichend klein (typischerweise unter Alpha = 0,05), kann die Nullhypothese keiner Unterschiede über die Zeit verworfen werden.
- Der p-Wert ist mit 1.06e-07 ( = 0.000000106) sehr klein.
- Die Nullhypothese kann verworfen werden: Es können Unterschiede über die Zeit beobachetet werden.
Die entscheidende Frage ist nun, zwischen welchen der drei Zeitpunkte ein Unterschied existiert. Hierzu braucht es zwingend die nachfolgende post-hoc-Analyse.
4.2 Post-hoc-Analyse paarweise Vergleiche über die Zeit
Diese führt man mittels paarweisen t-Tests ("pairwise_t_test()" aus dem rstatix-Paket) durch. Erneut hilft uns der Pipe-Operator (%>%) dies unkompliziert zu berechnen.
Allerdings muss hierbei der p-Wert angepasst werden, da das mehrfache Testen auf dieselbe Stichprobe zu einem erhöhten Alphafehler führt. Dies hat wiederum zur Folge, dass die Wahrscheinlichkeit einen Fehler 1. Art zu begehen steigt. Ultimativ könnte das dazu führen, dass man die Nullhypothese fälschlicherweise ablehnt, also Unterschiede unterstellt, die nicht existieren.
Aber keine Angst, es existiert das Argument namens "p.adjust.method". Es gibt für p.adjust.method verschiedene Anpassungsmöglichkeiten. Zumeist wählt man die konservativste "bonferroni". Weitere Möglichkeiten sind "holm", "hochberg", "hommel", "BH", "BY", "fdr", "none". Letzteres ist aus o.g. Grund nicht zu empfehlen.
Der Code zum paarweisen Vergleich sowie dem Anpassen des p-Wertes ist Testvariable~Zeitvariable und muss zwingend das Argument "paired = TRUE" enthalten:
library(dplyr)
library(rstatix)
df %>%
pairwise_t_test(value~time, paired = TRUE,
p.adjust.method = "bonferroni") %>%
as.data.frame()
Als Ergebnis erhält man eine kleine Übersichtstabelle mit allen paarweisen t-Tests. Bei 3 Gruppen werden in 3 Tests alle Gruppen jeweils mit den anderen verglichen. Von Interesse ist hier v.a. die Spalte p.adj., welche die adjustierten p-Werte für den jeweiligen Test enthält.
.y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
1 value T0 T1 60 60 -3.999294 59 1.79e-04 3.58e-04 ***
2 value T0 T2 60 60 -5.741029 59 3.46e-07 1.04e-06 ****
3 value T1 T2 60 60 -1.813418 59 7.50e-02 7.50e-02 ns
In der obigen Tabelle kann man folgendes erkennen:
- Der Unterschied zwischen der Zeitpunkt 0 und Zeitpunkt 1 weist eine adjustierte Signifikanz von p = 3.58e-04 (= 0.000358) aus. Für diese beiden Zeitpunkte kann die Nullhypothese keines Unterschiedes demzufolge abgelehnt werden.
- Für den Unterschied zwischen Zeitpunkt 0 und Zeitpunkt 2 ist die adjustierte Signifikanz p = 1.04e-06 (= 0.00000104). Auch hier kann die Nullhypothese keines Unterschiedes verworfen werden. Die Alternativhypothese wird angenommen.
- Für den Unterschied zwischen Zeitpunkt 1 und Zeitpunkt 2 ist allerdings eine adjustierte Signifikanz von p = 7.50e-02 (= 0.075) zu erkennen. Die Nullhypothese (kein Unterschied) kann nicht verworfen werden. Die Alternativhypothese wird ebenfalls angenommen.
- Hinweis: Die Signifikanzindikatoren in der letzten Spalte beziehen sich auf einen Vergleich des adjustierten p-Wertes mit einem Alpha-Wert von 0.05. Das strikte Festhalten an dieser Grenze ist nicht mehr zeitgemäß (vgl. Wasserstein (2016)).
Im Ergebnis kann festgehalten werden, dass sich zwischen Zeitpunkt 0 und 1 eine Steigerung ergeben hat. Zwischen Zeitpunkt 0 und 2 war die Steigerung ebenfalls beobachtbar. Dadurch, das allerdings zwischen Zeitpunkt 1 und 2 kein Unterschied beobachtbar war, scheint der erzielte Effekt nicht (direkt) wieder zu verschwinden - eine weitere Steigerung wäre sachlogisch ohnehin unwahrscheinlich gewesen und im Rahmen der deskriptiven Voranalyse auch schon ausgeschlossen.
4.3 Effektstärke der paarweisen Vergleiche
Sollten mit paarweisen Vergleichen hinreichende Unterschiede über die Zeit beobachtet werden können, sind diese Unterschiede (= Effekte) zu quantifizieren.
Erneut wird mit dem Konzept des Piping gearbeitet, welches mit dem dplyr-Paket zur Anwendung kommt. Gleichzeitig braucht es die cohens_d()-Funktion des rstatix-Pakets. Der Code hierfür ist analog zur pairwise_t_test()-Funktion aus dem vorherigen Abschnitt - abzüglich p.adjust.method:
library(dplyr)
library(rstatix)
df %>%
cohens_d(value~time, paired = TRUE) %>%
as.data.frame()
Dies führt im Ergebnis zu folgendem Output:
.y. group1 group2 effsize n1 n2 magnitude
1 value T0 T1 -0.5163066 60 60 moderate
2 value T0 T2 -0.7411637 60 60 moderate
3 value T1 T2 -0.2341112 60 60 small
Hier sehen wir die Effektstärken in der Spalte "effsize" für alle drei paarweisen Vergleiche. Die Effektstärken sollten nur für paarweise Vergleiche abgelesen werden, für die ein hinreichend kleiner p-Wert existiert.
Im Beispiel war das:
- T0 vs. T1 mit einer Effektstärke von d = 0.52 und
- T0 vs. T2 mit einer Effektstärke von d = 0.74.
- T1 vs. T2 wird wegen des zu hohen p-Wertes ignoriert.
- Effektstärken werden immer als positive Werte berichtet.
- Es gibt Fachdisziplinen die, auch für nicht signifikante Effekte ein Reporting der Effektstärke wünschen.
In der Spalte "magnitude" wird die Größe des Effektes eingeordnet, dies erfolgt auf Basis von Cohen (1992): A Power Primer, S. 157:
- Ab 0,2 ist es ein schwacher Effekt,
- ab 0,5 ein mittlerer und
- ab 0,8 ein starker Effekt.
4.4 Effektstärke der ANOVA
Die Effektstärke f wird von der anova_test()-Funktion mit ausgegeben: das sog. Eta², auch generalized eta square. Allerdings kann es zusätzlich noch in f umgewandelt werden. Wie wir bereits festgestellt haben, interessiert uns zwar eher das Ergebnis der post-hoc-Analyse. Dennoch kann man den f-Wert berechnen, der sich aus Eta² ergibt, wie folgende Formel aus Cohen (1988), S. 284 zeigt:
![\[ f = \sqrt{\frac{\eta^2}{1-\eta^2}} \]](https://bjoernwalther.com/wp-content/ql-cache/quicklatex.com-26b0e4363f8a9a93a22800043d6aac50_l3.png "Rendered by QuickLaTeX.com")
