Professur Forschungsmethoden und Analyseverfahren in der Biomechanik
Datum
11. April 2024
1 Datensatz und Fragestellung
Wir nutzen den Datensatz MODALIS, der einige nicht normalverteilte Variablen enthält. Insbesondere die Reaktionsschnelligkeit stellt eine solche Variable dar, die zumeist linkssteile Verteilungen aufweist. Wir bilden einen Subdatensatz der Senioren (Alter über 60) und unterteilen diese Gruppe nochmals nach “jungen” (bis.65) und “alten” Senioren (ueber.65).
Hier wird eine neue Variable erzeugt, die den Wert bis.65 annimmt, wenn das Alter <= 65 Jahre ist, andernfalls den Wert ueber.65 (logische Prüfung nach dem if - else Schema, siehe Funktion ifelse()).
Anhand dieser Daten könnte man nun der Frage nach dem Gruppenunterschied zwischen jüngeren und älteren Senioren nachgehen. Die Nullhypothese \(H_0\) wäre in diesem Fall, dass sich beide Altersgruppen hinsichtlich der Reaktionsschnelligkeit nicht unterscheiden. In einem solchen Fall müsste man von gleichen Gruppenmittelwerten ausgehen - ein t-Test wäre eine Möglichkeit, dies anhand der Stichprobendaten zu prüfen.
Jedoch sind nicht alle Voraussetzungen für einen t-Test vollumfänglich erfüllt. Die Verteilungsschiefe wirft Zweifel auf, inwiefern ein Mittelwert repräsentativ für die jeweiligen Gruppen wäre. Wir prüfen die Daten auf Normalverteilung in beiden Gruppen:
Shapiro-Wilk normality test
data: subset(SEN$reaktionsschn, SEN$group == "ueber.65")
W = 0.97113, p-value = 0.02826
Die bereits visuell ersichtliche Verteilungsschiefe bestätigt sich im Shapiro-Wilk Test als signifikante Testergebnisse in beiden Gruppen. Wir sollten daher davon ausgehen, dass die Daten nicht normalverteilt sind, und in Konsequenz auch von einem klassischen t-Test absehen und stattdessen ein alternatives Verfahren zur Hypothesenprüfung eines Gruppenunterschieds wählen.
1.1 Exkurs: Verteilungsschiefe in MODALIS$reaktionsschn ist kein Artefakt kleiner Stichproben!
Achtung: Dass es sich hierbei nicht um ein Phänomen kleiner Stichproben und “zufällig” erzeugter Verteilungsschiefen handelt, kann man daran erkennen, dass diese Konstellation bei Nutzung des gesamten Datensatzes sogar noch deutlicher zum Vorschein kommt. Die folgenden Abbildung zeigt dieselbe Variable über alle Probanden inkl. Lage des Mittelwerts und des Medians als vertikale Linien.
Einzeichnen der vertikalen Linie für den Mittelwert
2
Einzeichnen der vertikalen Linie für den Median
Durch die Verteilungsschiefe liegen Mittelwert und Median sichtbar auseinander, der Mittelwert liegt näher am Maximum und weniger “im Zentrum” als der Median. Mittelwerte derartiger Variablen zwischen Gruppen zu vergleichen erscheint daher zumindest fragwürdig. Wir ersetzen ihn durch einen Vergleich der Mediane, welche robuster hinsichtlich der Verteilungsschiefe sind. Hierfür stünden bekannte Testverfahren zur Verfügung (z. B. der Mann-Whitney-U-Test). Aus Gründen der Illustration des Bootstraps umgehen wir diesen Test hier und erzeugen uns ein 95%-Konfidenzintervall für die Differenz der Mediane aus beiden Altersgruppen.
Nach \(10^{4}\) Iterationen stellt sich die Verteilung der Differenzen der Mediane wie folgt dar:
P2 <-ggplot(THETA, aes(x=Med.Diff))+geom_histogram(fill="dodgerblue", color="white", bins=11)+theme_bw()+labs(x=bquote(theta~": Differenz der Mediane"), y="Häufigkeit")plot(P2)
Wir erkennen eine Verschiebung der Verteilung in den negativen Bereich, was aufgrund der Berechnungsweise der Differenz und der Bedeutung der gemessenen Größe darauf hindeutet, dass die Gruppe der älteren Senioren höhere Werte und damit zumindest tendenziell etwas schlechtere Reaktionszeitwerte erreicht. Das Zentrum dieser Verteilung liegt im Mittel bei -1.11 und im Median bei -1.12 cm Unterschied zwischen den Gruppen.
2.1 Ersatz für Signifikanztest: Wahrscheinlichkeiten anhand der Häufigkeitsverteilung ermitteln
Das Äquivalent zu einem Signifikanztest können wir dadurch erzeugen, indem wir beurteilen, wie wahrscheinlich die Nulldifferenz (Äquivalent zur \(H_0\) im konventionellen Signifikanztest) in unserer erhaltenen Verteilung von THETA (\(\theta\)) ist. Liegt die Nulldifferenz im Bereich der mittleren 95% der Verteilung, so besteht kein Grund davon auszugehen, dass sie unwahrscheinlich sei. Läge die Nulldifferenz außerhalb der mittleren 95% unserer Verteilung, wäre sie hingegen hinreichend unwahrscheinlich, um die Schlussfolgerung für einen Gruppenunterschied zu begründen. Ein Blick auf die Verteilung zeigt schon, dass letzteres womöglich nicht der Fall sein wird.
Wollen wir “einseitig testen” und hätten eine bergündete einseitige Hypothese, die eine Verschiebung der Verteilung beispielsweise in den negativen Wertebereich nach sich zöge, wäre ein einseitiges \(\alpha\) von 5% ausreichend. D. h. wir würden in einem solchen Fall – da wir aufgrund der Berechnungsweise tendenziell negative Werte der Gruppendifferenz erwarten – das 95. Perzentil (am oberen, rechten Rand) als kritische Grenze für die Lage der Nulldifferenz suchen (da die gesamte Verteilung ja tendenzell nach links von der Null abweicht).
Wollten wir aufgrund einer ungerichteten Hypothese zweiseitig testen, müssten wir auf beiden Seiten der Verteilung den kritischen Wert bei jeweils 2.5% der verbleibenden Fälle ziehen. Dies übersetzt sich in das 2.5 und 97.5-Perzentil der erhaltenen Verteilung von THETA.
# Einseitig bei 5%, Annahme bis.65 < ueber.65Q.EINSEIT <-quantile(THETA$Med.Diff, probs =0.95)# Zweiseitig bei insgesamt 5% (2.5% auf jeder Seite)Q.ZWEISEIT <-quantile(THETA$Med.Diff, probs =c(0.025, 0.975))
Ein Vergleich der Obergrenze (einseitig: Q.EINSEIT= 0.75 oder beiden zweiseitigen Werten Q.ZWEISEIT= -3.5, 1.25) mit dem Wert 0 (Nulldifferenz, kein Unterschied der Mediane) zeigt, dass wir in keinem der Fälle eine auf dem 5%-Niveau signifikante Differenz der Mediane vorfinden, da die Null noch unterhalb der einseitigen Obergrenze, bzw. innerhalb der zweiseitigen Ober- und Untergrenze liegt. Auf Basis der Stichprobendaten können wir also keine Evidenz für einen Unterschied der beiden Altersgruppen finden.
3 Für Fortgeschrittene: Bootstrap mit dem Paket boot
Achtung
Die Nutzung des Pakets boot setzt voraus, dass man sich in den Umgang mit selbst geschriebenen Funktionen einarbeitet. Wer dazu keine Energie aufwenden will, bleibt besser bei der händischen Lösung mittels Schleife.
Die Installation des Pakets vorausgesetzt (install.packages("boot")) können wir die obige Prozedur auch in alternativer Weise durchführen, und anstatt eine Schleifenkonstruktion zu schreiben eine entsprechende Funktion in den Code integrieren. Dadurch wird der Code nicht weniger, aber bei komplexeren Vorhaben etwas übersichtlicher.
3.1 Einführendes Beispiel: Mittelwertschätzung per Bootstrap
Um das Paket zunächst in seiner Funktionsweise zu illustrieren, beginnen wir mit der Schätzung eines Mittelwerts durch die Bootstrap-Prozedur (was angesichts der Verteilungsschiefe der Daten mit Vorsicht zu genießen ist). In einem zweiten Schritt illustrieren wir etwas kompliziertere Vorhaben wie die obige Berechnung einer Differenz der Mediane zweier Gruppen.
Das Kernstück der Prozedur bildet die Funktion MWcalc, in welcher die eigentlichen Bestandteile der Bootstrap-Iterationen untergebracht werden. Innerhalb dieser Funktion muss also die Schätzung von \(\theta\) erfolgen – im Beispiel ist dies der Mittelwert des Bootstrap Samples. Wir definieren eine solche Funktion wie folgt:
Die Funktionsargumente data und indices werden im Kontext der Funktion boot (aus dem gleichnamigen Paket boot) erzeugt, in welche wir MWcalc() einbetten:
data enthät die zu bootstrappenden Daten
indices enthält die – wie der Name sagt – Indices der Datenreihen, welche durch Zufallsziehung (mit Zurücklegen) die neue Stichprobenzusammensetzung des Bootstrap Samples definiert.
Diese Funktion MWcalc() wird von boot() iterativ aufgerufen und erzeugt in jedem Durchlauf ein Boostrap-Sample, sowie die daraus berechneten Größe MW. Dieses Iterieren und “Sammeln” der jeweiligen Ausprägungen von MW erledigen wir mit der Funktion boot() wie folgt:
Die Liste S.boot enthält neben zahlreichen(!) weiteren Ausgaben das eigentlich interessante Ergebnis des Bootstrap: eine Sammlung der iterativ erzeugten Mittelwerte der Bootstrap samples. Diese findet sich in S.boot$t. Damit könnte man sich (ähnlich wie oben) ein Konfidenzintervall händisch errechnen, aber zumindest eine Verteilung darstellen:
hist(S.boot$t)
Man muss das KI nun nicht händisch errechnen, das Paket boot enthält dazu sehr fortgeschrittene Funktionen. Dazu nutzt man die Funktion boot.ci(), welche mit dem Argument type=... unterschiedliche Arten von KIs erzeugen kann. Diese so erzeugbaren KIs sind nahezu allesamt mit anderen (aufwändigeren) Methoden berechnet als wir dies oben “von Hand” erledigt haben. Am ehesten vergleichbar zu unserer obigen Methode ist die “Perzentilmethode” type="perc". Interessierte finden die für die exakte Berechnung relevante Literatur in der Hilfe von R zur Funktion boot.ci(). Aus praktischer Perspektive sind die Unterschiede zwischen den Arten der KI meist weniger gravierend. Unterschiedliche Arten der KI erzeugen in der Regel keine anderen Ergebnisse bzw. Interpretationen, und resultieren meist lediglich in geringfügig anders ausgeprägten Konfidenzgrenzen.
# Ein an Perzentilen orientiertes 95%-KI:KI <-boot.ci(S.boot, type="perc")
Die Hilfe zu boot.ci() verrät, dass sich die eigentlichen Konfidenzgrenzen in der Ausgabe KI$percent an 4. und 5. Stelle des Vektors befinden. Zusammen mit dem Mittelwert (er findet sich unter KI$t0) lässt sich das KI zusammenfassen:
c(UG = KI$percent[4], MW = KI$t0, OG = KI$percent[5])
UG MW OG
20.51397 21.42877 22.37709
3.2 Anwendung von boot() auf obiges Beispiel: Differenz der Mediane
Um das obige Beispiel der Differenz der Mediane unter Einsatz des Pakets boot zu replizieren, müssen wir die selbst definierte Funktion MWcalc() durch eine Funktion ersetzen, welche die Differenz der Mediane berechnet. Da wir mit zwei Gruppen agieren, müssen wir diese in der Funktion entsprechend abbilden und im Kontext von boot() als sogenannte strata definieren.
KI <-boot.ci(S2.boot, type="perc")c(KI$perc[4], KI$t0, KI$perc[5])
[1] -3.500 -1.000 1.375
Damit erhalten wir dieselben Konfidenzgrenzen wie oben “von Hand” (Q.ZWEISEIT), und natürlich ein ebenfalls “nicht signifikantes” Ergebnis im Gruppenvergleich der Mediane.
