Lerneinheit 2: Gruppenvergleiche verbundene Stichproben (t-Test/Wilcoxon)

Szenario & Forschungsfrage

Hintergrund: Wir betrachten den Datensatz sleep. Er enthält Daten über die zusätzliche Schlafdauer (in Stunden) von 10 Patienten, die zwei verschiedene Schlafmittel verabreicht bekamen. Wichtig: Jeder Patient hat beide Medikamente getestet (Messwiederholung).

Mehr Infos über den Datensatz:

Frage: Wirkt Schlafmittel 2 stärker schlaffördernd als Schlafmittel 1?

Wir folgen dem Workflow der Datenanalyse: [Link]


1a. Daten einlesen & die “Paarung” verstehen

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("sleep")

# Struktur prüfen: Welche Datentypen liegen vor?
str(sleep)
'data.frame':   20 obs. of  3 variables:
 $ extra: num  0.7 -1.6 -0.2 -1.2 -0.1 3.4 3.7 0.8 0 2 ...
 $ group: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
 $ ID   : Factor w/ 10 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...

Das Besondere an diesem Datensatz: Wir sehen die Spalten ID und group sind Faktor-Variablen.

Die Zeile, in der ID = 1 und group = 1steht, gehört zur selben Person wie die Zeile ID = 1 und group = 2. Das ist das Kennzeichen für verbundene Daten [Link] (Paired Samples).

Die Spalte extra beschreibt eine numerische Variable und ist entsprechend intervall-skaliert.


1b. Datenaufbereitung (Data Cleaning) & Differenzberechnung

Bezug Glossar: Datenqualität sichern [Link]

Echte Daten sind nicht immer sauber. Wir prüfen auf fehlende Werte (NA) und bereinigen [Link] den Datensatz falls nötig.

# Gibt es fehlende Werte in der Spalte 'extra'? 
table(is.na(sleep$extra))

FALSE 
   20 

Wir haben glücklicherweise keine Fehlwerte und wir müssen keine extra Datenbereinigung durchführen. Da wir hier eine verbundene Stichprobe haben wäre diese Datenbereinigung auch etwas komplexer. Wir dürften nicht einfach einzelne Fehlwerte aus dem Datensatz entfernen, da wir damit den verbund von Wertepaaren zerstören würden. Das konkrete Prozedere zur Bereinigung ist über den Link [Link] nachzuverfolgen.

Um verbundene Tests sauber durchzuführen (und Voraussetzungen zu prüfen), ist es oft hilfreich, die Differenz zwischen den beiden Messzeitpunkten zu berechnen.

Logik: Differenz = Medikament 2 - Medikament 1.

# Wir extrahieren die Werte für Gruppe 1 und Gruppe 2
# (Dabei vertrauen wir darauf, dass der Datensatz nach ID sortiert ist)
werte_med1 <- sleep$extra[sleep$group == 1]
werte_med2 <- sleep$extra[sleep$group == 2]

# Wir berechnen die Differenz pro Patient
differenzen <- werte_med2 - werte_med1

# Kurzer Blick auf die Differenzen
summary(differenzen)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.00    1.05    1.30    1.58    1.70    4.60 

Interpretation: Ein positiver Wert bedeutet, dass der Patient mit Medikament 2 länger geschlafen hat.

Die Übersicht der Differenzen spricht schonmal sehr für die aufgestellte Forschungsfrage. Nichtsdestotrotz muss dieser Zusammenhang statistisch quantifiziert werden.


2. Explorative Datenanalyse (Visualisierung)

Bezug Glossar: Grafische Darstellung

Bevor wir zum testen kommen, schauen wir uns die Verteilungen an. Bei verbundenen Stichproben interessieren wir uns primär für die Verteilung der Differenzen, die wir hier mit einem Boxplot [Link] gut darstellen können.

# Boxplot der Differenzen (liegt der Median bei 0?)
boxplot(differenzen, 
        main = "Differenz der Schlafdauer (Medikament 2 - 1)", 
        ylab = "Stunden", col = "lightgreen")
abline(h = 0, col = "red", lty = 2) # Rote Linie bei 0 zur Orientierung

Visuelle Interpretation: Die Box liegt fast komplett über der roten 0-Linie. Das deutet ebenfalls darauf hin, dass Medikament 2 im Schnitt besser wirkt.


3. Voraussetzungen prüfen

Bezug Glossar: Verteilungs-Checks

Hier passiert ein häufiger Fehler: Wir testen nicht Gruppe 1 und Gruppe 2 separat auf Normalverteilung. Wir testen nur die Differenz auf Normalverteilung. Wir prüfen das zuerst grafisch mit einem Histogramm [Link] und anschließend mit einem Shapiro-Wilk-Test [Link], da dieser uns zusätzlich eine statistisch kräftige Aussage zur Normalverteilungsannahme der einzelnen Stichproben gibt, als der rein visuelle Vergleich.

# Histogramm der Differenzen
hist(differenzen, main = "Verteilung der Differenzen", 
     xlab = "Differenz (Std)", col = "lightgreen")

# Shapiro-Wilk-Test auf die Differenzen
shapiro.test(differenzen)

    Shapiro-Wilk normality test

data:  differenzen
W = 0.82987, p-value = 0.03334

Ergebnisauswertung: Das Histogramm zeigt schon starke Abweichungen zu einer Glockenkurve. Allerdings ist die Stichprobe mit nur 10 Elementen auch sehr gering, was eine klare Visualisierung erschwert.

Der Shapiro-Wilk-Test liefert einen p-Wert von 0.033, was knapp unter dem Signifikanzniveau von 0.05 liegt. Daher wird die Normalverteilungsannahme hier verworfen.

Streng genommen ist der t-Test hier also nicht ideal, weshalb wir später den Wilcoxon-Test betrachten. Da wir in diesem Lernmaterial allerdings die Durchführung des Tests einmal demonstrieren wollen, führen wir ihn dennoch einmal durch.

Hinweis: Bei Stichproben mit mehr als 30 Elementen sind Verletzungen generell eher unproblematisch. Das trifft aber auf unser Beispiel hier nicht zu.


4a Die Analyse - parametrisch - t-Test für verbundene Stichproben

Bezug Glossar: Parametrische Verfahren (Mittelwerte)

Wir berechnen einen t-Test für abhängige Stichproben [Link], da pro Person zwei zeitverzögerte Messwerte verglichen werden [Link: Unabhängig vs. Verbunden].

Dafür nutzen wir, wie auch im unverbundenen Fall die Funktion t.test, müssen aber zwingend das Argument paired = TRUE setzen.

# Durchführung des t-Tests für verbundene Stichproben
test_ergebnis <- t.test(werte_med2,werte_med1, paired = T)

# Ergebnis ausgeben
test_ergebnis

    Paired t-test

data:  werte_med2 and werte_med1
t = 4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
 0.7001142 2.4598858
sample estimates:
mean difference 
           1.58 

Nun schauen wir abschließend auf den Output. Dabei ist in erster Linie der p-Wert [Link] relevant zur Beantwortung unserer Anfangshypothese. Dieser Wert liegt weit unter 0.05.

Zusätzlich werden am Ende des Outputs noch der Mittelwerte der Differenzen zwischen Medikation 2 - Mediaktion 1 angegeben. Wir können hier erkennen, dass dieser Wert weit genug von 0 entfernt liegt, als weit entfernt von Gleicheit der Gruppen.

Ergebnisformulierung: Die Analyse zeigt einen signifikanten Unterschied in der zusätzlichen Schlafdauer zwischen den beiden Medikationsgruppen (t = 4.06, p = .003).

Um die statistische Aussage zu präzisieren und konkret unsere Forschungsfrage zu beantworten, führen wir den Test einseitig (gerichtet) durch. Wir prüfen, ob Medikation 2 signifikant längere Schlafzeiten erzielt. Da im Funktionsaufruf werte_med2 als erstes genannt wird, prüfen wir, ob “med2 > med1” ist (alternative = "greater")

# Durchführung des einseitigen t-Tests für unverbundene Stichproben
test_ergebnis_oneway <- t.test(werte_med2,werte_med1, paired = T,
                               alternative = "greater")

# Ergebnis ausgeben
test_ergebnis_oneway

    Paired t-test

data:  werte_med2 and werte_med1
t = 4.0621, df = 9, p-value = 0.001416
alternative hypothesis: true mean difference is greater than 0
95 percent confidence interval:
 0.8669947       Inf
sample estimates:
mean difference 
           1.58 

Ergebnisformulierung: Die Analyse bestätigt, dass die zusätzliche Schlafdauer bei Medikament 2 signifikant höher ist als bei Medikament 1 (t = 4.06, p = .001).


4b Die Analyse - nicht-parametrisch - Wilcoxon-Vorzeichen-Rang-Test

Bezug Glossar: Parametrische Verfahren (Mittelwerte)

Wie wir unter Kapitel 3 sehen konnten, waren die Voraussetzungen für den t-Test verletzt. Deshalb greifen wir zur sichereren, nicht-parametrischen Alternative. Im Falle von zwei verbundenen Stichproben ist diese der Wilcoxon(-Vorzeichen-Rang)-Test [Link] (in R als wilcox.test implementiert).

# Durchführung des Wilcoxon-Vorzeichen-Rang-Test
test_ergebnis_npar <- wilcox.test(werte_med2, werte_med1,
                                  paired = T, alternative = "greater", exact = T)                            
Warning in wilcox.test.default(werte_med2, werte_med1, paired = T, alternative
= "greater", : kann bei Bindungen keinen exakten p-Wert Berechnen
Warning in wilcox.test.default(werte_med2, werte_med1, paired = T, alternative
= "greater", : kann den exakten p-Wert bei Nullen nicht berechnen
# Exact wird auf TRUE gesetzt, wenn weniger als 40 Beobachtungen 
# vorliegen, wie es hier der Fall ist

# Ergebnis ausgeben
test_ergebnis_npar

    Wilcoxon signed rank test with continuity correction

data:  werte_med2 and werte_med1
V = 45, p-value = 0.004545
alternative hypothesis: true location shift is greater than 0

Auswertung: Auch der nicht-parametrische Test bestätigt das Ergebnis. Der p-Wert liegt wieder weit unter dem Signifikanzniveau von 0.05.

Ergebnisformulierung: Der Mann-Whitney-U-Test bestätigt, dass die zusätzliche Schlafdauer bei Medikament 2 signifikant höher ist als bei Medikament 1 (V = 45, p = .005).

Dieses Ergebnis hatten wir bereits erwartet nach den ersten visuellen Vergleichen und der Betrachtung der berechneten Differenzen. Wir könnten die Forschungsfrage noch etwas verschärfen mit diesem Vorwissen, dass Medikament 2 offensichtlich besser wirkt als Medikament 1.

fortgeschrittene Forschungsfrage: Ist der Schlaf durch Schlafmittel 2 um mindestens eine Stunde mehr verlängert als durch Schlafmittel 1?

Dazu führen wir wieder einen einseitigen Wilcoxon-Test aus und erweitern den Funktionsaufruf um das Argument mu = 1 (standardmäßig ist mu = 0). Den Parameter kann man ebenso bei t.test verwenden.

# Durchführung des Wilcoxon-Vorzeichen-Rang-Test mit 'mu = 1'
wilcox.test(werte_med2, werte_med1,
            paired = T, alternative = "greater", exact = T, mu = 1)
Warning in wilcox.test.default(werte_med2, werte_med1, paired = T, alternative
= "greater", : kann bei Bindungen keinen exakten p-Wert Berechnen

    Wilcoxon signed rank test with continuity correction

data:  werte_med2 and werte_med1
V = 44.5, p-value = 0.04609
alternative hypothesis: true location shift is greater than 1

Auch für diese Hypothese lehnt der Test auf einem Signifikanzniveau von 0.05 ab, allerdings deutlich knapper.

Ergebnisformulierung: Der Test bestätigt statistisch signifikant, dass die Verlängerung der Schlafdauer durch Medikament 2 im Vergleich zu Medikament 1 mehr als eine Stunde beträgt (V = 44.5, p = .046).

Hinweis: Der p-Wert ist mit 0.02 größer als vorher (0.005), da die Hürde (“mehr als 1 Stunde”) schwerer zu nehmen ist als die ursprüngliche Hürde (“mehr als 0 Stunden”).


Übungsabschnitt

Szenario: Ein Fitnessstudio bietet einen 4-wöchigen Kurs an. Wir haben das Gewicht von 9 Teilnehmenden vor und nach dem Kurs gemessen.

Frage: Hat der Kurs zu einer signifikanten Gewichtsabnahme geführt?

Der Datensatz heißt diaet_daten.

Wende nun das gerade dargestellte Prozedere auf diese Frage an.

1. Datenaufbereitung

# Prüfe die Struktur der Variablen
str(diaet_daten)
'data.frame':   9 obs. of  3 variables:
 $ Teilnehmer: Factor w/ 9 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9
 $ Vorher    : num  85 92 78 88 75 102 95 82 79
 $ Nachher   : num  83 89 79 84 73 98 92 80 78

Die Variable Teilnehmer ist ein Faktor. Die anderen beiden Variablen enthalten die gemessenen Gewichte Vorher bzw. Nachher. Wir haben also keine einzelne numerische Variable, die durch zwei Faktor-Variablen die Teilnehmenden klassifiziert, wie im Beispiel davor, sondern einen Faktor mit zwei numerischen Variablen. Den Extraktionsschritt aus Abschnitt 1b können wir somit überspringen und direkt aus den Spalten die Differenzen bilden. Vorab aber noch die Prüfung auf Fehlwerte.

# Gibt es Fehlwerte? 
table(is.na(diaet_daten))

FALSE 
   27 
# Nein, der Datensatz ist sauber.

# Abschließender Schritt der Datenaufbereitung --> Differenz berechnen (Nachher - Vorher)
diff_gewicht <- diaet_daten$Nachher - diaet_daten$Vorher

summary(diff_gewicht)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 -4.000  -3.000  -2.000  -2.222  -2.000   1.000 

Logik-Check: Wenn man abnimmt, sollte das Gewicht “Nachher” kleiner sein als “Vorher”. Die Differenz (Nachher - Vorher) sollte also negativ sein. Das Summary zeigt: Fast alle Werte liegen hier unter 0. Das sieht gut aus für das Fitnessstudio.

2. Explorative Datenanalyse

# Visualisierung --> Boxplot der Differenzen
boxplot(diff_gewicht, 
        main = "Gewichtsveränderung", 
        ylab = "kg", col = "lightgreen")
abline(h = 0, col = "red", lty = 2) # Rote Linie bei 0 zur Orientierung

Visuelle Interpretation: Bis auf ein Wert liegen viele unterhalb der roten Linie, heißt fast alle Teilnehmenden haben in dem Zeitintervall abgenommen.

3. Voraussetzungen prüfen

# grafisch --> Histogramm der Differenzen
hist(diff_gewicht, main = "Verteilung der Differenzen", 
     xlab = "Differenz (kg)", col = "lightgreen")

# Die Verteilung sieht eher rechtsschief aus. Der Stichprobenumfang ist aber auch sehr gering, worunter die Aussagekraft des Histogramms deutlich leidet.
# statistisch --> Shapiro-Wilk-Test auf die Differenzen
shapiro.test(diff_gewicht)

    Shapiro-Wilk normality test

data:  diff_gewicht
W = 0.90611, p-value = 0.2896

Auswertung: Der p-Wert beträgt 0.29. Das ist deutlich größer als 0.05. Das bedeutet: Wir behalten die Nullhypothese bei. Die Differenzen sind normalverteilt. Damit sind die Voraussetzungen für den parametrischen t-Test erfüllt und wir ziehen ihn der nicht-parametrische Alternative vor.

4. parametrischen oder nicht-parametrischen Mittelwertsvergleich

Da wir eine gerichtete Hypothese haben (“Hat der Kurs zum abnehmen geführt?”), testen wir einseitig. Wir übergeben der Funktion erst die Werte von Nachher und danach von Vorher. Hypothese: Nachher < Vorher. Daher nutzen wir alternative = "less".

# t-Test mit paired = TRUE und alternative = "less"
t.test(diaet_daten$Nachher, diaet_daten$Vorher, 
       paired = TRUE, 
       alternative = "less")

    Paired t-test

data:  diaet_daten$Nachher and diaet_daten$Vorher
t = -4.264, df = 8, p-value = 0.001373
alternative hypothesis: true mean difference is less than 0
95 percent confidence interval:
      -Inf -1.253105
sample estimates:
mean difference 
      -2.222222 

Ergebnisformulierung: Die Analyse bestätigt eine signifikante Gewichtsabnahme durch den Kurs. Das Gewicht nach dem Kurs ist signifikant niedriger als davor (t = -4.26, p = .001). Im Durchschnitt nahmen die Teilnehmer ca. 2,2 kg ab.