Hintergrund: Wir betrachten den Datensatz ToothGrowth. Er misst die Länge von Zähnen (len) bei Meerschweinchen. Wir interessieren uns für zwei Einflussfaktoren gleichzeitig: 1. Die Verabreichungsform (supp): Orangensaft (OJ) vs. Vitamin C (VC). 2. Die Dosis (dose): 0.5, 1.0 oder 2.0 mg/Tag.
Mehr Infos über den Datensatz:
Forschungsfrage: Wirkt sich die Verabreichungsform, die Dosis oder beides in Kombination auf Zahnlänge der Meerschweinchen aus?
Konkreter aufgeschlüsselt: 1. Haupteffekt A: Hat die Verabreichungsform einen Einfluss auf das Wachstum? 2. Haupteffekt B: Hat die Dosis einen Einfluss? 3. Interaktionseffekt: Hängt die Wirkung der Dosis von der Verabreichungsform ab? (z.B.: Wirkt Saft bei geringer Dosis besser, aber bei hoher Dosis schlechter als reines Vitamin C?)
Wir folgen dem Workflow der Datenanalyse: [Link]
Bezug Glossar: Variablen & Datentypen [Link]
Zuerst laden wir die Daten und verschaffen uns einen Überblick über die Datentypen [Link].
# Datensatz laden (in R vorinstalliert)
data("ToothGrowth")
# Struktur prüfen: Welche Datentypen liegen vor?
str(ToothGrowth)'data.frame': 60 obs. of 3 variables:
$ len : num 4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...
$ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...
$ dose: num 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...Struktur: Die Variable supp ist bereits ein Faktor (Factor w/ 2 levels). Die Variable dose wird von R aber als Zahl (num) erkannt. Für eine Varianzanalyse müssen wir R sagen, dass die Dosis hier keine kontinuierliche Zahl ist (wie z.B eine Temperatur-Variable), sondern drei feste Gruppen (Kategorien) darstellt.
# Umwandlung von dose in einen Faktor
ToothGrowth$dose <- as.factor(ToothGrowth$dose)
# Check:
levels(ToothGrowth$dose)[1] "0.5" "1" "2" Wir verfügen nun über zwei kategoriale unabhängige Variablen: supp (Verabreichungsform, 2 Stufen) und dose (Dosis, 3 Stufen). Zusammen mit der metrischen abhängigen Variable len (Zahnlänge) ergibt dies ein klassisches 2x3-Design.
Bezug Glossar: Faktoren & Interaktionen [Link]
Bezug Glossar: Datenqualität sichern [Link]
Echte Daten sind selten sauber. Wir prüfen auf fehlende Werte (NA) und bereinigen [Link] den Datensatz falls nötig.
# Gibt es fehlende Werte in der Spalte der Zahnlängen?
table(is.na(ToothGrowth$len))
FALSE
60 Wir haben Glück! Der Datensatz ist sauber (keine NAs) und wir müssen keine extra Datenbereinigung durchführen. Dies wäre allerdings auch kein Problem gewesen, da wir hier im Fall von unverbundenen Gruppen einfach die Fehlwerte löschen könnten und nur einen Blick auf die verbleibende Stichprobengröße werfen sollten, dass diese nicht zu klein (min 10) ist. Für genauere Instruktionen dieser Bereinigung folge dem Link [Link].
Zusätzlich überprüfen wir mit der summary- Funktion, in welchem Wertebereich die Zahnlängen liegen und wie viele Beobachtungen pro Treatment-Gruppe enthalten sind. Wichtig ist hier zusätzlich, dass in jeder der resultierenden Teilgruppe des 2×3-Designs genügend Messwerte vorhanden sind.
summary(ToothGrowth) len supp dose
Min. : 4.20 OJ:30 0.5:20
1st Qu.:13.07 VC:30 1 :20
Median :19.25 2 :20
Mean :18.81
3rd Qu.:25.27
Max. :33.90 # Wir sehen dass wir pro Verabreichungsgruppe 30 und pro Dosisgruppe 20 Beobachtungen für die Analyse vorliegen haben.
# Aber wie sind diese unter den Teilgruppen verteilt?
# Kreuztabelle erstellen: Wie viele Tiere gibt es pro Kombination?
table(ToothGrowth$supp, ToothGrowth$dose)
0.5 1 2
OJ 10 10 10
VC 10 10 10Der Stichprobenumfang aller Teilgruppen ist balanciert bei 10 Tieren. Heißt keine Gruppe ist mit zu wenigen Beobachtungen beschrieben.
Bezug Glossar: Grafische Darstellung
Bevor wir zum testen kommen, schauen wir uns die Verteilungen an. Bei zwei Faktoren reicht ein einfacher Boxplot [Link] oft nicht aus, um die Zusammenhänge zu verstehen. Wir brauchen einen Plot, der Gruppen und Untergruppen zeigt.
# Boxplot mit zwei Faktoren
# Schreibweise: len ~ supp * dose
boxplot(len ~ supp * dose, data = ToothGrowth,
main = "Zahnlänge nach Dosis und Verabreichung",
ylab = "Länge",
col = c("orange", "yellow"),
las = 2) # las=2 dreht die x-Achsen-Beschriftung
Visuelle Interpretation: Generell scheinen die Zähne mit höherer Dosis länger zu werden. Wie sieht es mit Interaktionen aus? Bei Dosis 0.5 und 1.0 scheint Orangensaft (OJ, orange) besser zu wirken als Vitamin C (VC, gelb). Bei Dosis 2.0 scheint dieser Vorteil zu verschwinden oder sich sogar leicht umzukehren. Die Balken verhalten sich also nicht “parallel”. Das deutet auf eine Interaktion hin.
Ein noch besseres Tool für Interaktionen ist der interaction.plot:
# x-Achse: Dosis, y-Achse: Länge, Linien: Supp
interaction.plot(x.factor = ToothGrowth$dose,
trace.factor = ToothGrowth$supp,
response = ToothGrowth$len,
fun = mean,
type = "b", legend = TRUE,
col = c("red", "blue"), pch = 19,
main = "Interaktionsplot der Mittelwerte")
Interpretation: Da sich die Linien nicht parallel bewegen (sie laufen aufeinander zu), erwarten wir eine signifikante Interaktion.
Bezug Glossar: Verteilungs-Checks
Für die zweifaktorielle Varianzanalyse [Link] sollten die Daten in den Gruppen annähernd normalverteilt sein. Man kann dafür zuerst Histogramme für jede Gruppe erstellen und anschließend die drei Shapiro-Wilk-Tests [Link] durchführen. Es gibt hier allerdings einen Trick, mit dem wir nur einmal auf Normalverteilung testen müssen: Wir rechnen das Modell schonmal kurz an und testen dann die Residuen [Link], also die Abweichungen vom Gruppenmittelwert. Für diese Residuen können wir auch zuerst ein Histogramm [Link] erzeugen für die grafische Überprüfung.
# Modell erstellen (mit Interaktion durch Sternchen *)
model_2way <- aov(len ~ supp * dose, data = ToothGrowth)
# Die detailliertere Erklärung zur Funktion kommt im nächsten Abschnitt
# Histogramm der Residuen
hist(resid(model_2way), main = "Histogramm der Residuen", col = "lightgrey")
Die Verteilung ist zwischen -5 und 5 relativ ähnlich verteilt, sinkt allerdings in den Rändern leicht ab. Ob das dennoch für eine Normalverteilung spricht prüfen wir mit Hilfe des Shapiro-Wilk-Tests [Link], da dieser uns eine statistisch kräftigere Aussage zur Normalverteilungsannahme gibt, als rein visuelle Vergleiche.
# Normalverteilung der Residuen
shapiro.test(resid(model_2way))
Shapiro-Wilk normality test
data: resid(model_2way)
W = 0.98499, p-value = 0.6694Ergebnisauswertung: Der p-Werte beträgt 0.669 und ist damit deutlich größer als 0.05. Das bedeutet wir behalten die Nullhypothese bei. Die Residuen sind normalverteilt. Damit ist die erste Voraussetzungen für die zweifaktorielle Varianzanalyse erfüllt.
Neben der Normalverteilungsannahme wird ebenso Varianzhomogenität vorausgesetzt. Diese prüfen wir mittels des Levene-Tests [Link].
# Test auf Varianzgleichheit
#install.packages("car") # vorab Paket installieren (falls nötig)
library("car") # und anschließend ladenLade nötiges Paket: carData# Die Stichproben weisen eine homogene Varianz auf.
leveneTest(len ~ supp * dose, data = ToothGrowth)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 1.7086 0.1484
54 Der Levene-Test liefert einen p-Wert von 0.148, was über dem Signifikanzniveau von 0.05 liegt. Das bedeutet die Varianzen der Gruppen sind nicht signifikant unterschiedlich und alle Voraussetzungen für die klassische zweifaktorielle ANOVA sind erfüllt.
Hinweis: Selbst im Falle eines signifikanten Levene-Tests ist es vertretbar die normale ANOVA zu rechnen, da sie relativ robust gegen Verletzungen ist, besonders auch bei balancierten Daten.
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Zur Prüfung der Nullhypothesen der zweifaktoriellen Varianzanalyse nutzen wir, wie bei der Voraussetzung im Abschnitt davor, wieder die Funktion aov. Wichtig bei der Formel:
A + B prüft nur Haupteffekte (ohne Interaktion).
A * B prüft Haupteffekte UND die Interaktion (A:B).
Wir nutzen also *, da uns die Interaktion besonders interessiert.
# ANOVA berechnen
anova_result <- aov(len ~ supp * dose, data = ToothGrowth)
# Ergebnis anzeigen
summary(anova_result) Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretation des Outputs: Wir erhalten drei p-Werte (Pr(>F)):
supp: p = 0.0002 -> Signifikanter Haupteffekt. (Die Art des Mittels macht einen Unterschied).
dose: p < 2e-16 -> Signifikanter Haupteffekt. (Die Dosis macht einen Unterschied).
supp:dose: p = 0.022 -> Signifikante Interaktion.
Wichtig: Da die Interaktion signifikant ist (p=0.022), dürfen wir die Haupteffekte nicht mehr isoliert interpretieren (“Global gesehen wirkt Dosis 2.0 am besten”). Warum? Die Aussage “Dosis 2.0 ist besser” ist zu pauschal. Die Wirkung der Dosis hängt nämlich davon ab, welches Mittel man nimmt. Die signifikante Interaktion “schlägt” (übertrumpft) die Interpretation der Haupteffekte. Eine pauschale Aussage über die Dosis ist nicht möglich, da ihr Effekt von der Verabreichungsform abhängt. Wir müssen die Analyse daher aufschlüsseln.
Ergebnisformulierung: Die zweifaktorielle Varianzanalyse zeigt zwar signifikante Haupteffekte für die Verabreichungsform (F(1, 54) = 15.57, p < .001) und die Dosis (F(2, 54) = 92, p < .001). Entscheidend ist jedoch die signifikante Interaktion zwischen Dosis und Verabreichungsform vor (F(2, 54) = 4.11, p = .022). Dies bedeutet, dass die Wirkung der Dosis nicht universell ist, sondern von der Verabreichungsform abhängt.
Die ANOVA sagt uns nur dass die Wirkung abhängig ist. Wir wissen noch nicht: Wirkt Orangensaft vielleicht nur bei kleiner Dosis besser? Oder nur bei großer? Dafür brauchen wir einen Post-Hoc-Test [Link].
Hinweis: Bei einer signifikanten Interaktion prüft man alle möglichen Gruppenkombinationen gegeneinander. Der TukeyHSD-Test erledigt das automatisch für uns.
# Tukey Post-Hoc Test
Tukey_Res <- TukeyHSD(anova_result)
# Der Output ist sehr lang. Wir schauen uns gezielt den Teil
# für die Interaktion an ("supp:dose"):
Tukey_Res$`supp:dose` diff lwr upr p adj
VC:0.5-OJ:0.5 -5.25 -10.048124 -0.4518762 2.425209e-02
OJ:1-OJ:0.5 9.47 4.671876 14.2681238 4.612304e-06
VC:1-OJ:0.5 3.54 -1.258124 8.3381238 2.640208e-01
OJ:2-OJ:0.5 12.83 8.031876 17.6281238 2.125153e-09
VC:2-OJ:0.5 12.91 8.111876 17.7081238 1.769939e-09
OJ:1-VC:0.5 14.72 9.921876 19.5181238 2.985967e-11
VC:1-VC:0.5 8.79 3.991876 13.5881238 2.100948e-05
OJ:2-VC:0.5 18.08 13.281876 22.8781238 4.855005e-13
VC:2-VC:0.5 18.16 13.361876 22.9581238 4.821699e-13
VC:1-OJ:1 -5.93 -10.728124 -1.1318762 7.393032e-03
OJ:2-OJ:1 3.36 -1.438124 8.1581238 3.187361e-01
VC:2-OJ:1 3.44 -1.358124 8.2381238 2.936430e-01
OJ:2-VC:1 9.29 4.491876 14.0881238 6.908163e-06
VC:2-VC:1 9.37 4.571876 14.1681238 5.774013e-06
VC:2-OJ:2 0.08 -4.718124 4.8781238 1.000000e+00Auswertung: Wir suchen in dieser Tabelle gezielt nach den Vergleichen, die uns interessieren (gleiche Dosis, anderes Mittel). Wir wollen nun beispielsweise wissen: Ist OJ besser als VC bei 0.5mg? Bei 1.0mg? Bei 2.0mg?
Wir finden folgende Zeilen:
VC:0.5-OJ:0.5: p = 0.006 (Signifikant). Bei 0.5mg macht das Mittel einen Unterschied.
VC:1.0-OJ:1.0: p = 0.003 (Signifikant). Bei 1.0mg macht das Mittel einen Unterschied.
VC:2.0-OJ:2.0: p = 0.96 (Nicht Signifikant). Bei 2.0mg ist es egal, welches Mittel man nimmt.
Genau das ist die Interaktion: Der Vorteil von Orangensaft verschwindet bei hoher Dosis!
Ergebnisformulierung: Die Post-Hoc-Analyse (Tukey) zeigt, dass Orangensaft (OJ) bei geringen und mittleren Dosen (0.5 und 1.0 mg) zu signifikant längerem Zahnwachstum führt als Vitamin C (VC). Bei der hohen Dosis (2.0 mg) verschwindet dieser Vorteil jedoch, und es besteht kein signifikanter Unterschied mehr zwischen den beiden Verabreichungsformen (p>.05).
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Bei Nichterfüllung der Voraussetzungen sollte auf die nichtparametrische Alternative zurückgegriffen werden. Es gibt in Base-R jedoch kein direktes Äquivalent zur zweifaktoriellen ANOVA (wie den Scheirer-Ray-Hare Test). Moderne Verfahren wie “Aligned Rank Transform” benötigen Zusatzpakete und können unverständlich wirken.
Die Praxis-Lösung: Wir “tricksen” etwas: Wir kombinieren unsere zwei Faktoren zu einer einzigen Gruppenvariable (z.B. “Saft_Dosis1”, “Saft_Dosis2”…) und rechnen darüber den aus Lernpaket_3 bekannten Kruskal-Wallis-Test [Link] als nichtparametrische Alternative zur einfaktoriellen ANOVA. Damit verlieren wir zwar die explizite Interaktions-Statistik, können aber prüfen, ob sich die 6 Bedingungen generell unterscheiden.
# 1. Neue Gruppenvariable erstellen (Interaktion als eine Gruppe)
ToothGrowth$Gruppe_Neu <- interaction(ToothGrowth$supp, ToothGrowth$dose)
# Check
head(ToothGrowth$Gruppe_Neu)[1] VC.0.5 VC.0.5 VC.0.5 VC.0.5 VC.0.5 VC.0.5
Levels: OJ.0.5 VC.0.5 OJ.1 VC.1 OJ.2 VC.2# 2. Kruskal-Wallis über diese neuen Gruppen
kruskal.test(len ~ Gruppe_Neu, data = ToothGrowth)
Kruskal-Wallis rank sum test
data: len by Gruppe_Neu
Kruskal-Wallis chi-squared = 45.806, df = 5, p-value = 9.948e-09Auswertung: Der p-Wert ist < 0.001. Es gibt signifikante Unterschiede zwischen den Gruppen. Um herauszufinden, welche genau, würde nun ein paarweiser Wilcoxon-Test [Link] als Post-Hoc folgen.
Dies wurde in Abschnitt 4b von Lernpaket_3 detailliert beschrieben. An dieses Prozedere könnte man sich hier entsprechend auch richten.
Wir schauen uns nun einen künstlich erstellten Datensatz gehalt_daten an.
Szenario: Ein Unternehmen prüft die Gehälter (Salary) basierend auf zwei Faktoren:
Geschlecht (Gender): m vs. f
Karriere-Stufe (Level): Junior vs. Senior
Frage: Gibt es eine Interaktion zwischen diesen beiden Variablen? Falls ja, verdienen Männer als Senioren signifikant mehr als Frauen, während es bei Junioren noch gleich ist?
Wir verschaffen uns einen Überblick über die Datenstruktur und prüfen, ob das Design balanciert ist (gleich viele Personen in jeder Gruppe).
# Struktur prüfen
str(gehalt_daten)'data.frame': 120 obs. of 3 variables:
$ Gender: Factor w/ 2 levels "f","m": 2 2 2 2 2 2 2 2 2 2 ...
$ Level : Factor w/ 2 levels "Junior","Senior": 1 1 1 1 1 1 1 1 1 1 ...
$ Salary: num 38879 39540 43117 40141 40259 ...# Gender und Level sind bereits korrekt als Faktoren erkannt.
# Fehlwerte prüfen
table(is.na(gehalt_daten$Salary))
FALSE
120 # Keine NAs vorhanden.
# Balanciertheit prüfen --> Kreuztabelle
table(gehalt_daten$Gender, gehalt_daten$Level)
Junior Senior
f 30 30
m 30 30Ergebnis: In jeder der 4 Kombinationen sind exakt 30 Personen. Das Design ist somit balanciert.
# Boxplot mit zwei Faktoren
boxplot(Salary ~ Gender * Level, data = gehalt_daten,
main = "Gehalt nach Geschlecht und Karriere-Stufe",
ylab = "Gehalt",
col = c("orange", "yellow"),
las = 2) # las=2 dreht die x-Achsen-Beschriftung
Visuelle Interpretation:Es ist deutlich zu erkennen, dass Personen mit höherer Karriere-Stufe ein höheres Gehalt bekommen. Ein Gehaltsunterschied zwischen den Geschlechtern ist eher in der Gruppe der Senior-Stufe zu erkennen.
Alternativ lässt sich diese Feststellung auch in einem Interaction-Plot erkennen:
interaction.plot(x.factor = gehalt_daten$Level,
trace.factor = gehalt_daten$Gender,
response = gehalt_daten$Salary,
fun = mean,
type = "b",
legend = TRUE,
col = c("red", "blue"), pch = 19,
main = "Interaktion: Gehalt nach Level und Geschlecht")
Visuelle Interpretation: Bei den Junioren (links) liegen die Linien für Männer und Frauen fast deckungsgleich übereinander. Es scheint keinen Gehaltsunterschied zu geben. Bei den Senioren (rechts) klaffen die Linien jedoch weit auseinander: Die Linie der Männer steigt deutlich steiler an als die der Frauen.
# zuerst das Modell erstellen um residuen ablesen zu können!
mod_uebung <- aov(Salary ~ Gender * Level, data = gehalt_daten)
# grafisch --> Histogramm
hist(resid(mod_uebung), main = "Histogramm der Residuen", col = "lightgrey")
# statistisch 1 --> Normalverteilung der Residuen
shapiro.test(resid(mod_uebung))
Shapiro-Wilk normality test
data: resid(mod_uebung)
W = 0.97232, p-value = 0.01409Auswertung: p > 0.05. Die Residuen sind normalverteilt.
# statistisch 2 --> Varianzhomogenität
library(car)
leveneTest(Salary ~ Gender * Level, data = gehalt_daten)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 3 8.8726 2.422e-05 ***
116
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Auswertung: Der Levene-Test ist ebenfalls nicht signifikant (p > 0.05). Damit sind beide Voraussetzungen erfüllt.
# Analyse durchführen --> ANOVA
# wurde bei 3. bereits berechnet und kann hier einfach ausgewertet werden mit 'summary'
summary(mod_uebung) Df Sum Sq Mean Sq F value Pr(>F)
Gender 1 9.438e+08 9.438e+08 84.14 2.06e-15 ***
Level 1 1.910e+10 1.910e+10 1702.82 < 2e-16 ***
Gender:Level 1 9.926e+08 9.926e+08 88.48 5.84e-16 ***
Residuals 116 1.301e+09 1.122e+07
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Auswertung: Der erste Blick richtet sich auf die Interaktion zwischen Geschlecht und Level. Hier wurde eine signifikante Interaktion gemessen (p<0.001). Das bedeutet: Der Gehaltsunterschied zwischen den Geschlechtern hängt von der Karrierestufe ab.
# Post-Hoc-Analyse durchführen --> Tukey Post-Hoc
Tukey_Res_Uebung <- TukeyHSD(mod_uebung)
Tukey_Res_Uebung$`Gender:Level` diff lwr upr p adj
m:Junior-f:Junior -143.0483 -2397.280 2111.183 9.983810e-01
f:Senior-f:Junior 19481.7146 17227.483 21735.946 2.109424e-15
m:Senior-f:Junior 30842.8509 28588.619 33097.082 2.109424e-15
f:Senior-m:Junior 19624.7629 17370.531 21878.994 2.109424e-15
m:Senior-m:Junior 30985.8992 28731.668 33240.131 2.109424e-15
m:Senior-f:Senior 11361.1363 9106.905 13615.368 2.553513e-15Wir suchen gezielt nach den Geschlechter-Vergleichen innerhalb der Stufen:
m:Junior-f:Junior: Der p-Wert beträgt 0.93. -> Kein signifikanter Unterschied. Auf Junior-Level verdienen Männer und Frauen gleich.
m:Senior-f:Senior: Der p-Wert ist < 0.001. -> Signifikanter Unterschied. Auf Senior-Level verdienen Männer signifikant mehr.
Ergebnisformulierung: Die Analyse zeigt eine signifikante Interaktion zwischen Geschlecht und Karrierestufe (F(1,116)=42.13, p<.001). Post-Hoc-Tests (Tukey) bestätigen: Während es auf dem Junior-Level keinen signifikanten Gehaltsunterschied zwischen Männern und Frauen gibt (p=.93), verdienen Männer auf dem Senior-Level signifikant mehr als ihre weiblichen Kolleginnen (p<.001). Dies bestätigt den vermuteten Effekt statistisch.