Hintergrund: Wir betrachten den Datensatz PlantGrowth. Er enthält das Gewicht (weight) von Pflanzen unter drei verschiedenen Bedingungen (group): 1. ctrl: Kontrollgruppe 2. trt1: Behandlung 1 3. trt2: Behandlung 2
Mehr Infos sehen wir über:
Frage: Unterscheidet sich das Pflanzengewicht signifikant zwischen den drei Gruppen? Hinweis: Da wir hier mehr als zwei Gruppen vergleichen, dürfen wir nicht einfach viele t-Tests rechnen (Fehlerkumulierung! [Link]). Wir nutzen die Varianzanalyse (ANOVA).
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("PlantGrowth")
# Struktur prüfen: Welche Datentypen liegen vor?
str(PlantGrowth)'data.frame': 30 obs. of 2 variables:
$ weight: num 4.17 5.58 5.18 6.11 4.5 4.61 5.17 4.53 5.33 5.14 ...
$ group : Factor w/ 3 levels "ctrl","trt1",..: 1 1 1 1 1 1 1 1 1 1 ...# Welche Gruppen gibt es genau?
levels(PlantGrowth$group)[1] "ctrl" "trt1" "trt2"Struktur: Wir haben eine metrische Zielvariable (weight) und eine kategoriale Gruppenvariable (group) mit 3 Stufen (> 2). Das ist der klassische Fall für eine einfaktorielle ANOVA.
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 weight?
table(is.na(PlantGrowth$weight))
FALSE
30 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 (>5) sein sollte. Für genauere Instruktionen dieser Bereinigung folge dem Link [Link].
Zusätzlich überprüfen wir mit der summary- Funktion, in welchem Wertebereich die Gewichtswerte liegen und wie viele Beobachtungen pro Treatment-Gruppe enthalten sind.
summary(PlantGrowth) weight group
Min. :3.590 ctrl:10
1st Qu.:4.550 trt1:10
Median :5.155 trt2:10
Mean :5.073
3rd Qu.:5.530
Max. :6.310 Wir sehen dass wir pro Gruppe 10 Beobachtungen für die Analyse vorliegen haben.
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 von Gruppen.
# Boxplot: Gewichtswerte (y) abhängig von der Behandlungsgruppe (x)
boxplot(weight ~ group, data = PlantGrowth,
main = "Pflanzengewicht nach Gruppe",
xlab = "Behandlungsgruppe", ylab = "Trockengewicht",
col = c("grey", "lightblue", "orange"))
Visuelle Interpretation: trt1 scheint im Schnitt leichtere Pflanzen zu haben als die Kontrolle. trt2 scheint schwerere Pflanzen zu haben. Die Frage ist: Sind diese Unterschiede groß genug, um statistisch signifikant zu sein?
Bezug Glossar: Verteilungs-Checks
Für die einfaktorielle 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 kurz erstellen (aov ist der Befehl für ANOVA)
model_temp <- aov(weight ~ group, data = PlantGrowth)
# Histogramm der Residuen
hist(resid(model_temp), main = "Histogramm der Residuen", col = "lightgrey")
Die Verteilung sieht relativ symmetrisch aus, was dafür spricht, dass die Voraussetzung hier gilt. Um das konkret zu prüfen führen wir einen Shapiro-Wilk-Test [Link] durch, da dieser uns eine statistisch kräftigere Aussage zur Normalverteilungsannahme gibt, als rein visuelle Vergleiche.
# Test auf Normalverteilung
# H0: Die Stichprobe stammt aus einer Normalverteilung.
# Shapiro-Test auf die Residuen
shapiro.test(resid(model_temp))
Shapiro-Wilk normality test
data: resid(model_temp)
W = 0.96607, p-value = 0.4379Ergebnisauswertung: Der p-Werte beträgt 0.438 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 einfaktorielle Varianzanalyse erfüllt.
Neben der Normalverteilungsannahme wird für 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(weight ~ group, data = PlantGrowth)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 2 1.1192 0.3412
27 Der Levene-Test liefert einen p-Wert von 0.341, was deutlich über dem Signifikanzniveau von 0.05 liegt. Das bedeutet die Varianzen der Gruppen sind nicht signifikant unterschiedlich und wir dürfen die klassische ANOVA rechnen
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Wir berechnen die ANOVA ohne Messwiederholung [Link], da die einzelnen Treatment-Gruppen statistisch nicht voneinander abhängen [Link: Unabhängig vs. Verbunden].
Die Nullhypothese der ANOVA lautet: Alle Gruppenmittelwerte sind gleich.
# ANOVA berechnen
anova_model <- aov(weight ~ group, data = PlantGrowth)
# Ergebnis anzeigen
summary(anova_model) 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
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nun schauen wir abschließend auf den Output. Dabei ist in erster Linie der p-Wert [Link] in der Spalte Pr(>F) relevant zur Beantwortung unserer Anfangshypothese. Dieser Wert liegt unter 0.05.
Hinweis: Der ANOVA-Output selbst zeigt uns nicht die Mittelwerte der Gruppen (die Spalte Mean Sq steht für ‘Mean Squares’, also die erklärten Varianzen, nicht für das durchschnittliche Gewicht). Um zu sehen, welche Gruppe die höchsten zentralen Werte hat, könnten wir nochmal auf unseren Boxplots aus 2. zurückschauen.
Ergebnisformulierung: Die Analyse zeigt einen signifikanten Unterschied in den Pflanzengewichten zwischen den Treatment-Gruppen (F(2,27) = 4.85, p = .016).
ABER: Die ANOVA sagt uns nicht, wer sich von wem unterscheidet. Ist Trt1 anders als Ctrl? Oder Trt1 anders als Trt2? Dafür brauchen wir einen Post-Hoc-Test [Link]. Der Standard ist der “Tukey HSD”.
# Post-Hoc Test (Tukey Honest Significant Differences)
TukeyHSD(anova_model) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = weight ~ group, data = PlantGrowth)
$group
diff lwr upr p adj
trt1-ctrl -0.371 -1.0622161 0.3202161 0.3908711
trt2-ctrl 0.494 -0.1972161 1.1852161 0.1979960
trt2-trt1 0.865 0.1737839 1.5562161 0.0120064Detail-Auswertung: Wir schauen auf die p-Werte (p adj) der Paarvergleiche:
trt1-ctrl: p = 0.391 (Kein signifikanter Unterschied)
trt2-ctrl: p = 0.198 (Kein signifikanter Unterschied)
trt2-trt1: p = 0.012 (Signifikant!)
Ergebnisformulierung: Die Varianzanalyse zeigt einen signifikanten Haupteffekt der Behandlung (F(2, 27) = 4.85, p = .016). Der Post-Hoc-Test nach Tukey offenbart, dass sich jedoch nur die Gruppen “Treatment 1” und “Treatment 2” signifikant voneinander unterscheiden (p = .012). Die Behandlungen unterscheiden sich nicht signifikant von der Kontrollgruppe.
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Hätten wir in Kapitel 3 eine Verletzung der Voraussetzungen festgestellt, wäre der Kruskal-Wallis-Test die geeignete nicht-parametrische Alternative. Dieser ist robust gegen Verletzungen der Normalverteilung und basiert auf Rängen.
# Durchführung des Kruskal-Wallis-Tests
test_ergebnis_npar <- kruskal.test(weight ~ group, data = PlantGrowth)
# Ergebnis ausgeben
test_ergebnis_npar
Kruskal-Wallis rank sum test
data: weight by group
Kruskal-Wallis chi-squared = 7.9882, df = 2, p-value = 0.01842Auswertung: Auch der nicht-parametrische Test bestätigt das Ergebnis signifikanter Unterschiede. Der p-Wert liegt mit 0.018 wieder unter dem Signifikanzniveau von 0.05.
Post-Hoc für Kruskal-Wallis: 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(PlantGrowth$weight, PlantGrowth$group,
p.adjust.method = "bonferroni")Warning in wilcox.test.default(xi, xj, paired = paired, ...): kann bei
Bindungen keinen exakten p-Wert Berechnen
Pairwise comparisons using Wilcoxon rank sum test with continuity correction
data: PlantGrowth$weight and PlantGrowth$group
ctrl trt1
trt1 0.596 -
trt2 0.189 0.027
P value adjustment method: bonferroni Auswertung: Auch hier sehen wir in der p-Wert-Tabelle, dass der kleinste p-Wert (0.027) zwischen trt1 und trt2 liegt und die anderen beiden Tests die Hypothese der Gleichheit in den zentralen Tendenzen nicht verwerfen.
Ergebnisformulierung: Der Kruskal-Wallis-Test zeigt einen signifikanten Haupteffekt der Behandlung (chisq = 7.99, p = .018). Der paarweise Wilcoxon-Test als Post-Hoc-Überprüfung offenbart, dass sich jedoch nur die Gruppen “Treatment 1” und “Treatment 2” signifikant voneinander unterscheiden (p = .027). Die Behandlungen unterscheiden sich nicht signifikant von der Kontrollgruppe.
Szenario: Der Datensatz chickwts vergleicht das Gewicht von Küken bei 6 verschiedenen Futterdiäten.
Forschungsfrage: Gibt es signifikante Unterschiede bei den Kükengewichten bei vreschiedenen Futtermitteln? Und falls ja, Unterscheiden sich sunflower (Sonnenblumen) signifikant von linseed (Leinsamen)?
# Datensatz laden (in R vorinstalliert)
data("chickwts")
# Struktur prüfen:
str(chickwts)'data.frame': 71 obs. of 2 variables:
$ weight: num 179 160 136 227 217 168 108 124 143 140 ...
$ feed : Factor w/ 6 levels "casein","horsebean",..: 2 2 2 2 2 2 2 2 2 2 ...# wir haben wieder eine metrische Zielvariable 'weight' und eine kategoriale Gruppenvariable 'feed'
# Welche Gruppen gibt es genau?
levels(chickwts$feed)[1] "casein" "horsebean" "linseed" "meatmeal" "soybean" "sunflower"# Jetzt wissen wir welche 6 Futtermittel miteinander verglichen werden.
# Zur Vollständigkeit des Prozedere prüfen wir auf vollständigkeit der Daten.
table(is.na(chickwts$weight))
FALSE
71 # Wie hier erwartbar war gibt es keine Fehlwerte.
# Wie viele Beobachtungen liegen pro Treatment-Gruppe vor?
summary(chickwts) weight feed
Min. :108.0 casein :12
1st Qu.:204.5 horsebean:10
Median :258.0 linseed :12
Mean :261.3 meatmeal :11
3rd Qu.:323.5 soybean :14
Max. :423.0 sunflower:12 # Die Anzahl der Hühner pro Futtermittel sind nicht identisch und im Bereich zwischen 10 und 14.
# Das ist aber vollkommen unproblematisch, da die Analyse hier keine gleiche Stichprobengröße voraussetzt.# Visualisierung --> Boxplot
boxplot(weight ~ feed, data = chickwts,
main = "Kükengewicht nach Futterart", col = "gold")
Visuelle Interpretation: sunflower und casein scheinen die höchsten Gewichte zu erzielen. horsebean wirkt am schlechtesten. Es ist offensichtlich, dass es Unterschiede in den Gewichten gibt. Aber sind diese Unterschiede auch signifikant?
# zuerst ANOVA Modell erstellen
model_chick <- aov(weight ~ feed, data = chickwts)
# grafisch --> Histogramm der Residuen
hist(resid(model_chick), main = "Histogramm der Residuen", col = "lightgrey")
Die Verteilung sieht relativ symmetrisch aus, was dafür spricht, dass die Voraussetzung hier gilt.
# statistisch 1 --> Test auf Normalverteilung der Residuen
shapiro.test(resid(model_chick))
Shapiro-Wilk normality test
data: resid(model_chick)
W = 0.98616, p-value = 0.6272# p > 0.05 -> Normalverteilung nicht abgelehnt.
# statistisch 2 --> Varianzhomogenität
library("car") # Sicherstellen, dass das Paket geladen ist
leveneTest(weight ~ feed, data = chickwts)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 5 0.7493 0.5896
65 # p > 0.05 -> Varianzhomogenität nicht abgelehnt.Ergebnisauswertung: Die Testergebnisse sprechen nicht gegen die Voraussetzungen der Normalverteilung der Residuen und der Varainzhomogenität. Daher kann im folgenden die einfaktorielle ANOVA berechnet werden mit anschließenden post-hoc-Tests.
# Analyse durchführen --> ANOVA
# wurde bei 3. bereits berechnet und kann hier einfach ausgewertet werden mit 'summary'
summary(model_chick) Df Sum Sq Mean Sq F value Pr(>F)
feed 5 231129 46226 15.37 5.94e-10 ***
Residuals 65 195556 3009
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ergebnisformulierung: Die Analyse zeigt einen deutlich signifikanten Unterschied in den Hühnergewichten zwischen den Futterarten (F(5,65) = 15.37, p < .001).
# Post-Hoc-Analyse durchführen --> Tukey-Test
TukeyHSD(model_chick) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = weight ~ feed, data = chickwts)
$feed
diff lwr upr p adj
horsebean-casein -163.383333 -232.346876 -94.41979 0.0000000
linseed-casein -104.833333 -170.587491 -39.07918 0.0002100
meatmeal-casein -46.674242 -113.906207 20.55772 0.3324584
soybean-casein -77.154762 -140.517054 -13.79247 0.0083653
sunflower-casein 5.333333 -60.420825 71.08749 0.9998902
linseed-horsebean 58.550000 -10.413543 127.51354 0.1413329
meatmeal-horsebean 116.709091 46.335105 187.08308 0.0001062
soybean-horsebean 86.228571 19.541684 152.91546 0.0042167
sunflower-horsebean 168.716667 99.753124 237.68021 0.0000000
meatmeal-linseed 58.159091 -9.072873 125.39106 0.1276965
soybean-linseed 27.678571 -35.683721 91.04086 0.7932853
sunflower-linseed 110.166667 44.412509 175.92082 0.0000884
soybean-meatmeal -30.480519 -95.375109 34.41407 0.7391356
sunflower-meatmeal 52.007576 -15.224388 119.23954 0.2206962
sunflower-soybean 82.488095 19.125803 145.85039 0.0038845Wir können jetzt einzeln sehen, welche Gruppen sich signifikant unterscheiden. Zur Beantwortung der anfänglichen Forschungsfrage interessier uns die Zeile sunflower-linseed. Hier liegt der p-Wert unter dem Signifikanzniveau von 0.05, womit wir einen signifikanten Unterschied im mittleren Gewicht bestätigen können.
Ergebnisformulierung: Die Varianzanalyse zeigt einen signifikanten Haupteffekt der Futtermittel-Wahl (F(5, 65) = 15.37, p < .001). Der Post-Hoc-Test nach Tukey offenbart, dass Sonnenblumenfutter zu signifikant höherem Gewicht führt als Leinsamenfutter (p < .001).