Verwenden von lmer für lineare Mischeffektmodelle mit wiederholten Messungen

41

EDIT 2: Ursprünglich dachte ich, ich müsste eine Zweifaktor-ANOVA mit wiederholten Messungen für einen Faktor durchführen, aber jetzt denke ich, dass ein lineares Mischeffektmodell für meine Daten besser funktioniert. Ich glaube, ich weiß fast, was passieren muss, aber ich bin immer noch durch einige Punkte verwirrt.

Die Experimente, die ich analysieren muss, sehen folgendermaßen aus:

  • Die Probanden wurden einer von mehreren Behandlungsgruppen zugeordnet
  • Die Messungen jedes Probanden wurden an mehreren Tagen durchgeführt
  • Damit:
    • Das Subjekt ist in der Behandlung verschachtelt
    • Die Behandlung wird mit dem Tag gekreuzt

(Jedes Subjekt ist nur einer Behandlung zugeordnet, und an jedem Tag werden Messungen an jedem Subjekt vorgenommen.)

Mein Datensatz enthält folgende Informationen:

  • Betreff = Sperrfaktor (Zufallsfaktor)
  • Tag = innerhalb des Subjekts oder Faktor für wiederholte Messungen (fester Faktor)
  • Behandlung = zwischen Subjektfaktor (fester Faktor)
  • Obs = gemessene (abhängige) Variable

UPDATE OK, also habe ich mit einem Statistiker gesprochen, aber er ist ein SAS-Benutzer. Er meint, dass das Modell sein sollte:

Behandlung + Tag + Subjekt (Behandlung) + Tag * Subjekt (Behandlung)

Offensichtlich unterscheidet sich seine Notation von der R-Syntax, aber dieses Modell soll berücksichtigen:

  • Behandlung (fix)
  • Tag (fest)
  • die Behandlung * Day Interaktion
  • Betreff in Behandlung verschachtelt (zufällig)
  • Tag gekreuzt mit "Subject within Treatment" (zufällig)

Ist dies also die richtige Syntax?

m4 <- lmer(Obs~Treatment*Day + (1+Treatment/Subject) + (1+Day*Treatment/Subject), mydata)

Ich bin besonders besorgt darüber, ob der Tag, der mit dem Teil "Thema innerhalb der Behandlung" gekreuzt wurde, richtig ist. Ist jemand, der mit SAS vertraut ist oder der zuversichtlich ist, dass er versteht, was in seinem Modell vor sich geht, und in der Lage ist, zu kommentieren, ob mein trauriger Versuch der R-Syntax passt?

Hier sind meine früheren Versuche, ein Modell zu erstellen und eine Syntax zu schreiben (siehe Antworten und Kommentare):

m1 <- lmer(Obs ~ Treatment * Day + (1 | Subject), mydata)

Wie gehe ich damit um, dass das Thema in der Behandlung verschachtelt ist? Wie unterscheidet m1sich von:

m2 <- lmer(Obs ~ Treatment * Day + (Treatment|Subject), mydata)
m3 <- lmer(Obs ~ Treatment * Day + (Treatment:Subject), mydata)

und sind m2und m3gleichwertig (und wenn nicht, warum)?

Muss ich nlme anstelle von lme4 verwenden, um die Korrelationsstruktur (wie correlation = corAR1) anzugeben ? Gemäß den Wiederholten Messungen ist für eine Analyse mit wiederholten Messungen eines Faktors die Kovarianzstruktur (die Art der Korrelationen zwischen Messungen desselben Subjekts) wichtig.

Als ich versuchte, eine ANOVA mit wiederholten Messungen durchzuführen, hatte ich mich für eine SS vom Typ II entschieden. Ist dies immer noch relevant und wenn ja, wie gehe ich vor, um dies zu spezifizieren?

Hier ist ein Beispiel, wie die Daten aussehen:

mydata <- data.frame(
  Subject  = c(13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 
               34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 
               19, 20, 21, 22, 23, 24, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
               40, 62, 63, 64, 65, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 
               29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 62, 63, 64, 65), 
  Day       = c(rep(c("Day1", "Day3", "Day6"), each=28)), 
  Treatment = c(rep(c("B", "A", "C", "B", "C", "A", "A", "B", "A", "C", "B", "C", 
                      "A", "A", "B", "A", "C", "B", "C", "A", "A"), each = 4)), 
  Obs       = c(6.472687, 7.017110, 6.200715, 6.613928, 6.829968, 7.387583, 7.367293, 
                8.018853, 7.527408, 6.746739, 7.296910, 6.983360, 6.816621, 6.571689, 
                5.911261, 6.954988, 7.624122, 7.669865, 7.676225, 7.263593, 7.704737, 
                7.328716, 7.295610, 5.964180, 6.880814, 6.926342, 6.926342, 7.562293, 
                6.677607, 7.023526, 6.441864, 7.020875, 7.478931, 7.495336, 7.427709, 
                7.633020, 7.382091, 7.359731, 7.285889, 7.496863, 6.632403, 6.171196, 
                6.306012, 7.253833, 7.594852, 6.915225, 7.220147, 7.298227, 7.573612, 
                7.366550, 7.560513, 7.289078, 7.287802, 7.155336, 7.394452, 7.465383, 
                6.976048, 7.222966, 6.584153, 7.013223, 7.569905, 7.459185, 7.504068, 
                7.801867, 7.598728, 7.475841, 7.511873, 7.518384, 6.618589, 5.854754, 
                6.125749, 6.962720, 7.540600, 7.379861, 7.344189, 7.362815, 7.805802, 
                7.764172, 7.789844, 7.616437, NA, NA, NA, NA))
phosphoreliert
quelle

Antworten:

18

Ich denke, dass Ihr Ansatz richtig ist. Das Modell m1gibt für jedes Subjekt einen eigenen Schnittpunkt an. Das Modell m2fügt für jedes Motiv eine eigene Neigung hinzu. Ihre Piste verläuft über mehrere Tage, da die Probanden nur an einer Behandlungsgruppe teilnehmen. Wenn Sie das Modell m2wie folgt schreiben, ist es offensichtlicher, dass Sie für jedes Motiv einen eigenen Schnittpunkt und eine eigene Steigung modellieren

m2 <- lmer(Obs ~ Treatment * Day + (1+Day|Subject), mydata)

Dies ist äquivalent zu:

m2 <- lmer(Obs ~ Treatment + Day + Treatment:Day + (1+Day|Subject), mydata)

Dh die Hauptwirkungen der Behandlung, des Tages und der Wechselwirkung zwischen beiden.

Ich denke, Sie brauchen sich keine Gedanken über das Verschachteln zu machen, solange Sie die Probanden-IDs in den Behandlungsgruppen nicht wiederholen. Welches Modell richtig ist, hängt wirklich von Ihrer Forschungsfrage ab. Gibt es Grund zu der Annahme, dass die Neigungen der Probanden zusätzlich zum Behandlungseffekt variieren? Sie können beide Modelle ausführen und mit vergleichen, anova(m1,m2)um festzustellen, ob die Daten eines von beiden unterstützen.

Ich bin nicht sicher, was Sie mit Modell ausdrücken möchten m3? Die Verschachtelungssyntax verwendet eine /, z (1|group/subgroup).

Ich glaube nicht, dass Sie sich wegen der Autokorrelation mit einer so geringen Anzahl von Zeitpunkten Sorgen machen müssen.

user12719
quelle
Das ist nicht richtig. Die Behandlung ist eine Variable der Ebene 2 und kann nicht in Subjekten verschachtelt werden.
Patrick Coulombe
Informationen zur Autokorrelation und zur Anzahl der Zeitpunkte: In diesem Beispiel werden nur drei Daten angezeigt. Meine realen Daten enthalten jedoch Beobachtungen an acht verschiedenen Tagen. Daher denke ich, dass dies wahrscheinlich ein Problem sein wird. Irgendwelche Ideen, wie man das einfügt?
Phosphorelated
1
Außerdem bin ich jetzt ziemlich verwirrt über das Verschachteln; unterscheidet sich (1 + Behandlung | Subjekt) von (1 + Behandlung / Subjekt)? Was bedeutet das "|" im Klartext? Ich verstehe die Erklärungen, die ich gelesen habe, nicht.
Phosphorelated
Hallo. Was ist hier ein "separater Hang für jedes Thema"? weil subject eine Faktorvariable ist, keine stetige Variable.
22.
12

Ich fühle mich nicht wohl genug, um Ihr Problem mit automatisch korrelierten Fehlern (oder die verschiedenen Implementierungen in lme4 vs. nlme) zu kommentieren, aber ich kann mit dem Rest sprechen.

Ihr Modell m1ist ein Zufallsschnittstellenmodell, bei dem Sie die Interaktion zwischen Behandlung und Tag auf verschiedenen Ebenen berücksichtigt haben (die Wirkung des Tages kann zwischen den Behandlungsgruppen variieren). Um zu ermöglichen, dass sich die Änderung über die Zeit zwischen den Teilnehmern unterscheidet (dh um einzelne Unterschiede in der Änderung über die Zeit explizit zu modellieren), müssen Sie auch zulassen, dass der Effekt des Tages zufällig ist . Dazu würden Sie Folgendes angeben:

m2 <- lmer(Obs ~ Day + Treatment + Day:Treatment + (Day | Subject), mydata)

In diesem Modell:

  • Der Abschnitt, wenn die vorhergesagte Punktzahl für die Behandlungsreferenzkategorie am Tag = 0 ist
  • Der Koeffizient für Tag ist die vorhergesagte zeitliche Änderung für jede Erhöhung der Tage um 1 Einheit für die Referenzkategorie der Behandlung
  • Die Koeffizienten für die beiden Dummy-Codes für die Behandlungsgruppen (automatisch von R erstellt) sind die vorhergesagte Differenz zwischen jeder verbleibenden Behandlungsgruppe und der Referenzkategorie am Tag = 0
  • Die Koeffizienten für die beiden Interaktionsterme sind der Unterschied in der Auswirkung der Zeit (Tag) auf die vorhergesagten Scores zwischen der Referenzkategorie und den verbleibenden Behandlungsgruppen.

Sowohl die Abschnitte als auch die Auswirkung des Tages auf die Punktzahl sind zufällig (jedes Subjekt darf an Tag = 0 eine andere vorhergesagte Punktzahl und eine andere lineare Änderung über die Zeit haben). Die Kovarianz zwischen Abschnitten und Steigungen wird ebenfalls modelliert (sie dürfen kovarieren).

Wie Sie sehen, ist die Interpretation der Koeffizienten für die beiden Dummy-Variablen an Tag = 0 gebunden. Sie werden Ihnen mitteilen, ob der prognostizierte Score am Tag = 0 für die Referenzkategorie signifikant von den beiden verbleibenden Behandlungsgruppen abweicht. Daher ist es wichtig, wo Sie sich entscheiden, Ihre Tagesvariable zu zentrieren. Wenn Sie am ersten Tag zentrieren, geben die Koeffizienten Aufschluss darüber, ob die vorhergesagte Punktzahl für die Referenzkategorie am ersten Tag erheblich von der vorhergesagten Punktzahl der beiden verbleibenden Gruppen abweicht. Auf diese Weise können Sie feststellen, ob es bereits Unterschiede zwischen den Gruppen gibt . Wenn Sie am 3. Tag zentrieren, geben die Koeffizienten an, ob die vorhergesagte Punktzahl für die Referenzkategorie am 3. Tag erreicht wurdeunterscheidet sich signifikant von der vorhergesagten Punktzahl der beiden verbleibenden Gruppen. Auf diese Weise können Sie am Ende des Eingriffs feststellen, ob es Unterschiede zwischen den Gruppen gibt .

Beachten Sie schließlich, dass die Probanden nicht in der Behandlung verschachtelt sind. Ihre drei Behandlungen sind keine zufälligen Ebenen einer Population von Ebenen, auf die Sie Ihre Ergebnisse verallgemeinern möchten. Wie Sie bereits erwähnt haben, sind Ihre Ebenen festgelegt, und Sie möchten Ihre Ergebnisse nur auf diese Ebenen verallgemeinern. (Ganz zu schweigen davon, dass Sie keine Mehrebenenmodellierung verwenden sollten, wenn Sie nur 3 Einheiten der oberen Ebene haben; siehe Maas & Hox, 2005.) Stattdessen ist die Behandlung ein Prädiktor der Ebene 2, dh ein Prädiktor, der über Tage hinweg einen einzelnen Wert annimmt (Level-1-Einheiten) für jedes Fach. Daher wird es lediglich als Prädiktor in Ihr Modell aufgenommen.

Referenz:
Maas, CJM & Hox, JJ (2005). Ausreichende Stichprobengrößen für die mehrstufige Modellierung. Methodik: European Journal of Research Methods für die Verhaltens- und Sozialwissenschaften , 1 , 86-92.

Patrick Coulombe
quelle
1
Sie kann von lmer nicht geschätzt werden, da die Anzahl der zufälligen Effekte und die verbleibende Varianz wahrscheinlich nicht identifizierbar sind.
Shuguang,
Die Formelstruktur in der Antwort ist korrekt. Um den von @Shuguang erwähnten Fehler zu überschreiben, müssten Sie hinzufügen ...,control=lmerControl(check.nobs.vs.nRE="ignore"). Weitere Erklärungen von Ben Bolker finden Sie unter diesem Link .
NiuBiBang
Schöne Erklärungen. Können Sie uns bitte etwas näher erläutern, warum "Probanden nicht in der Behandlung verschachtelt sind" und warum Sie keinen Fehlerbegriff + (Behandlung | Betreff) erstellen und warum nicht nur (1 | Betreff) oder sogar (1 | Behandlung * Tag) )?
22.
Technisch Sie könnten Nest Themen innerhalb Behandlung, aber wenn der Prädiktor ist eine , die die gleiche , egal wäre , wie oft man das Experiment lief, sollte es eine feste (nicht zufällig) Wirkung. Faktoren, die sich bei jeder Durchführung des Experiments ändern würden , wie z. B. die individuellen Merkmale des Probanden - z. B. der Anfangswert oder die eigenwillige Reaktion auf Änderungen der Behandlung im Zeitverlauf -, sind zufällige Effekte. (1 + Day|Subject)bedeutet ein Modell mit zufälligen Steigungen, bei dem der Anfangswert (Intercept) und die Änderungsrate des Ergebnisses für jedes Subjekt unterschiedlich sind.
Llewmills