Mehrfachvergleiche an einem Mixed-Effects-Modell

31

Ich versuche, einige Daten mit einem gemischten Effektmodell zu analysieren. Die von mir gesammelten Daten repräsentieren das Gewicht einiger Jungtiere unterschiedlichen Genotyps im Zeitverlauf.

Ich verwende den hier vorgeschlagenen Ansatz: https://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/

Insbesondere verwende ich Lösung # 2

Also ich habe sowas

require(nlme)
model <- lme(weight ~ time * Genotype, random = ~1|Animal/time, 
         data=weights)    
av <- anova(model)

Jetzt möchte ich einige mehrfache Vergleiche haben. Mit multcompkann ich:

require(multcomp)
comp.geno <- glht(model, linfct=mcp(Genotype="Tukey"))
print(summary(comp.geno))

Und natürlich könnte ich das auch mit der Zeit tun.

Ich habe zwei Fragen:

  1. Wie kann ich mcpdie Interaktion zwischen Zeit und Genotyp anzeigen?
  2. Wenn ich renne, glhtbekomme ich folgende Warnung:

    covariate interactions found -- default contrast might be inappropriate

    Was heißt das? Kann ich das ignorieren? Oder was soll ich tun, um das zu vermeiden?

BEARBEITEN: Ich habe dieses PDF gefunden , in dem steht:

Da es in diesem Fall nicht möglich ist, die interessierenden Parameter automatisch zu bestimmen, generiert mcp () in multcomp standardmäßig nur Vergleiche für die Haupteffekte, wobei Kovariaten und Interaktionen ignoriert werden . Seit Version 1.1-2 kann angegeben werden, dass über Interaktionsterme und Kovariaten mit den Argumenten interact_average = TRUE und covariate_average = TRUE gemittelt werden soll, während Versionen älter als 1.0-0 automatisch über Interaktionsterme gemittelt werden. Wir empfehlen den Benutzern jedoch, die gewünschten Kontraste manuell zu schreiben.Dies sollte immer dann erfolgen, wenn Zweifel an der Messung der Standardkontraste bestehen, was normalerweise bei Modellen mit Interaktionsbegriffen höherer Ordnung der Fall ist. Wir verweisen auf Hsu (1996), Kapitel ~ 7, und Searle (1971), Kapitel ~ 7.3, für weitere Diskussionen und Beispiele zu diesem Thema.

Ich habe keinen Zugang zu diesen Büchern, aber vielleicht hat jemand hier?

nico
quelle
Schauen Sie sich InvivoStat ( invivostat.co.uk ) an, es sollte die Art von Analyse durchführen, nach der Sie suchen

Antworten:

29

Wenn timeund Genotypebeide kategoriale Prädiktoren sind, und Sie daran interessiert sind, alle Zeit / Genotyp-Paare miteinander zu vergleichen, können Sie einfach eine Interaktionsvariable erstellen und Tukey-Kontraste verwenden:

weights$TimeGeno <- interaction(weigths$Time, weights$Geno)
model <- lme(weight ~ TimeGeno, random = ~1|Animal/time, data=weights) 
comp.timegeno <- glht(model, linfct=mcp(TimeGeno="Tukey")) 

Wenn Sie an anderen Kontrasten interessiert sind, können Sie die Tatsache nutzen, dass das linfctArgument eine Koeffizientenmatrix für die Kontraste annehmen kann. Auf diese Weise können Sie genau die Vergleiche einrichten, die Sie möchten.

BEARBEITEN

In den Kommentaren taucht Bedenken auf, dass sich das mit dem TimeGenoPrädiktor ausgestattete Modell von dem mit dem Time * GenotypePrädiktor ausgestatteten Originalmodell unterscheidet . Dies ist nicht der Fall , die Modelle sind gleichwertig. Der einzige Unterschied besteht in der Parametrisierung der festen Effekte, die so eingerichtet ist, dass die Verwendung der glhtFunktion vereinfacht wird .

Ich habe einen der integrierten Datasets (es hat Diet anstelle von Genotype) verwendet, um zu demonstrieren, dass die beiden Ansätze die gleiche Wahrscheinlichkeit, die gleichen vorhergesagten Werte usw. haben:

> # extract a subset of a built-in dataset for the example
> data(BodyWeight)
> ex <- as.data.frame(subset(BodyWeight, Time %in% c(1, 22, 44)))
> ex$Time <- factor(ex$Time)
> 
> #create interaction variable
> ex$TimeDiet <- interaction(ex$Time, ex$Diet)
    > 
    > model1 <- lme(weight ~ Time * Diet, random = ~1|Rat/Time,  data=ex)    
    > model2 <- lme(weight ~ TimeDiet, random = ~1|Rat/Time, data=ex)    
    > 
    > # the degrees of freedom, AIC, BIC, log-likelihood are all the same 
    > anova(model1, model2)
           Model df      AIC      BIC    logLik
    model1     1 12 367.4266 387.3893 -171.7133
    model2     2 12 367.4266 387.3893 -171.7133
    Warning message:
    In anova.lme(model1, model2) :
      fitted objects with different fixed effects. REML comparisons are not meaningful.
    > 
    > # the second model collapses the main and interaction effects of the first model
    > anova(model1)
                numDF denDF   F-value p-value
    (Intercept)     1    26 1719.5059  <.0001
    Time            2    26   28.9986  <.0001
    Diet            2    13   85.3659  <.0001
    Time:Diet       4    26    1.7610  0.1671
    > anova(model2)
                numDF denDF   F-value p-value
    (Intercept)     1    24 1719.5059  <.0001
    TimeDiet        8    24   29.4716  <.0001
    > 
    > # they give the same predicted values
    > newdata <- expand.grid(Time=levels(ex$Time), Diet=levels(ex$Diet))
    > newdata$TimeDiet <- interaction(newdata$Time, newdata$Diet)
> newdata$pred1 <- predict(model1, newdata=newdata, level=0)
    > newdata$pred2 <- predict(model2, newdata=newdata, level=0)
> newdata
  Time Diet TimeDiet   pred1   pred2
1    1    1      1.1 250.625 250.625
2   22    1     22.1 261.875 261.875
3   44    1     44.1 267.250 267.250
4    1    2      1.2 453.750 453.750
5   22    2     22.2 475.000 475.000
6   44    2     44.2 488.750 488.750
7    1    3      1.3 508.750 508.750
8   22    3     22.3 518.250 518.250
9   44    3     44.3 530.000 530.000

Der einzige Unterschied ist, dass die Hypothesen leicht zu testen sind. Zum Beispiel ist es im ersten Modell einfach zu testen, ob die beiden Prädiktoren interagieren, im zweiten Modell gibt es keinen expliziten Test dafür. Andererseits ist die gemeinsame Wirkung der beiden Prädiktoren im zweiten Modell leicht zu testen, nicht jedoch im ersten. Die anderen Hypothesen sind überprüfbar, es ist nur mehr Arbeit, sie aufzustellen.

Aniko
quelle
3
glhtverwendet die im lme-Modell angegebenen Freiheitsgrade. Ich bin mir nicht sicher, ob diese Freiheitsgrade angemessen sind ...?
Stéphane Laurent
2
Ich bin auch gespannt, wie das am besten geht. Dieser Ansatz führt jedoch zu Effekten eines anderen Modells - eines, das im Wesentlichen nur eine Interaktion testet. Das zweite Modell enthält überhaupt nicht die beiden möglichen Haupteffekte. Dies scheint keine geeignete Methode zur Überprüfung der Auswirkungen im ersten Modell zu sein.
Marcus Morrisey
@Aniko, ich habe gerade darüber nachgedacht, zwei kategoriale Variablen in eine zu kombinieren, aber ich habe gezögert, weil nur eine der Variablen im Betreff liegt, die andere zwischen. Können Sie bestätigen, dass dies nicht zutrifft? Mir ist aufgefallen, dass Sie in dem Beispiel behalten, Animal/timewas jetzt nicht einer der Faktoren ist. Ist mir das wirklich understandso?
toto_tico
@toto_tico, ich habe die Antwort bearbeitet, um zu zeigen, dass das zweite Modell dem ersten Modell entspricht.
Aniko
3
@toto_tico, ich habe dir ein reproduzierbares Beispiel gegeben. Warum versuchst du nicht zu all.equal(resid(model1), resid(model2))sehen, dass sie gleich sind, bevor du etwas anderes errätst? Nur die Parametrisierung der festen Effekte ist unterschiedlich. TimeDietist kein reiner Interaktionsbegriff und nicht gleichbedeutend mit Time:Diet, sondern vielmehr mit Time + Diet + Time:Diet.
Aniko