Hintergrund: Wir untersuchen die Konzentrationsfähigkeit mittels Reaktionstests von 10 Probanden zu drei verschiedenen Tageszeiten: 1. morgens (08:00 Uhr) 2. mittags (13:00 Uhr) 3. abends (20:00 Uhr).
Da dieselben Personen zu allen drei Zeitpunkten getestet wurden, handelt es sich um verbundene Stichproben (Messwiederholung).
Frage: Unterscheidet sich die Reaktionszeit (in ms) signifikant je nach Tageszeit?
Hinweis: Da wir hier mehr als zwei Gruppen vergleichen, dürfen wir nicht einfach viele t-Tests rechnen (Fehlerkumulierung! [Link]). Wir nutzen die Varianzanalyse mit Messwiederholung (rmANOVA).
Wir folgen dem Workflow der Datenanalyse: [Link]
Bezug Glossar: Variablen & Datentypen [Link]
Der Datensatz den wir analysieren wollen hat den Namen reaktion_daten. Dieser wurde für unser Lern-Szenario zunächst künstlich erstellt.
Wir prüfen die Struktur:
# Struktur prüfen: Welche Datentypen liegen vor?
str(reaktion_daten)'data.frame': 39 obs. of 3 variables:
$ ID : Factor w/ 13 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
$ Zeit : Factor w/ 3 levels "morgens","mittags",..: 1 2 3 1 2 3 1 2 3 1 ...
$ Reaktion: num 350 450 400 340 460 390 360 440 410 340 ...# Welche Gruppen gibt es genau?
levels(reaktion_daten$Zeit)[1] "morgens" "mittags" "abends" Struktur: Wir haben eine metrische Zielvariable (Reaktion) und zwei kategoriale Faktoren ID (die Versuchspersonen) und Zeit (der Messzeitpunkt) mit 3 Stufen (> 2). Das ist der klassische Fall für eine einfaktorielle ANOVA mit Messwiederholung.
Bezug Glossar: Datenqualität sichern [Link]
Wir prüfen zunächst auf fehlende Werte (NA). Bei der Untersuchung von Daten die auf die Zuverlässigkeit der Probanden angewiesen sind, kommt es schnell zu Messlücken, die für nachfolgende Analysen problematisch sein könnten.
# Gibt es fehlende Werte in der Spalte Reaktion?
table(is.na(reaktion_daten$Reaktion))
FALSE TRUE
35 4 Der Datensatz enthält 4 Fehlwerte. Diese können im ungünstigsten Fall auf 4 verschiedene Personen verteilt sein oder einzelne Personen haben mehr als einen Fehlwert. Dies identifizieren wir wie folgt:
NA_IDs <- unique(reaktion_daten$ID[is.na(reaktion_daten$Reaktion)])
NA_IDs[1] 4 9 13
Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13Wir erfahren, dass die 3 Probanden mit der ID 4,9 und 13 mindestens einen Fehlerwert haben. Da wir hier Daten mit Messwiederholung analysieren wäre es falsch nur die einzelnen NA’s zu löschen. Das würde das Zusammenhangsmuster der Daten zerstören. Das einzelne löschen ist nur im unverbundenen Fall möglich. Hier benötigen wir sogenannte Complete Cases und müssen entsprechend alle Messwerte eines Probanden löschen, wenn er mindestens einen Fehlwert hat. Diese Listwise Deletion geht einfach wie folgt:
reaktion_clean <- subset(reaktion_daten, !(ID %in% NA_IDs))Durch das Löschen kompletter Probanden dünnen wir den Datensatz aus, was zu einem Verlust an Power [Link] führen kann. Eine Alternative wäre die Imputation [Link] (das künstliche Auffüllen fehlender Werte). Dies ist jedoch statistisch komplex und kann bei falscher Anwendung Ergebnisse verfälschen, weshalb wir hier beim Ausschluss bleiben.
Abschließend verschaffen wir uns noch einen statistischen Überblick mit der summary- Funktion.
summary(reaktion_clean) ID Zeit Reaktion
1 : 3 morgens:10 Min. :325.0
2 : 3 mittags:10 1st Qu.:356.2
3 : 3 abends :10 Median :402.5
5 : 3 Mean :402.7
6 : 3 3rd Qu.:448.8
7 : 3 Max. :475.0
(Other):12 Wir sehen, dass wir für jede der 10 Personen (ID) exakt 3 Beobachtungen (Zeit) für die Analyse vorliegen haben. Das Design ist also “balanciert”.
Bezug Glossar: Grafische Darstellung
Bevor wir zum testen kommen, schauen wir uns die Verteilungen an. Ein Boxplot [Link] eignet sich perfekt für den Vergleich.
# Boxplot: Reaktion (y) abhängig von der Zeit (x)
boxplot(Reaktion ~ Zeit, data = reaktion_daten,
main = "Reaktionszeit nach Tageszeit",
ylab = "Reaktionszeit (ms)",
col = c("lightblue", "orange", "lightgreen"))
Visuelle Interpretation: Die Reaktionszeiten scheinen mittags deutlich höher (also langsamer) zu sein als morgens. Abends liegen die Werte dazwischen. Die Boxen sind relativ kompakt, was auf geringe Streuung hindeutet. Die Frage ist: Sind diese Unterschiede auch groß genug, um statistisch signifikant zu sein?
Bezug Glossar: Verteilungs-Checks
Die einfaktorielle Varianzanalyse mit Messwiederholungen [Link] hat zwei Haupt-Voraussetzungen: 1. Normalverteilung der Testvariable je Messzeitpunkt. 2. Sphärizität: Die Varianzen der Differenzen zwischen allen Paaren von Zeitpunkten müssen gleich sein. Also eine deutlich komplexere Version der Varianzhomogenität, die bei dem ANOVA-Pendant ohne Messwiederholung vorausgesetzt wird.
Zum Check für die Normalverteilung nutzen wir hier visuell den QQ-Plot [Link]. Hinweis: Wir nutzen das Paket ggpubr für schönere Diagnostik-Plots.
library(ggpubr)Lade nötiges Paket: ggplot2# QQ-Plot für jeden Zeitpunkt einzeln
ggqqplot(reaktion_daten, "Reaktion", facet.by = "Zeit")Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_qq()`).Warning: Removed 4 rows containing non-finite outside the scale range
(`stat_qq_line()`).
Removed 4 rows containing non-finite outside the scale range
(`stat_qq_line()`).
Interpretation: Die Punkte liegen in allen drei Gruppen schön auf der Referenzlinie und innerhalb des grauen Konfidenzbereichs. Das spricht für eine Normalverteilung. Wir gehen davon aus, dass Voraussetzung 1 erfüllt ist. Um diesen visuellen Eindruck statistisch abzusichern, führen wir für jeden Zeitpunkt einen Shapiro-Wilk-Test durch.
shapiro.test(subset(reaktion_clean$Reaktion, reaktion_clean$Zeit == "morgens"))
Shapiro-Wilk normality test
data: subset(reaktion_clean$Reaktion, reaktion_clean$Zeit == "morgens")
W = 0.97455, p-value = 0.9295shapiro.test(subset(reaktion_clean$Reaktion, reaktion_clean$Zeit == "mittags"))
Shapiro-Wilk normality test
data: subset(reaktion_clean$Reaktion, reaktion_clean$Zeit == "mittags")
W = 0.97798, p-value = 0.9535shapiro.test(subset(reaktion_clean$Reaktion, reaktion_clean$Zeit == "abends"))
Shapiro-Wilk normality test
data: subset(reaktion_clean$Reaktion, reaktion_clean$Zeit == "abends")
W = 0.94219, p-value = 0.5776Das visuelle Ergebnis wird hiermit klar bestätigt.
Der Test auf Sphärizität [Mauchly-Test] erfolgt bei der rmANOVA in R meist direkt im Rahmen der Modellberechnung. Wir schauen uns das Ergebnis gleich im Output an.
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Wir berechnen die ANOVA mit Messwiederholung [Link], da die Messungen von denselben Personen stammen. [Link: Unabhängig vs. Verbunden].
Die Nullhypothese der ANOVA lautet: Die mittleren Reaktionszeiten sind zu allen drei Zeitpunkten gleich.
Wir nutzen die Funktion anova_test aus dem Paket rstatix, da sie uns automatisch den Mauchly-Test und eventuelle Korrekturen liefert.
library(rstatix)
Attache Paket: 'rstatix'Das folgende Objekt ist maskiert 'package:stats':
filter# Berechnung der rmANOVA
# dv = Dependent Variable (Reaktion), wid = Probanden ID, within = Innersubjekt-Faktor (Zeit)
rep_anova <- anova_test(data = reaktion_daten, dv = Reaktion, wid = ID, within = Zeit)Warning: NA detected in rows: 12,26,38,39.
Removing this rows before the analysis.# Blick auf die Sphärizität (Mauchly-Test)
rep_anova$`Mauchly's Test for Sphericity` Effect W p p<.05
1 Zeit 0.623 0.15 Check Sphärizität: Der p-Wert des Mauchly-Tests ist 0.15 (> 0.05). Die Voraussetzung wird nicht verworfen. Wir können die “normale” ANOVA ohne Korrektur interpretieren.
Wäre die Sphärizität verletzt (p<0.05), müssten wir die Freiheitsgrade korrigieren, um den Fehler 1. Art nicht zu erhöhen. Die Standardmethode hierfür ist die Greenhouse-Geisser-Korrektur. Das rstatix-Paket ist sehr komfortabel: Es wendet diese Korrektur automatisch an, wenn nötig. Du musst dich darum also nicht explizit kümmern, aber zur Auswertung sollte die Korrektur erwähnt werden.
# Ausgabe der ANOVA-Tabelle
get_anova_table(rep_anova, correction = "auto") # Bei Verletzung: "GG", sonst: "none"ANOVA Table (type III tests)
Effect DFn DFd F p p<.05 ges
1 Zeit 2 18 195.028 6.32e-13 * 0.941Interpretation des Outputs:
F: Der F-Wert ist mit 195 extrem hoch.
p: Der p-Wert liegt weit unter 0.001.
ges: Das generalisierte Eta-Quadrat liegt bei 0.94. Das ist ein extrem starker Effekt [Link], was bedeutet, dass fast die gesamte Variation in den Daten durch die Tageszeit erklärt werden kann.
Ergebnisformulierung: Die Analyse zeigt somit schonmal einen hochsignifikanten Unterschied in der Reaktionszeit zwischen den Tageszeiten (F(2,18) = 195, p < .001).
Aber wer unterscheidet sich von wem? Die ANOVA sagt uns nicht, zwischen welchen der drei Zeitpunkte Unterschiede existieren. Ist die Reaktion mittags schneller als abends? Oder morgens verschieden zu abends? Dafür brauchen wir einen Post-Hoc-Test [Link]. Der Standard ist der paarweise t-Test mit Messwiederholung und Bonferroni-Korrektur [Link].
# Paarweise t-Tests mit Bonferroni-Korrektur
pairwise.t.test(reaktion_daten$Reaktion, reaktion_daten$Zeit,
paired = TRUE,
p.adjust.method = "bonferroni")
Pairwise comparisons using paired t tests
data: reaktion_daten$Reaktion and reaktion_daten$Zeit
morgens mittags
mittags 1.6e-08 -
abends 3.0e-06 4.6e-07
P value adjustment method: bonferroni Auswertung: Alle Vergleiche (morgens-mittags, morgens-abends, mittags-abends) haben p-Werte < 0.05. Jede Tageszeit unterscheidet sich signifikant von jeder anderen.
Ergebnisformulierung: Die Varianzanalyse mit Messwiederholung zeigt einen signifikanten Haupteffekt der Tageszeit (F(2, 18) = 195, p < .001). Paarweise t-Tests (Bonferroni-korrigiert) bestätigen signifikante Unterschiede zwischen allen drei Zeitpunkten (alle p < .001).
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Hätten wir vorab eine Verletzung der Normalverteilung festgestellt (z.B. stark gekrümmte QQ-Plots), wäre der Friedman-Test [Link] die geeignete nicht-parametrische Alternative. Er arbeitet mit Rängen.
friedman.test(reaktion_daten$Reaktion, reaktion_daten$Zeit, reaktion_daten$ID)
Friedman rank sum test
data: reaktion_daten$Reaktion, reaktion_daten$Zeit and reaktion_daten$ID
Friedman chi-squared = 20, df = 2, p-value = 4.54e-05Auswertung: Auch der nicht-parametrische Test bestätigt das Ergebnis signifikanter Unterschiede. Der p-Wert liegt weit unter 0.001 und damit unter dem Signifikanzniveau von 0.05.
Post-Hoc für Friedman: Hier nutzt man meist den paarweisen Wilcoxon-Test mit p-Wert-Korrektur [Link] (z.B. nach Bonferroni), um Fehler durch multiples Testen [Link] zu vermeiden.
# Paarweiser Wilcoxon-Test
pairwise.wilcox.test(reaktion_daten$Reaktion, reaktion_daten$Zeit,
paired = T,
p.adjust.method = "bonferroni")Warning in wilcox.test.default(xi, xj, paired = paired, ...): kann bei
Bindungen keinen exakten p-Wert Berechnen
Warning in wilcox.test.default(xi, xj, paired = paired, ...): kann bei
Bindungen keinen exakten p-Wert Berechnen
Warning in wilcox.test.default(xi, xj, paired = paired, ...): kann bei
Bindungen keinen exakten p-Wert Berechnen
Pairwise comparisons using Wilcoxon signed rank test with continuity correction
data: reaktion_daten$Reaktion and reaktion_daten$Zeit
morgens mittags
mittags 0.0105 -
abends 0.0063 0.0127
P value adjustment method: bonferroni Auswertung: Auch hier sehen wir in der p-Wert-Tabelle, dass alle 3 Tests die Nullhypothese verwerfen und entsprechend ein Unterschied in der Reaktionsgeschwindigkeit zwischen allen drei Gruppen vorliegt.
Ergebnisformulierung: Der Friedman-Test bestätigt, dass sich die Reaktionsgeschwindigkeiten je nach Tageszeit signifikant unterscheiden (chisq = 20, p < .001). Der paarweise Wilcoxon-Test mit Bonferroni-Korektur als Post-Hoc-Überprüfung offenbart, dass sich auch alle drei Gruppen untereinander signifikant voneinander unterscheiden.
Szenario: Wir testen eine neue Lern-App. 30 Studierende machen einen Vokabeltest zu drei Zeitpunkten: - Start (Basislinie) - Woche2 (Nach 2 Wochen Training) - Woche4 (Nach 4 Wochen Training)
Forschungsfrage: Unterscheidet sich der Lernerfolg (gemessen an Scores) signifikant zwischen den Messzeiten?
Der Datensatz heißt lern_daten.
# Struktur prüfen:
str(lern_daten)'data.frame': 90 obs. of 3 variables:
$ ID : Factor w/ 30 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Zeit : Factor w/ 3 levels "Start","Woche2",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Score: num 45.4 45.1 41.4 47.1 40.4 ...# wir haben wieder eine metrische Zielvariable 'Score' und zwei kategoriale Gruppenvariablen 'ID' und 'Zeit'
# Welche Gruppen gibt es genau?
levels(lern_daten$Zeit)[1] "Start" "Woche2" "Woche4"# Die Gruppen wurden im Szenario beschrieben
# Zur Vollständigkeit des Prozedere prüfen wir auf vollständigkeit der Daten.
table(is.na(lern_daten$Score))
FALSE
90 # Wie hier erwartbar war gibt es keine Fehlwerte.
# statistischer Überblick --> 'summary'-Output.
summary(lern_daten) ID Zeit Score
1 : 3 Start :30 Min. :18.77
2 : 3 Woche2:30 1st Qu.:44.62
3 : 3 Woche4:30 Median :53.62
4 : 3 Mean :53.74
5 : 3 3rd Qu.:60.79
6 : 3 Max. :94.88
(Other):72 # Die Infos über die Anzahl von 30 Personen mit je 3 Zeitmesswerten war schon der Szenarienbeschreibung zu entnehmen
# Die Score-Verteilung schauen wir uns im nächsten Schritt genauer an# Boxplot: Lernscore (y) abhängig vom Messzeitpunkt (x)
boxplot(Score ~ Zeit, data = lern_daten,
main = "Lernscore nach Messzeitpunkt",
ylab = "Lernscore",
col = c("lightblue", "orange", "lightgreen"))
Visuelle Interpretation: Die Scores steigen im Zeitverlauf im allgemeinen an, wobei die Varianz der Scores nach 4 Wochen deutlich größer ist als in den Zeitgruppen davor. Aber sind diese Unterschiede auch signifikant?
1. Normalverteilung:
# grafisch --> QQ-Plot für jeden Zeitpunkt einzeln
library(ggpubr)
ggqqplot(lern_daten, "Score", facet.by = "Zeit")
Visuelle Interpretation: Die Punkte folgen abgesehen von der Gruppe Woche4 eher weniger der Referenzgeraden. Besonders an den Randbereichen erkennen wir bei Start und Woche2 systematische Abweichungen. Die Punkte liegen allerdings noch größtenteils innerhalb des grauen Konfidenzbereichs. Zur besseren Einschätzung blicken wir auf die Normalverteilungstests:
# statistisch 1 --> Shapiro-Wilk-Test auf NV
shapiro.test(subset(lern_daten$Score, lern_daten$Zeit == "Start"))
Shapiro-Wilk normality test
data: subset(lern_daten$Score, lern_daten$Zeit == "Start")
W = 0.91157, p-value = 0.0163shapiro.test(subset(lern_daten$Score, lern_daten$Zeit == "Woche2"))
Shapiro-Wilk normality test
data: subset(lern_daten$Score, lern_daten$Zeit == "Woche2")
W = 0.92161, p-value = 0.02954shapiro.test(subset(lern_daten$Score, lern_daten$Zeit == "Woche4"))
Shapiro-Wilk normality test
data: subset(lern_daten$Score, lern_daten$Zeit == "Woche4")
W = 0.96594, p-value = 0.4347Hiermit bekommen wir bestätigt, dass nicht alle Gruppen einer Normalverteilung folgen, da für Start und Woche2 die p-Werte < 0.05 sind. Als Konsequenz der Voraussetzungs-Verletzung werden wir die Hypothese nicht mit der normalen ANOVA testen, sondern mit dem nichtparametrischen Friedman-Test.
2. Sphärizität: Obwohl schon feststeht, dass wir nicht mit der normalen ANOVA rechnen schauen wir uns zur Vollständigkeit das Ergebnis des Mauchly-Tests auf Sphärizität an, der Teil des ANOVA-Outputs im Paket `rstatix`` ist.
# statistisch 2 --> Sphärizität (Mauchly-Test)
library(rstatix)
rep_anova <- anova_test(data = lern_daten, dv = Score, wid = ID, within = Zeit)
rep_anova$`Mauchly's Test for Sphericity` Effect W p p<.05
1 Zeit 0.435 8.66e-06 *Check Sphärizität: Der p-Wert des Mauchly-Tests ist < 0.05. Die Voraussetzung wird also ebenfalls verworfen.
Wieso? Bei den Lernkurven verhält es sich hier so, dass die Messungen, zu Beginn auf einem ähnlichen Score liegen und der Lernerfolg beim Großteil der Personen zu Beginn gleichmäßig anwächst. Dieser Lerneffekt lässt bei manchen Personen nach mehr als 2 Wochen offensichtlich nach, wohingegen andere sich weiter verbessern. Demnach steigt die Spanne zwischen den Lernstarken und den Stagnierenden und genau dieser Variabilitätswandel verletzt die Sphärizität.
Wäre die NV-Annahme vorher nicht verletzt würde get_anova_table(rep_anova, correction = "GG") zum Ergebnis führen. Im nächsten Schritt folgt aber die in diesem Sachverhalt korrekte Methode des Friedman-Tests.
# Analyse durchführen --> Friedman-Test
friedman.test(lern_daten$Score, lern_daten$Zeit, lern_daten$ID)
Friedman rank sum test
data: lern_daten$Score, lern_daten$Zeit and lern_daten$ID
Friedman chi-squared = 56.267, df = 2, p-value = 6.051e-13Auswertung: Der Test bestätigt das Ergebnis signifikanter Unterschiede. Der p-Wert liegt weit unter 0.001 und damit unter dem Signifikanzniveau von 0.05.
# Post-Hoc-Analyse durchführen --> Paarweiser Wilcoxon-Test
pairwise.wilcox.test(lern_daten$Score, lern_daten$Zeit,
paired = T,
p.adjust.method = "bonferroni")
Pairwise comparisons using Wilcoxon signed rank exact test
data: lern_daten$Score and lern_daten$Zeit
Start Woche2
Woche2 5.6e-09 -
Woche4 5.6e-09 1.8e-07
P value adjustment method: bonferroni Auswertung: Auch hier sehen wir in der p-Wert-Tabelle, dass alle 3 Tests die Nullhypothese klar verwerfen und entsprechend ein Unterschied in den Scores zwischen allen drei Gruppen vorliegt.
Ergebnisformulierung: Der Friedman-Test zeigt einen signifikanten Unterschied zwischen den Scores in der Lernapp zu verschiedenen Messzeitpunkten (chisq = 41.04, p < .001). Der paarweise Wilcoxon-Test mit Bonferroni-Korektur als Post-Hoc-Überprüfung offenbart, dass sich auch alle drei Gruppen untereinander signifikant voneinander unterscheiden (alle p < .001).