Que faire avec une corrélation d'effets aléatoires égale à 1 ou -1?

L'occurrence pas si rare lorsqu'il s'agit de modèles mixtes maximaux complexes (estimation de tous les effets aléatoires possibles pour des données et un modèle donnés) est parfaite (+1 ou -1) ou une corrélation presque parfaite entre certains effets aléatoires. Aux fins de la discussion, observons le modèle et le résumé de modèle suivants

Model: Y ~ X*Cond + (X*Cond|subj)

# Y = logit variable  
# X = continuous variable  
# Condition = values A and B, dummy coded; the design is repeated 
#             so all participants go through both Conditions  
# subject = random effects for different subjects  

Random effects:
 Groups  Name             Variance Std.Dev. Corr             
 subject (Intercept)      0.85052  0.9222                    
         X                0.08427  0.2903   -1.00            
         CondB            0.54367  0.7373   -0.37  0.37      
         X:CondB          0.14812  0.3849    0.26 -0.26 -0.56
Number of obs: 39401, groups:  subject, 219

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)       2.49686    0.06909   36.14  < 2e-16 ***
X                -1.03854    0.03812  -27.24  < 2e-16 ***
CondB            -0.19707    0.06382   -3.09  0.00202 ** 
X:CondB           0.22809    0.05356    4.26 2.06e-05 ***

La raison supposée de ces corrélations parfaites est que nous avons créé un modèle trop complexe pour les données dont nous disposons. Le conseil commun qui est donné dans ces situations est (par exemple, Matuschek et al., 2017; papier ) de fixer les coefficients sur-paramétrisés à 0, car de tels modèles dégénérés ont tendance à réduire la puissance. Si nous observons un changement marqué des effets fixes dans un modèle réduit, nous devons l'accepter; s'il n'y a pas de changement, il n'y a aucun problème à accepter l'original.

Supposons cependant que nous ne nous intéressions pas seulement aux effets fixes contrôlés pour les RE (effets aléatoires), mais aussi à la structure des RE. Dans le cas donné, il serait théoriquement judicieux de supposer que Interceptla pente Xa une corrélation négative non nulle. Plusieurs questions suivent:

1. Que faire dans de telles situations? Doit-on rapporter la corrélation parfaite et dire que nos données ne sont pas "assez bonnes" pour estimer la "vraie" corrélation? Ou devrions-nous signaler le modèle de corrélation 0? Ou devrions-nous peut-être essayer de mettre une autre corrélation à 0 dans l'espoir que celle "importante" ne sera plus parfaite? Je ne pense pas qu'il y ait de réponses 100% correctes ici, je voudrais surtout entendre vos opinions.

2. Comment écrire le code qui fixe la corrélation de 2 effets aléatoires spécifiques à 0, sans influencer les corrélations entre les autres paramètres?

Matrices de covariance à effet aléatoire singulières

L'obtention d'une estimation de corrélation à effet aléatoire de +1 ou -1 signifie que l'algorithme d'optimisation a atteint "une limite": les corrélations ne peuvent pas être supérieures à +1 ni inférieures à -1. Même s'il n'y a pas d'erreurs ou d'avertissements de convergence explicites, cela indique potentiellement certains problèmes de convergence car nous ne nous attendons pas à ce que de vraies corrélations se situent à la frontière. Comme vous l'avez dit, cela signifie généralement qu'il n'y a pas suffisamment de données pour estimer tous les paramètres de manière fiable. Matuschek et al. 2017 dit que dans cette situation, le pouvoir peut être compromis.

Une autre façon d'atteindre une limite consiste à obtenir une estimation de la variance de 0: Pourquoi est-ce que j'obtiens une variance nulle d'un effet aléatoire dans mon modèle mixte, malgré certaines variations dans les données?

Les deux situations peuvent être vues comme l'obtention d'une matrice de covariance dégénérée d'effets aléatoires (dans votre exemple, la matrice de covariance de sortie est ); une variance nulle ou une corrélation parfaite signifie que la matrice de covariance n'est pas de rang complet et [au moins] une de ses valeurs propres est nulle. Cette observation suggère immédiatement qu'il existe d' autres façons plus complexes d'obtenir une matrice de covariance dégénérée: on peut avoir une matrice de covariance sans zéros ni corrélations parfaites mais néanmoins déficiente en rang (singulier). Bates et al. Modèles mixtes parcimonieux 20154 × $4\times 4$$4\times 4$(préimpression non publiée) recommande d'utiliser l'analyse en composantes principales (ACP) pour vérifier si la matrice de covariance obtenue est singulière. Si tel est le cas, ils suggèrent de traiter cette situation de la même manière que les situations singulières ci-dessus.

Alors que faire?

S'il n'y a pas suffisamment de données pour estimer tous les paramètres d'un modèle de manière fiable, alors nous devrions envisager de simplifier le modèle. En prenant votre exemple de modèle, X*Cond + (X*Cond|subj)il existe différentes façons de le simplifier:

1. Supprimez l'un des effets aléatoires, généralement la corrélation de premier ordre:

   X*Cond + (X+Cond|subj)

2. Débarrassez-vous de tous les paramètres de corrélation:

   X*Cond + (X*Cond||subj)

   Mise à jour: comme le note @Henrik, la ||syntaxe ne supprimera les corrélations que si toutes les variables à gauche sont numériques. Si des variables catégorielles (telles que Cond) sont impliquées, on devrait plutôt utiliser son afexpackage pratique (ou des solutions de contournement manuelles lourdes). Voir sa réponse pour plus de détails.

3. Débarrassez-vous de certains des paramètres de corrélation en divisant le terme en plusieurs, par exemple:

   X*Cond + (X+Cond|subj) + (0+X:Cond|subj)

4. Contraignez la matrice de covariance d'une manière spécifique, par exemple en définissant une corrélation spécifique (celle qui a atteint la limite) à zéro, comme vous le suggérez. Il n'y a aucun moyen intégré lme4pour y parvenir. Voir la réponse de @ BenBolker sur SO pour une démonstration de la façon d'y parvenir via un piratage intelligent.

Contrairement à ce que vous avez dit, je ne pense pas que Matuschek et al. 2017 recommande spécifiquement # 4. L'essentiel de Matuschek et al. 2017 et Bates et al. 2015 semble être que l'on part du modèle maximal à la Barr et al. 2013 , puis diminue la complexité jusqu'à ce que la matrice de covariance soit pleine. (De plus, ils recommandent souvent de réduire encore plus la complexité, afin d'augmenter la puissance.) Mise à jour: En revanche, Barr et al. recommander de réduire la complexité UNIQUEMENT si le modèle n'a pas convergé; ils sont prêts à tolérer des matrices de covariance singulières. Voir la réponse de @ Henrik.

Si l'on est d'accord avec Bates / Matuschek, alors je pense que c'est bien d'essayer différentes façons de diminuer la complexité afin de trouver celui qui fait le travail en faisant "le moins de dégâts". En regardant ma liste ci-dessus, la matrice de covariance d'origine a 10 paramètres; # 1 a 6 paramètres, # 2 a 4 paramètres, # 3 a 7 paramètres. Il est impossible de dire quel modèle supprimera les corrélations parfaites sans les adapter.

Mais que faire si ce paramètre vous intéresse?

La discussion ci-dessus traite la matrice de covariance à effet aléatoire comme un paramètre de nuisance. Vous vous posez une question intéressante de savoir quoi faire si vous êtes spécifiquement intéressé par un paramètre de corrélation que vous devez «abandonner» afin d'obtenir une solution complète et significative.

Notez que la fixation du paramètre de corrélation à zéro ne donnera pas nécessairement des BLUP ( ranef) non corrélés; en fait, ils pourraient même ne pas être affectés du tout (voir la réponse de @ Placidia pour une démonstration ). Une option consisterait donc à examiner les corrélations des BLUP et à en rendre compte.

Une autre option, peut-être moins attrayante, serait d'utiliser le traitement subjectcomme un effet fixe Y~X*cond*subj, d'obtenir les estimations pour chaque sujet et de calculer la corrélation entre eux. Cela équivaut à exécuter des Y~X*condrégressions distinctes pour chaque sujet séparément et à obtenir les estimations de corrélation à partir d'eux.

Voir également la section sur les modèles singuliers dans la FAQ sur les modèles mixtes de Ben Bolker:

Il est très fréquent que les modèles mixtes sur-ajustés entraînent des ajustements singuliers. Techniquement, la singularité signifie que certains des paramètres (décomposition de Cholesky de variance-covariance) correspondant aux éléments diagonaux du facteur de Cholesky sont exactement nuls, ce qui est le bord de l'espace réalisable, ou de manière équivalente que la matrice de variance-covariance a un zéro des valeurs propres (c.-à-d. est semi-définie positive plutôt que définie positive), ou (presque de manière équivalente) que certaines des variances sont estimées à zéro ou certaines des corrélations sont estimées à +/- 1.$\theta$

Je suis d'accord avec tout ce qui est dit dans la réponse d'amibe, qui fournit un excellent résumé de la discussion actuelle sur cette question. J'essaierai d'ajouter quelques points supplémentaires et, sinon, je ferai référence au polycopié de mon récent cours sur les modèles mixtes qui résume également ces points.

La suppression des paramètres de corrélation (options 2 et 3 dans la réponse de l'amibe) via ||ne fonctionne que pour les covariables numériques dans lmeret non pour les facteurs. Ceci est discuté en détail avec le code de Reinhold Kliegl .

Cependant, mon afexpackage fournit la fonctionnalité pour supprimer la corrélation également entre les facteurs if argument expand_re = TRUEdans l'appel à mixed()(voir aussi fonction lmer_alt()). Il le fait essentiellement en mettant en œuvre l'approche discutée par Reinhold Kliegl (c'est-à-dire en transformant les facteurs en covariables numériques et en spécifiant la structure des effets aléatoires sur ceux-ci).

Un exemple simple:

library("afex")
data("Machines", package = "MEMSS") # same data as in Kliegl code

# with correlation:
summary(lmer(score ~ Machine + (Machine  | Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr       
#  Worker   (Intercept) 16.6405  4.0793              
#           MachineB    34.5467  5.8776    0.48      
#           MachineC    13.6150  3.6899   -0.37  0.30
#  Residual              0.9246  0.9616              
# Number of obs: 54, groups:  Worker, 6

## crazy results:
summary(lmer(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name        Variance Std.Dev. Corr     
#  Worker   (Intercept)  0.2576  0.5076            
#  Worker.1 MachineA    16.3829  4.0476            
#           MachineB    74.1381  8.6103   0.80     
#           MachineC    19.0099  4.3600   0.62 0.77
#  Residual              0.9246  0.9616            
# Number of obs: 54, groups:  Worker, 6

## as expected:
summary(lmer_alt(score ~ Machine + (Machine  || Worker), data=Machines))
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  16.600   4.0743  
#  Worker.1 re1.MachineB 34.684   5.8894  
#  Worker.2 re1.MachineC 13.301   3.6471  
#  Residual               0.926   0.9623  
# Number of obs: 54, groups:  Worker, 6

Pour ceux qui ne le savent pas afex, la principale fonctionnalité des modèles mixtes est de fournir des valeurs de p pour les effets fixes, par exemple:

(m1 <- mixed(score ~ Machine + (Machine  || Worker), data=Machines, expand_re = TRUE))
# Mixed Model Anova Table (Type 3 tests, KR-method)
# 
# Model: score ~ Machine + (Machine || Worker)
# Data: Machines
#    Effect      df        F p.value
# 1 Machine 2, 5.98 20.96 **    .002
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1

summary(m1)  
# [...]
# Random effects:
#  Groups   Name         Variance Std.Dev.
#  Worker   (Intercept)  27.4947  5.2435  
#  Worker.1 re1.Machine1  6.6794  2.5845  
#  Worker.2 re1.Machine2 13.8015  3.7150  
#  Residual               0.9265  0.9626  
# Number of obs: 54, groups:  Worker, 6
# [...]

Dale Barr de Barr et al. (2013), l'article est plus prudent lorsqu'il recommande de réduire la structure à effets aléatoires que ce qui est présenté dans la réponse d'Amoeba. Dans un récent échange sur Twitter, il a écrit:

• "La réduction du modèle introduit un risque inconnu d'anticonservativité, et doit être effectuée avec prudence, le cas échéant." et
• "Ma principale préoccupation est que les gens comprennent les risques associés à la réduction du modèle et que la minimisation de ce risque nécessite une approche plus conservatrice que celle couramment adoptée (par exemple, chaque pente testée à 0,05)."

La prudence est donc de mise.

En tant que critique, je peux également vous expliquer pourquoi nous, les Bates et al. (2015), le document n'a pas été publié. Moi et les deux autres examinateurs (qui ont signé, mais qui resteront sans nom ici) ont critiqué l'approche PCA (semble sans principes et rien ne prouve qu'elle soit supérieure en termes de pouvoir). En outre, je pense que les trois ont critiqué le fait que le document ne se concentre pas sur la question de savoir comment spécifier la structure à effets aléatoires, mais essaie également d'introduire des GAMM. Ainsi, l'article de Bates et al (2015) s'est transformé en Matuschek et al. (2017) , qui aborde la question de la structure à effets aléatoires avec des simulations et Baayen et al. (2017) article présentant les GAMM.

Mon examen complet de Bates et al. brouillon peut être trouvé ici . IIRC, les autres examens avaient en quelque sorte des points principaux similaires.

Moi aussi, j'ai eu ce problème lors de l'utilisation de l'estimation du maximum de vraisemblance - seulement j'utilise l'algorithme Goldstein IGLS tel qu'implémenté via le logiciel MLwiN et non LME4 dans R. Cependant, dans chaque cas, le problème a été résolu lorsque je suis passé à l'estimation MCMC en utilisant le même Logiciel. J'ai même eu une corrélation supérieure à 3 qui s'est résolue lorsque j'ai changé d'estimation. En utilisant IGLS, la corrélation est calculée après estimation comme la covariance divisée par le produit de la racine carrée du produit des variances associées - et cela ne tient pas compte de l'incertitude dans chacune des estimations constituantes.

Le logiciel IGLS ne «sait» pas que la covariance implique une corrélation et calcule simplement les estimations d'une fonction de variance constante, linéaire, quadratique, etc. En revanche, l'approche MCMC est construite sur l'hypothèse d'échantillons d'une distribution normale multivariée qui correspond à des variances et des covariances avec de bonnes propriétés et une propagation d'erreur totale afin que l'incertitude dans l'estimation des covariances soit prise en compte dans l'estimation des variances et vice versa.

MLwin appartient à la chaîne d'estimation MCMC avec les estimations IGLS et la matrice de covariance à variance définie non négative pourrait devoir être modifiée en changeant la covariance à zéro au début avant de commencer l'échantillonnage.

Pour un exemple travaillé, voir

Développer des modèles à plusieurs niveaux pour analyser la contextualité, l'hétérogénéité et le changement à l'aide de MLwiN 3, Volume 1 (mis à jour en septembre 2017); Le volume 2 est également sur RGate

https://www.researchgate.net/publication/320197425_Vol1Training_manualRevisedSept2017

Annexe au chapitre 10