Mise à jour de mai 2017 : En fin de compte, un lof de ce que j'ai écrit ici est un peu malicieux . Certaines mises à jour sont effectuées tout au long de la publication.
Je suis d'accord avec ce qui a déjà été dit par Ben Bolker (merci pour les remerciements à afex::mixed()
) mais permettez-moi d'ajouter quelques réflexions plus générales et spécifiques sur cette question.
Concentrez-vous sur les effets fixes par rapport aux effets aléatoires et sur la façon de communiquer les résultats
Pour le type de recherche expérimentale représenté dans l'exemple de jeu de données de Jonathan Baron, la question importante est généralement de savoir si un facteur manipulé a ou non un effet global. Par exemple, trouvons-nous un effet principal global ou une interaction de Task
? Un point important est que dans ces ensembles de données, tous les facteurs sont généralement sous contrôle expérimental complet et assignés au hasard. Par conséquent, l'intérêt se porte généralement sur les effets fixes.
En revanche, les composantes des effets aléatoires peuvent être considérées comme des paramètres de «nuisance» qui capturent la variance systématique (c.-à-d. Les différences interindividuelles dans la taille de l'effet) qui ne sont pas nécessairement importantes pour la question principale. De ce point de vue, la suggestion d'utiliser la structure à effets aléatoires maximales préconisée par Barr et al. suit un peu naturellement. Il est facile d'imaginer qu'une manipulation expérimentale n'affecte pas tous les individus de la même manière et nous voulons contrôler cela. D'un autre côté, le nombre de facteurs ou de niveaux n'est généralement pas trop important, de sorte que le risque de surapprentissage semble relativement faible.
Par conséquent, je suivrais la suggestion de Barr et al. et spécifier une structure d'effets aléatoires maximale et rapporter les tests des effets fixes comme mes principaux résultats. Pour tester les effets fixes, je suggérerais également d'utiliser afex::mixed()
car il rapporte des tests d'effets ou de facteurs (au lieu de tester des paramètres) et calcule ces tests d'une manière quelque peu raisonnable (par exemple, utilise la même structure d'effets aléatoires pour tous les modèles dans lesquels un l'effet unique est supprimé, utilise des contrastes de somme à zéro, propose différentes méthodes pour calculer les valeurs de p , ...).
Qu'en est-il des données d'exemple
Le problème avec les données d'exemple que vous avez fournies est que pour cet ensemble de données, la structure d'effets aléatoires maximale conduit à un modèle sursaturé car il n'y a qu'un seul point de données par cellule du plan:
> with(df, table(Valence, Subject, Task))
, , Task = Cued
Subject
Valence Faye Jason Jim Ron Victor
Neg 1 1 1 1 1
Neu 1 1 1 1 1
Pos 1 1 1 1 1
, , Task = Free
Subject
Valence Faye Jason Jim Ron Victor
Neg 1 1 1 1 1
Neu 1 1 1 1 1
Pos 1 1 1 1 1
Par conséquent, lmer
s'étouffe sur la structure maximale des effets aléatoires:
> lmer(Recall~Task*Valence + (Valence*Task|Subject), df)
Error: number of observations (=30) <= number of random effects (=30) for term
(Valence * Task | Subject); the random-effects parameters and the residual variance
(or scale parameter) are probably unidentifiable
Malheureusement, il n'y a à ma connaissance aucun moyen convenu de résoudre ce problème. Mais laissez-moi en esquisser et en discuter:
Une première solution pourrait être de supprimer la pente aléatoire la plus élevée et de tester les effets de ce modèle:
require(afex)
mixed(Recall~Task*Valence + (Valence+Task|Subject), df)
Effect F ndf ddf F.scaling p.value
1 Task 6.56 1 4.00 1.00 .06
2 Valence 0.80 2 3.00 0.75 .53
3 Task:Valence 0.42 2 8.00 1.00 .67
Cependant, cette solution est un peu ad-hoc et pas trop motivée.
Mise à jour de mai 2017: c'est l'approche que j'approuve actuellement. Voir cet article de blog et le brouillon du chapitre dont je suis co-auteur , section "Structures d'effets aléatoires pour les conceptions ANOVA traditionnelles".
Une solution alternative (et qui pourrait être considérée comme préconisée par la discussion de Barr et al.) Pourrait être de toujours supprimer les pentes aléatoires pour le plus petit effet. Cela pose cependant deux problèmes: (1) Quelle structure d'effets aléatoires utilisons-nous pour déterminer quel est le plus petit effet et (2) R hésite à supprimer un effet d'ordre inférieur tel qu'un effet principal si des effets d'ordre supérieur tels qu'un l'interaction de cet effet est présente (voir ici ). En conséquence, il faudrait configurer cette structure d'effets aléatoires à la main et passer la matrice de modèle ainsi construite à l'appel lmer.
Une troisième solution pourrait être d'utiliser une paramétrisation alternative de la partie effets aléatoires, à savoir celle qui correspond au modèle RM-ANOVA pour ces données. Malheureusement (?), lmer
Ne permet pas les "variances négatives", donc ce paramétrage ne correspond pas exactement au RM-ANOVA pour tous les ensembles de données , voir la discussion ici et ailleurs (par exemple ici et ici ). Le "lmer-ANOVA" pour ces données serait:
> mixed(Recall~Task*Valence + (1|Subject) + (1|Task:Subject) + (1|Valence:Subject), df)
Effect F ndf ddf F.scaling p.value
1 Task 7.35 1 4.00 1.00 .05
2 Valence 1.46 2 8.00 1.00 .29
3 Task:Valence 0.29 2 8.00 1.00 .76
Compte tenu de tous ces problèmes, je n'utiliserais tout simplement pas lmer
pour ajuster des ensembles de données pour lesquels il n'y a qu'un seul point de données par cellule de la conception, à moins qu'une solution plus convenue pour le problème de la structure d'effets aléatoires maximale ne soit disponible.
Au lieu de cela, je voudrais que l' on puisse également utiliser l'ANOVA classique. Utiliser l'un des wrappers car::Anova()
dans afex
les résultats serait:
> aov4(Recall~Task*Valence + (Valence*Task|Subject), df)
Effect df MSE F ges p
1 Valence 1.44, 5.75 4.67 1.46 .02 .29
2 Task 1, 4 4.08 7.35 + .07 .05
3 Valence:Task 1.63, 6.52 2.96 0.29 .003 .71
A noter que afex
permet désormais aussi de renvoyer le modèle équipé avec aov
lequel on peut passer pour lsmeans
des tests post-hoc (mais pour le test des effets ceux rapportés par car::Anova
sont encore plus raisonnables):
> require(lsmeans)
> m <- aov4(Recall~Task*Valence + (Valence*Task|Subject), df, return = "aov")
> lsmeans(m, ~Task+Valence)
Task Valence lsmean SE df lower.CL upper.CL
Cued Neg 11.8 1.852026 5.52 7.17157 16.42843
Free Neg 10.2 1.852026 5.52 5.57157 14.82843
Cued Neu 13.0 1.852026 5.52 8.37157 17.62843
Free Neu 11.2 1.852026 5.52 6.57157 15.82843
Cued Pos 13.6 1.852026 5.52 8.97157 18.22843
Free Pos 11.0 1.852026 5.52 6.37157 15.62843
Confidence level used: 0.95
lmer
ne permet pas des variances négatives mais bien évidemment des corrélations négatives entre les composantes de variance.