Modele wielokrotne porównań mieszanych dla interakcji między predyktorem ciągłym i kategorycznym

11

Chciałbym użyć, lme4aby dopasować regresję efektów mieszanych i multcompobliczyć porównania parami. Mam złożony zestaw danych z wieloma ciągłymi i kategorycznymi predyktorami, ale moje pytanie można zademonstrować na podstawie wbudowanego ChickWeightzestawu danych jako przykładu:

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

Timejest ciągły i Dietma charakter kategoryczny (4 poziomy), a na dietę przypada wiele piskląt. Wszystkie pisklęta zaczynały mniej więcej taką samą wagę, ale ich dieta (może) wpływać na tempo wzrostu, więc Dietprzechwytywanie powinno być (mniej więcej) takie samo, ale nachylenie może być inne. Mogę uzyskać porównania par dla efektu przechwytywania w Dietnastępujący sposób:

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

i rzeczywiście nie różnią się one znacząco, ale jak mogę wykonać analogiczny test Time:Dietefektu? Samo wstawienie terminu interakcji mcppowoduje błąd:

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.
źródło
Ma Time*Diet, co jest tylko uproszczeniem Time + Diet + Time:Diet. Użycie anova(m)lub summary(m)potwierdza, że ​​termin interakcji znajduje się w modelu.
Dan M.

Odpowiedzi:

8

Domyślnie lmertraktuje poziom odniesienia predyktora jakościowego jako wartość bazową i szacuje parametry dla innych poziomów. Otrzymujesz więc parę porównań w domyślnym wyjściu i możesz uzyskać inne, używając releveldo zdefiniowania nowego poziomu odniesienia i ponownego dopasowania modelu. Ma to tę zaletę, że pozwala korzystać z porównań modeli lub MCMC w celu uzyskania wartości p, ale nie poprawia wielu porównań (chociaż można później zastosować własną korektę).

Aby użyć multcomp, musisz zdefiniować matrycę kontrastu. Każdy wiersz w macierzy kontrastu reprezentuje wagi efektów, które otrzymujesz w domyślnym wyjściu modelu, zaczynając od punktu przecięcia. Jeśli więc chcesz uzyskać efekt, który jest już zawarty w podstawowym wyjściu, po prostu umieść „1” w pozycji odpowiadającej temu efektowi. Ponieważ oszacowania parametrów odnoszą się do wspólnego poziomu odniesienia, można uzyskać porównania między dowolnymi dwoma innymi poziomami, ustawiając wagę jednego na „-1”, a drugiego „1”. Oto, jak to by działało dla Time:Dietwarunków w ChickWeightprzykładzie:

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))

Zastrzegający emptor: To podejście wydaje się wykorzystywać normalne przybliżenie, aby uzyskać wartości p, co jest nieco antykonserwatywne, a następnie stosuje pewną korektę dla wielu porównań. Rezultatem jest to, że ta metoda zapewnia łatwy dostęp do tylu oszacowań parametrów par i błędów standardowych, ile chcesz, ale wartości p mogą być lub nie być takie, jakie chcesz.

(Dzięki Scott Jackson z r-ling-lang-L za pomoc w tym)

Dan M.
źródło