1 Datenstruktur: abhängige Daten

Wir beginnen die Erarbeitung anhand eines Minimalbeispiels, und setzen die dort gewonnenen Erkenntnisse dann im Datensatz MODALIS um.

Die Struktur abhängiger Daten stellt sich in vielen Datensätzen (und Softwarepaketen wie beispielsweise SPSS) so dar, dass wiederholte Messungen in Spalten nebeneinander stehen (müssen). Diese Anordnung entspricht einem “breiten” Datenformat. Unser Minimalbeispiel verdeutlicht dies: pro Proband gibt es eine Zeile, alle Wiederholunsgmessungen derselben Person stehen in Spalten.

load("schlaf_breit.RData")
head(schlaf.b)

  Proband -20h -24h -28h +1h
1     P01   19   13   12  15
2     P02   21   15   14  17
3     P03   17   18   15  15
4     P04   20   17   16  17
5     P05   23   16   13  16

boxplot(schlaf.b[,2:5])

Der Datensatz enthält Scores aus Konzentrationstests (höherer Wert \(\hat{=}\) bessere Leistung) über wiederholte Messungen nach Schlafentzug. In der improvisierten und Grafik (wir werden weiter unten eine bessere Grafik dazu erstellen!) ist das zu erwartenden Bild erkennbar, dass die Leistung mit zunehmendem Schlafentzug sinkt und nach einer Stunde Erholungsschlaf wieder leicht ansteigt.

Bei R ist eine solche Datenanordnung für abhängige Daten nicht erforderlich – im Gegenteil. Es wäre sogar vorteilhaft, wenn diese Daten in eine “lange” Schreibweise überführt würden, da die meisten Funktionen und Pakete in R eine solche, faktororientierte Schreibweise voraussetzen. Zunächst ist also die Aufgabe, die Daten aus einem breiten in ein langes Format zu überführen.

Im Minimalbeispiel könnte man das händisch tun, in größeren Datensätzen wie MODALIS wird das aber schnell zu unübersichtlich und vor allem fehleranfällig. Um ein breites in ein langes Format zu überführen, steht die Funktion melt() aus dem Paket reshape2 zur Verfügung. Man erspart man sich damit händische Tipparbeit und vermeidet Fehler. Eine sehr schöne Anleitung zum Umwandeln von breit zu lang und umgekehrt findet sich hier.

library(reshape2)
schlaf.l <- melt(schlaf.b, id.vars = "Proband",
                 measure.vars = c("-20h", "-24h", "-28h", "+1h"),
                 variable.name = "Schlafentzug",
                 value.name = "Score")

head(schlaf.l, n=12)

   Proband Schlafentzug Score
1      P01         -20h    19
2      P02         -20h    21
3      P03         -20h    17
4      P04         -20h    20
5      P05         -20h    23
6      P01         -24h    13
7      P02         -24h    15
8      P03         -24h    18
9      P04         -24h    17
10     P05         -24h    16
11     P01         -28h    12
12     P02         -28h    14

Wir sehen, dass die ursprüngliche Anordnung durch das lange Datenformat ersetzt wird, bei welchem die Spalte Schlafentzug die Faktorstufen der Messwiederholung enthält, und die Spalte Score den eigentlichen Messwert. Daraus resultiert, dass jeder Proband mehrere Datenzeilen in diesem data.frame erhält, nämlich eine für jede Messung. Eine farbcodierte Übersicht der Messdaten auf Ebene der individuellen Verläufe des Scores erhalten wir mit ggplot2:

library(ggplot2)

P <- ggplot(schlaf.l, aes(x=Schlafentzug, y=Score, group=Proband, color=Proband))+
  geom_point(shape=19, size=2)+
  geom_line()+
  theme_bw()
plot(P)

Achtung

Damit das funktioniert, muss man das Attribut group im intialen Mapping (aes) setzen. Dadurch wird definiert, welche der Datenpunkte jeweils mit Linien verbunden werden sollen.

2 Hypothesenprüfung: Schlafentzug verursacht Leistungsabfall

In bisherigen ALMs sind wir davon ausgegangen, dass wir die Streuung innerhalb der Gruppen nicht erklären können. Daher wurden diese Streuungen immer als “Error” modelliert und als unerklärte Varianz aufgefasst. In Datensätzen mit abhängigen Daten ist dies aber nicht mehr gleichermaßen der Fall:

P2 <- ggplot(schlaf.l, aes(x=Proband, y=Score, fill=Proband))+
  geom_boxplot()+
  theme_bw()
plot(P2)

Wir sehen, dass sich die Versuchspersonen systematisch voneinander unterscheiden. So erreicht P01 tendenziell eher schlechterere Werte als beispielsweise P04. Dies kann und muss man berücksichtigen, wir definieren diese individuelle Eigenschaft der Probanden als “random effect”. Das bedeutet, dass wir die mittlere Ausprägung der Scores einer Person insofern berücksichtigen, dass wir sie als Fehler, den eine bestimmte Person einbringt, trennen von Fehlern, die wir gar nicht kontrollieren können. Bei den Abweichungen der mittleren Scores über alle Messzeitpunkte zwischen den Probanden haben wir zumindest die Erklärung parat, dass dies eine individuelle Eigenschaft sein muss, die wir jedoch (aufgrund der Zufallsstichprobe) nicht kontrollieren, sondern bestenfalls von anderen, zufälligen Fehlern separieren können.

Das resultiert darin, dass wir jedem Proband unterstellen, dass der Effekt, den Schlafentzug auf die Scores hat, gleich ist, aber dieser Effekt auf unterschiedlichen Scoreniveaus stattfindet. Die Liniengrafik oben zeigt dies dadurch, dass viele parallele Linienverläufe sichtbar werden. Im Sinne eines linearen Modells wäre das ein konstanter Slope (d. h. eine konstante Steigung), aber ein variabler (zufälliger, “random”) Intercept. Ein solches random intercept Modell generieren wir mit dem Paket lme4 und der darin enthaltenen Funktion lmer(). Die darin enthaltene Notation (1|Proband) bedeutet, dass wir jedem Proband einen eigenen Intercept (Konstante im Modell) zugestehen, welcher durch die 1 repräsentiert ist.

library(lme4)

Lade nötiges Paket: Matrix

M2 <- lmer(Score ~ Schlafentzug + (1|Proband), data=schlaf.l)
summary(M2)

Linear mixed model fit by REML ['lmerMod']
Formula: Score ~ Schlafentzug + (1 | Proband)
   Data: schlaf.l

REML criterion at convergence: 69.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-1.8120 -0.5708 -0.1094  0.5440  1.7220 

Random effects:
 Groups   Name        Variance Std.Dev.
 Proband  (Intercept) 0.4583   0.677   
 Residual             2.5917   1.610   
Number of obs: 20, groups:  Proband, 5

Fixed effects:
                 Estimate Std. Error t value
(Intercept)        20.000      0.781  25.607
Schlafentzug-24h   -4.200      1.018  -4.125
Schlafentzug-28h   -6.000      1.018  -5.893
Schlafentzug+1h    -4.000      1.018  -3.929

Correlation of Fixed Effects:
            (Intr) Sch-24 Sch-28
Schlfntz-24 -0.652              
Schlfntz-28 -0.652  0.500       
Schlfntzg+1 -0.652  0.500  0.500

Mit diesem Modell gelingen Vergleiche zum Referenzzeitpunkt. Ergebnisse von Signifikanztests werden allerdings nicht ausgegeben. Wir werden in der kommenden Woche eine Methode kennen lernen, wie man (ähnlich einer ANOVA) testen kann, ob der Faktor Schlafentzug in der Gesamtschau einen signifikanten Effekt hat. Dieses Verfahren ist sehr anschaulich in den Tutorials von Bodo Winter beschrieben.

Gruppenvergleiche zwischen den Folgezeitpunkten (also den unterschiedlichen Schlafentzugsstufen oberhalb von 20h) sind wiederum nur nach Hinzunahme eines Mehrgruppenvergleichs möglich. Zwei Varianten, die nicht nur für Mixed Effects Modelle funktionieren, sind mit den Paketen emmeans¹ und multcompverfügbar, die im Tutorial zu unabhängigen Gruppen durchgeführte Methode TukeyHSD() ist bei Mixed Effects Modellen nicht direkt anwendbar.

library(emmeans)
emmeans(M2, list(pairwise ~ Schlafentzug), adjust = "tukey")

$`emmeans of Schlafentzug`
 Schlafentzug emmean    SE df lower.CL upper.CL
 -20h           20.0 0.781 15     18.3     21.7
 -24h           15.8 0.781 15     14.1     17.5
 -28h           14.0 0.781 15     12.3     15.7
 +1h            16.0 0.781 15     14.3     17.7

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$`pairwise differences of Schlafentzug`
 1               estimate   SE df t.ratio p.value
 (-20h) - (-24h)      4.2 1.02 12   4.125  0.0067
 (-20h) - (-28h)      6.0 1.02 12   5.893  0.0004
 (-20h) - (+1h)       4.0 1.02 12   3.929  0.0094
 (-24h) - (-28h)      1.8 1.02 12   1.768  0.3340
 (-24h) - (+1h)      -0.2 1.02 12  -0.196  0.9972
 (-28h) - (+1h)      -2.0 1.02 12  -1.964  0.2540

Degrees-of-freedom method: kenward-roger 
P value adjustment: tukey method for comparing a family of 4 estimates

2.1 Zwischenfazit

Aus den Ergebnissen des Mehrfachvergleichs erkennt man die fehlenden Signifikanzen der Unterschiede zwischen -24h und -28h sowie -28h und +1h. Damit unterscheiden sich also nicht alle Messzeitpunkte voneinander, sondern lediglich die die paarweisen Vergleiche zum Ausgangslevel von -20h.

Die erforderlichen Schritte zum Hypothesentest zwischen beliebigen Messzeitpunkten sind:

Daten aufbereiten und in ein langes Format bringen
Mixed Effect Modell anpassen (random intercept oder random slopes Modell (siehe nächstes Tutorial))
Mehrgruppenvegleiche mit emmeans rechnen.

Stand: 2024-04-18

Fußnoten

emmeans erscheint das einfacher zu benutzende Paket, daher verwenden wir ausschließlich dieses in diesem Tutorial. In einigen Fällen war bei Nutzung des Pakets emmeans ein händisches Nachinstallieren von pbkrtest erforderlich. Am besten schon bei der Installation den Befehl install.packages("emmeans", dependencies=TRUE) nutzen.↩︎