Valeur de p pour le terme d'interaction dans les modèles à effets mixtes utilisant lme4

J'analyse certaines données comportementales en utilisant lme4in R, principalement en suivant les excellents tutoriels de Bodo Winter , mais je ne comprends pas si je gère correctement les interactions. Pire encore, personne d'autre impliqué dans cette recherche n'utilise des modèles mixtes, donc je suis un peu à la dérive quand il s'agit de s'assurer que les choses vont bien.

Plutôt que de simplement lancer un appel à l'aide, j'ai pensé que je devrais faire de mon mieux pour interpréter le problème, puis demander vos corrections collectives. Quelques autres apartés sont:

• En écrivant, j'ai trouvé cette question , montrant que nlmedonner plus directement des valeurs de p pour les termes d'interaction, mais je pense qu'il est toujours valable de se poser en relation avec lme4.
• Livius'La réponse à cette question a fourni des liens vers de nombreuses lectures supplémentaires, que j'essaierai de parcourir au cours des prochains jours, je vais donc commenter tout progrès qui en découlera.

Dans mes données, j'ai une variable dépendante dv, une conditionmanipulation (0 = contrôle, 1 = condition expérimentale, ce qui devrait entraîner une situation plus élevée dv), et également une condition préalable, étiquetée appropriate: les essais codés 1pour cela devraient montrer l'effet, mais les essais codés 0pourraient non, car il manque un facteur crucial.

J'ai également inclus deux interceptions aléatoires, pour subjectet pour target, reflétant des dvvaleurs corrélées dans chaque sujet et dans chacun des 14 problèmes résolus (chaque participant a résolu à la fois un contrôle et une version expérimentale de chaque problème).

library(lme4)
data = read.csv("data.csv")

null_model        = lmer(dv ~ (1 | subject) + (1 | target), data = data)
mainfx_model      = lmer(dv ~ condition + appropriate + (1 | subject) + (1 | target),
                         data = data)
interaction_model = lmer(dv ~ condition + appropriate + condition*appropriate +
                              (1 | subject) + (1 | target), data = data)
summary(interaction_model)

Production:

## Linear mixed model fit by REML ['lmerMod']
## ...excluded for brevity....
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  subject  (Intercept) 0.006594 0.0812  
##  target   (Intercept) 0.000557 0.0236  
##  Residual             0.210172 0.4584  
## Number of obs: 690, groups: subject, 38; target, 14
## 
## Fixed effects:
##                                Estimate Std. Error t value
## (Intercept)                    0.2518     0.0501    5.03
## conditioncontrol               0.0579     0.0588    0.98
## appropriate                   -0.0358     0.0595   -0.60
## conditioncontrol:appropriate  -0.1553     0.0740   -2.10
## 
## Correlation of Fixed Effects:
## ...excluded for brevity.

L'ANOVA montre alors interaction_modelun ajustement significativement meilleur que celui mainfx_model, à partir duquel je conclus qu'il y a une interaction significative présente (p = 0,035).

anova(mainfx_model, interaction_model)

Production:

## ...excluded for brevity....
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## mainfx_model       6 913 940   -450      901                          
## interaction_model  7 910 942   -448      896  4.44      1      0.035 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

À partir de là, j'isole un sous-ensemble des données pour lesquelles l' appropriateexigence est satisfaite (c.-à-d. appropriate = 1), Et pour qu'il s'adapte à un modèle nul et à un modèle incluant conditioncomme effet, compare à nouveau les deux modèles en utilisant l'ANOVA, et voilà, je trouve conditionest un prédicteur significatif.

good_data = data[data$appropriate == 1, ]
good_null_model   = lmer(dv ~ (1 | subject) + (1 | target), data = good_data)
good_mainfx_model = lmer(dv ~ condition + (1 | subject) + (1 | target), data = good_data)

anova(good_null_model, good_mainfx_model)

Production:

## Data: good_data
## models:
## good_null_model: dv ~ (1 | subject) + (1 | target)
## good_mainfx_model: dv ~ condition + (1 | subject) + (1 | target)
##                   Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
## good_null_model    4 491 507   -241      483                          
## good_mainfx_model  5 487 507   -238      477  5.55      1      0.018 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Je ne vois pas grand chose à dire ici. Je pense que vous avez fait du bon travail.

Il y a plusieurs façons dont les gens ont discuté pour tester les effets et obtenir des valeurs de p pour les modèles complexes d'effets mixtes. Il y a un bon aperçu ici . Le mieux est d'utiliser des méthodes de calcul intensif (bootstrap ou méthodes bayésiennes), mais c'est plus avancé pour la plupart des gens. La deuxième meilleure méthode (et la plus pratique) consiste à utiliser un test de rapport de vraisemblance. C'est ce que anova()(techniquement ? Anova.merMod () ) fait. Il est important d'utiliser uniquement un test de rapport de vraisemblance des modèles qui correspondaient à la probabilité maximale complète , plutôt que la probabilité maximale restreinte(REML). En revanche, pour votre modèle final et pour l'interprétation, vous souhaitez utiliser REML. C'est déroutant pour beaucoup de gens. Dans votre sortie, nous voyons que vous ajustez vos modèles avec REML (c'est parce que l'option est définie TRUEpar défaut dans lmer(). Cela signifierait que votre test n'est pas valide, cependant, parce que c'est une erreur courante, anova.merMod()contient un refitargument qui par la valeur par défaut est définie sur TRUE, et vous ne l'avez pas modifiée. Ainsi, la prévoyance des développeurs de packages vous y a sauvé.

En ce qui concerne votre stratégie pour déballer l'interaction, ce que vous avez fait est bien. N'oubliez pas que l'interaction utilise toutes les données pour son test. Il est possible d'avoir une interaction significative mais aucun des tests stratifiés ne soit significatif, ce qui déroute certaines personnes. (Il ne semble cependant pas vous être arrivé.)

Je suis moi-même un novice et je suis les conseils de Zuur et al .. J'utilise à lmepartir du nlmepackage au lieu de lme4quand j'ai besoin d'ajouter une structure d'erreur hiérarchique à un modèle par ailleurs linéaire. Ma réponse est peut-être loin.

Deux commentaires:

(1) Je ne suis pas sûr qu'il soit logique de tester conditiondans le sous-ensemble appropriate==1uniquement. Si vous souhaitez obtenir des valeurs de p pour les effets principaux, vous pouvez utiliser à Anovapartir du package 'car':

require(car)
Anova(M,type="III")# type III Sum of Squares. M was fitted using method=REML

Si vous souhaitez résoudre l'interaction, vous pouvez exécuter des comparaisons par paires directement (?) Ou faire ce que vous avez fait mais sur les deux sous-ensembles (c'est-à-dire également avec le sous-ensemble où appropriate==0).

(2) Vous voudrez peut-être d'abord sélectionner votre structure d'erreur au lieu de supposer que (1 | subject) + (1 | target)c'est la meilleure structure d'erreur. D'après ce que vous avez écrit, je suppose qu'il conditions'agit d'un facteur intra-sujet, alors qu'il appropriates'agit soit d'un facteur inter-sujet, soit d'un facteur inter-cible. Vous pouvez ajouter des pentes pour les facteurs intra-sujet et / ou intra-cible, par exemple: dv ~ condition + appropriate + (1+condition | subject) + (1 | target)ajoute une pente aléatoire pour le facteur intra-sujet condition. Aucune pente n'est nécessaire pour les facteurs entre sujets / cibles.

À votre santé