Comparaisons multiples du modèle mixte pour l'interaction entre le prédicteur continu et le prédicteur catégorique

11

Je voudrais utiliser lme4pour ajuster une régression à effets mixtes et multcomppour calculer les comparaisons par paire. J'ai un ensemble de données complexe avec plusieurs prédicteurs continus et catégoriques, mais ma question peut être démontrée en utilisant l' ChickWeightensemble de données intégré comme exemple:

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

Timeest continu et Dietcatégorique (4 niveaux) et il y a plusieurs poussins par régime. Tous les poussins ont commencé à peu près au même poids, mais leur régime alimentaire (peut) affecter leur taux de croissance, donc les Dietinterceptions devraient être (plus ou moins) les mêmes, mais les pentes peuvent être différentes. Je peux obtenir les comparaisons par paire pour l'effet d'interception Dietcomme ceci:

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

et, en effet, ils ne sont pas significativement différents, mais comment puis-je faire le test analogue pour l' Time:Dieteffet? La simple saisie du terme d'interaction mcpproduit une erreur:

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.
la source
Il a Time*Diet, ce qui est juste une simplification de Time + Diet + Time:Diet. L'utilisation de anova(m)ou summary(m)confirme que le terme d'interaction se trouve dans le modèle.
Dan M.

Réponses:

8

Par défaut, lmertraite le niveau de référence d'un prédicteur catégorique comme référence et estime les paramètres pour les autres niveaux. Ainsi, vous obtenez des comparaisons par paires dans la sortie par défaut et vous pouvez obtenir les autres en utilisant relevelpour définir un nouveau niveau de référence et réajuster le modèle. Cela a l'avantage de vous permettre d'utiliser des comparaisons de modèles ou MCMC pour obtenir les valeurs de p, mais ne corrige pas les comparaisons multiples (bien que vous puissiez appliquer votre propre correction par la suite).

Pour l'utiliser multcomp, vous devez définir la matrice de contraste. Chaque ligne de la matrice de contraste représente les pondérations des effets que vous obtenez dans la sortie du modèle par défaut, en commençant par l'interception. Donc, si vous voulez un effet déjà inclus dans la sortie de base, il vous suffit de mettre un "1" à la position correspondant à cet effet. Étant donné que les estimations des paramètres sont relatives à un niveau de référence commun, vous pouvez obtenir des comparaisons entre deux autres niveaux en définissant le poids de l'un sur "-1" et de l'autre "1". Voici comment cela fonctionnerait pour les Time:Diettermes de l' ChickWeightexemple:

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

Caveat emptor: Cette approche semble utiliser l'approximation normale pour obtenir les valeurs de p, ce qui est quelque peu anti-conservateur, puis applique une correction pour les comparaisons multiples. Le résultat est que cette méthode vous donne un accès facile à autant d'estimations de paramètres par paire et d'erreurs standard que vous le souhaitez, mais les valeurs de p peuvent ou non être ce que vous voulez.

(Merci à Scott Jackson de r-ling-lang-L pour son aide)

Dan M.
la source