Gemischte Modell-Mehrfachvergleiche für die Interaktion zwischen kontinuierlichem und kategorialem Prädiktor

11

Ich möchte verwenden lme4, um eine Regression mit gemischten Effekten anzupassen und multcompdie paarweisen Vergleiche zu berechnen. Ich habe einen komplexen Datensatz mit mehreren kontinuierlichen und kategorialen Prädiktoren, aber meine Frage kann am Beispiel des integrierten ChickWeightDatensatzes demonstriert werden :

m <- lmer(weight ~ Time * Diet + (1 | Chick), data=ChickWeight, REML=F)

Timeist kontinuierlich und Dietist kategorisch (4 Stufen) und es gibt mehrere Küken pro Diät. Alle Küken begannen mit ungefähr dem gleichen Gewicht, aber ihre Ernährung (kann) ihre Wachstumsrate beeinflussen, so dass die DietAbschnitte (mehr oder weniger) gleich sein sollten, aber die Steigungen können unterschiedlich sein. Ich kann die paarweisen Vergleiche für den Abfangeffekt Dietwie folgt erhalten:

summary(glht(m, linfct=mcp(Diet = "Tukey")))

und tatsächlich unterscheiden sie sich nicht wesentlich, aber wie kann ich den analogen Test für den Time:DietEffekt durchführen? Das einfache Einfügen des Interaktionsterms mcpführt zu einem Fehler:

summary(glht(m, linfct=mcp('Time:Diet' = "Tukey")))
Error in summary(glht(m, linfct = mcp(`Time:Diet` = "Tukey"))) : 
  error in evaluating the argument 'object' in selecting a method for function
 'summary': Error in mcp2matrix(model, linfct = linfct) : 
Variable(s) ‘Time:Diet’ have been specified in ‘linfct’ but cannot be found in ‘model’!
r mixed-model multiple-comparisons Dan M.
quelle
Es hat Time*Diet, was nur eine Vereinfachung von ist Time + Diet + Time:Diet. Verwenden anova(m)oder summary(m)bestätigen, dass sich der Interaktionsterm im Modell befindet.
Dan M.

Antworten:

8

Standardmäßig wird lmerdie Referenzstufe eines kategorialen Prädiktors als Basislinie behandelt und die Parameter für die anderen Stufen geschätzt. Sie erhalten also einige paarweise Vergleiche in der Standardausgabe, und Sie können die anderen erhalten, indem Sie releveleinen neuen Referenzpegel definieren und das Modell neu anpassen. Dies hat den Vorteil, dass Sie Modellvergleiche oder MCMC verwenden können, um die p-Werte zu erhalten, korrigiert jedoch nicht mehrere Vergleiche (obwohl Sie anschließend Ihre eigene Korrektur anwenden können).

Zur Verwendung multcompmüssen Sie die Kontrastmatrix definieren. Jede Zeile in der Kontrastmatrix repräsentiert Gewichte für die Effekte, die Sie in der Standardmodellausgabe erhalten, beginnend mit dem Intercept. Wenn Sie also einen Effekt wünschen, der bereits in der Basisausgabe enthalten ist, setzen Sie einfach eine "1" an die Position, die diesem Effekt entspricht. Da sich die Parameterschätzungen auf ein gemeinsames Referenzniveau beziehen, können Sie Vergleiche zwischen zwei anderen Niveaus erhalten, indem Sie das Gewicht eines auf "-1" und des anderen auf "1" setzen. So würde das für die Time:DietBegriffe im ChickWeightBeispiel funktionieren :

contrast.matrix <- rbind("Time:Diet1 vs. Time:Diet2" =  c(0, 0, 0, 0, 0, 1, 0, 0),
                           "Time:Diet1 vs. Time:Diet3" =  c(0, 0, 0, 0, 0, 0, 1, 0),
                           "Time:Diet1 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, 0, 0, 1),
                           "Time:Diet2 vs. Time:Diet3" =  c(0, 0, 0, 0, 0, -1, 1, 0),
                           "Time:Diet2 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, -1, 0, 1),
                           "Time:Diet3 vs. Time:Diet4" =  c(0, 0, 0, 0, 0, 0, -1, 1))
summary(glht(m, contrast.matrix))

Vorsichtsmaßnahme: Dieser Ansatz scheint die normale Näherung zu verwenden, um p-Werte zu erhalten, was etwas antikonservativ ist, und wendet dann eine Korrektur für mehrere Vergleiche an. Das Ergebnis ist, dass Sie mit dieser Methode einfach auf so viele paarweise Parameterschätzungen und Standardfehler zugreifen können, wie Sie möchten. Die p-Werte können jedoch Ihren Wünschen entsprechen oder nicht.

(Danke an Scott Jackson von r-ling-lang-L für die Hilfe dabei)

Dan M.
quelle