load("MODALIS.RData")1 Beispieldaten
Auch für diese Einheit ziehen wir den Datensatz MODALIS heran. Die Variable pubertaet findet sich aber ausschließlich in der Vollversion mit allen Variablen.
2 ANOVA bei Gruppenvergleichen (Gruppierungen mit mehr als 2 Levels)
Sind wir an Gruppenvergleichen von zwei oder mehr Gruppen interessiert, wird üblicherweise eine ANOVA als Testinstrument herangezogen. Für die Beispieldaten aus MODALIS ergibt sich eine solche Datenkonstellation, wenn wir nach Sprungkraftunterschieden in Abhängigkeit der Pubertätsphase fragen. Die Daten stellen sich wie folgt dar (wir reduzieren den Datensatz zunächst auf die Gruppe der Jugendlichen):
J <- subset(MODALIS.full, gruppe=="J" & !is.na(pubertaet))
# Faktorvariante der Pubertätsvariable erstellen. Original behalten für Plot-Jitter!
J$pubertaet.fac <- as.factor(J$pubertaet)
# Halbtransparente Farbdefinition
farbe <- rgb(red=0.1, green= 0.3, blue=0.8, alpha=0.3)
# Gruppenmittel errechnen
MW <- tapply(J$sprungkraft, J$pubertaet, FUN=mean, na.rm=TRUE)
# Boxplot mit Punkten überlagert plus Mittelwerte in rot
boxplot(sprungkraft ~ pubertaet, data=J, xlab="Pubertät", ylab="Sprungkraft [cm]")
points(sprungkraft ~ jitter(pubertaet, 0.3), data=J, pch=21, bg=farbe, col=farbe)
points(x = c(1,2,3), y=MW, col="red", bg="red", pch = 22)
2.1 Berechnung der ANOVA (F-Test)
Es gibt zahlreiche Wege in R, eine ANOVA zu berechnen. Für unseren Fall voneinander unabhängiger Gruppen existiert ein einfacher Weg mit der Standardfunktion aov(). Sobald abhängige Daten ins Spiel kommen nutzen wir eine andere Funktion – damit befassen wir uns auch in einem anderen Tutorial.
AOV <- aov(sprungkraft ~ pubertaet.fac, data=J)
summary(AOV) Df Sum Sq Mean Sq F value Pr(>F)
pubertaet.fac 2 2533 1266.7 27.76 8.29e-12 ***
Residuals 307 14009 45.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
13 Beobachtungen als fehlend gelöschtDas hochsignifikante Ergebnis muss durch Post-hoc Tests spezifiziert werden, da der F-Test alleine keinen Aufschluss über die Ausprägung der Unterschiede zwischen den Gruppen ergibt.
PW.t <- pairwise.t.test(J$sprungkraft, J$pubertaet.fac, p.adjust.method = "bonf")
PW.t
Pairwise comparisons using t tests with pooled SD
data: J$sprungkraft and J$pubertaet.fac
1 2
2 0.013 -
3 7.0e-12 1.6e-07
P value adjustment method: bonferroni 3 ALM bei Gruppenvergleichen
3.1 Alternative zur ANOVA per ALM
Das ALM erstellen wir mit nahezu derselben Syntax wie die ANOVA. Die im Hintergrund ablaufende Dummy-Codierung müssen wir nicht beachten, R erledigt dies autormatisch, sobald wir eine kategoriale (Faktor-)Variable als Prädiktor verwenden.
LM <- lm(sprungkraft ~ pubertaet.fac, data=J)
summary(LM)
Call:
lm(formula = sprungkraft ~ pubertaet.fac, data = J)
Residuals:
Min 1Q Median 3Q Max
-18.3344 -4.8887 -0.5244 4.5827 19.3585
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 25.6916 0.7552 34.018 < 2e-16 ***
pubertaet.fac2 2.6728 0.9259 2.887 0.00417 **
pubertaet.fac3 8.0499 1.1014 7.309 2.35e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6.755 on 307 degrees of freedom
(13 Beobachtungen als fehlend gelöscht)
Multiple R-squared: 0.1531, Adjusted R-squared: 0.1476
F-statistic: 27.76 on 2 and 307 DF, p-value: 8.287e-12Die Funktion aov() kann aus einem linearen Modell auch eine ANOVA-Ausgabe erzeugen:
summary(aov(LM)) Df Sum Sq Mean Sq F value Pr(>F)
pubertaet.fac 2 2533 1266.7 27.76 8.29e-12 ***
Residuals 307 14009 45.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
13 Beobachtungen als fehlend gelöschtAuch daran erkennt man, dass ANOVA und lineare Modellierung per lm() äquivalent und inenander überführbar sind. Der einzige Unterschied in den Modellausgaben besteht darin, dass wir bei der ANOVA und den post-hoc Tests automatisch alle Gruppenvergleiche berechnen können, beim ALM hingegen nur jene Koeffizienten enthalten sind, die die Vergleiche zur Referenzgruppe repräsentieren.
3.2 Gezielte Gruppenvergleiche im ALM
Achtung
Eine den post-hoc Tests vergleichbare Korrektur der p-Werte (für Mehrfachvergleiche) findet im ALM zunächst nicht statt. Dies ist insbesondere bei kategorialen Prädiktoren (d. h. Gruppenvergleichen) mit zahlreichen Levels problematisch, insbesondere auch wenn alle Gruppenvergleiche von Interesse sind. Die Korrektur sollte daher händisch nachgeholt werden.
3.2.1 Weg 1: traditionell, aber kompliziert:
Die auf Mehrfachtestung korrigierten p-Werte erzeugt man aus den bereits unkorrigiert vorliegenden p-Werten mit der Funktion p.adjust() wie folgt:
SUMMARY.LM <- summary(LM)
# Modell 1: Pubertätsstufe 1 als Referenzgruppe
p.adjust(SUMMARY.LM$coefficients[,4], method="bonf") (Intercept) pubertaet.fac2 pubertaet.fac3
1.098402e-105 1.251329e-02 7.040831e-12 Auf diese Weise lassen sich die post-hoc Tests aus der ANOVA inkl. Bonferroni-Korrektur auch aus dem ALM erzeugen. Allerdings ist auch in diesem Fall der Vergleich zwischen den Gruppen 2 und 3 nicht vorhanden, da das ALM hierfür keinen Koeffizienten liefert. Man müsste daher zunächst umständlich das ALM umstellen (Faktorlevels umsortieren), um diesen Gruppenvergleich 2 vs. 3 gezielt zu berechnen und statistisch zu testen.
3.2.2 Weg 2: eleganter, aber Tukey only
Für p.adjust() steht eine Vielzahl an Korrekturmethoden zur Verfügung, siehe dazu die Hilfefunktion zu p.adjust(). Allerdings beinhalten diese Optionen nicht die ebenfalls populäre Methode der Tukey Honest Significant Difference (Korrektur nach Tukey). Stattdessen existiert die Funktion TukeyHSD(), welche sogar einfacher zu bedienen ist und direkt die Ausgabe eines ALM akzeptiert, wenn man dieses per aov() direkt als ANOVA Ausgabe formatiert. Damit entfallen lästige händische Korrekturvorgänge bzw. das Umstellen der Faktorlevels komplett:
TukeyHSD(aov(LM)) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = LM)
$pubertaet.fac
diff lwr upr p adj
2-1 2.672778 0.4920607 4.853494 0.0115891
3-1 8.049924 5.4559947 10.643854 0.0000000
3-2 5.377147 3.1063407 7.647953 0.0000002Die Ausgabe enthält neben der Mittelwertsdifferenz (diff) sowohl die obere und untere Konfidenzgrenze der Differenz (lower: lwr, upper: upr) sowie die korrigierten p-Werte.
Fazit
Auch bei Hypothesentests zu Unterschieden zwischen mehreren unabhängigen Gruppen kann man ein ALM einsetzen. Die paarweisen Gruppenvergleiche lassen sich auch daraus erzeugen, entweder mittels pairwise.t.test() oder TukeyHSD().
4 Multiple lineare Regression (mehrere metrische Prädiktoren)
Natürlich können wir in ähnlicher Form die Sprungkraft nicht nur aus einer mehrstufigen kategorialen Variable (pubertaet) vorhersagen, sondern dafür auch mehrere metrische Prädiktoren einsetzen. Damit erzeugen wir ebenso ein lineares Modell mit mehr als einer Prädiktorvariable (analog zu den Dummies für die Faktorstufen), jedoch verläuft dann die Interpretation der Koeffizienten (wie auch im Fall der einfachen linearen Regression) etwas anders.
In MODALIS finden wir Beispieldaten in Form der Vorhersage der Sprungkraft aus Daten der Aktionsschnelligkeit (Anzahl Fußtaps in 10 Sekunden) und der Handkraft (gemessen in kp). Wir prüfen zunächst, ob die beiden Prädiktorvariablen unkorreliert sind (d. h. unabhängig voneinander sind):
plot(J$aktionsschn, J$handkraft, xlab="Aktionsschnelligkeit [Taps]", ylab ="Handkraft [kp]")
text(x=25, y=50, label=paste0("r=",round(cor(J$aktionsschn,J$handkraft,
use="complete.obs"),3)))
Da die Prädiktoren nicht bedeutsam korreliert sind, können wir sie berechtigterweise als unabhängige Prädiktoren in das lineare Modell einbinden:
LM2 <- lm(sprungkraft ~ aktionsschn + handkraft, data=J)
summary(LM2)
Call:
lm(formula = sprungkraft ~ aktionsschn + handkraft, data = J)
Residuals:
Min 1Q Median 3Q Max
-24.0601 -4.1841 -0.4025 4.0587 15.4901
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 16.14149 2.57057 6.279 1.19e-09 ***
aktionsschn -0.04773 0.07806 -0.611 0.541
handkraft 0.54000 0.03928 13.746 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.719 on 301 degrees of freedom
(19 Beobachtungen als fehlend gelöscht)
Multiple R-squared: 0.3972, Adjusted R-squared: 0.3932
F-statistic: 99.17 on 2 and 301 DF, p-value: < 2.2e-16Im Ergebnis erkennen wir, dass die Aktionsschnelligkeit keinen signifikanten Beitrag zum Modell zu leisten scheint, der p-Wert ist zu groß, das Modell könnte womöglich ohne diesen Prädiktor auskommen. Im Modell ist sie durch den Koeffizient aktionsschn repräsentiert. Für jede Steigerung der Aktionsschnelligkeit um 1 Tap verändert sich die Vorhersage der Sprungkraft um ca. -0.048 cm. Diese Veränderung ist so klein, dass sie nicht sicher von Null abweicht und daher nicht signifikant wird.
Mit der Handkraft verhält es sich vergleichbar, wenngleich auch hochsignifikant: mit jeder Steigerung der Handkraft um 1 kp verändert sich die Sprunghöhe um ca. 0.54 cm.
Wie auch im Fall der einfachen linearen Regression interessiert die Ausprägung des Intercept hier nicht. Eine Vorhersage Sprunghöhe für Personen mit Aktionsschnelligkeit und Handkraft von Null ist schlichtweg irrelevant und kommt in der Praxis niemals vor.
5 Grafische Modellprüfung
Wie immer lassen sich bei allen Varianten des ALM die Voraussetzungen und Annahmen per plot() Funktion prüfen. Hier im Beispiel des ursprünglichen Modells (ANOVA mit 3 Gruppen):
plot(LM, which=c(1,2))
DIe Voraussetzungen sind erüllt, die Residualverteilung unkritisch. Daher steht einer Anwendung und Interpretation eines solchen Modells nichts im Weg.
Auch das Modell mit metrischen Prädiktoren kann auf diese Weise geprüft werden:
plot(LM2, which=c(1:2))
Auch hier erscheint die Verteilung der Residuale keinem Muster zu folgen oder anderweitige Auffälligkeiten zu beinhalten, so dass wir auch die Modellvoraussetzungen in diesem Fall als erfüllt ansehen können.
Fazit
Natürlich unterscheiden sich die Residualplots der beiden Modelle erheblich, da metrische Prädiktoren zahlreiche (stetige) Vorhersagewerte generieren, kategoriale Prädiktoren hingegen nur jene Anzahl an Abstufungen vorhersagen, die als Faktorlevels zur Verfügung stehen.
Ende | Stand: 2024-04-18