Qu'est-ce qu'une carence en rang et comment y remédier?

Une régression logistique à l'aide de lme4 se termine par

Error in mer_finalize(ans) : Downdated X'X is not positive definite. 

Une cause probable de cette erreur est apparemment une insuffisance de rang. Qu'est-ce qu'une déficience de rang et comment devrais-je y remédier?

Le manque de rang dans ce contexte indique que les informations contenues dans vos données sont insuffisantes pour estimer le modèle souhaité. Il provient de nombreuses origines. Je parlerai ici de la modélisation dans un contexte assez général, plutôt que d’une régression logistique explicite, mais tout s’applique toujours au contexte spécifique.

La carence peut provenir simplement de trop peu de données. En général, vous ne pouvez pas estimer de manière unique n paramètres avec moins de n points de données. Cela ne signifie pas que tout ce dont vous avez besoin est de n points, car s'il y avait du bruit dans le processus, vous obtiendriez des résultats plutôt médiocres. Vous avez besoin de plus de données pour aider l’algorithme à choisir une solution qui représentera l’ensemble des données, avec un minimum d’erreur. C'est pourquoi nous utilisons les outils des moindres carrés. De combien de données avez-vous besoin? On m'a toujours posé cette question dans une vie antérieure, et la réponse a été supérieure à ce que vous avez ou à ce que vous pouvez obtenir. :)

Parfois, vous pouvez avoir plus de données que nécessaire, mais certains (trop) de points sont des réplicats. La réplication est BON dans le sens où elle contribue à réduire le bruit, mais pas à augmenter le rang numérique. Ainsi, supposons que vous n’ayez que deux points de données. Vous ne pouvez pas estimer un modèle quadratique unique à travers les points. Un million de répétitions de chaque point ne vous permettra toujours pas de passer plus d'une ligne droite, à travers ce qui n'est encore qu'une paire de points. Essentiellement, la réplication n'ajoute pas de contenu d'informations. Tout ce que cela fait, c'est réduire le bruit aux endroits où vous avez déjà des informations.

Parfois, vous avez des informations aux mauvais endroits. Par exemple, vous ne pouvez pas adapter un modèle quadratique à deux dimensions si vous ne disposez que de points situés tous en ligne droite dans deux dimensions. En d'autres termes, supposons que vous ayez des points dispersés uniquement le long de la ligne x = y dans le plan et que vous souhaitiez ajuster un modèle à la surface z (x, y). Même avec des zillions de points (pas même des répliques), vous aurez suffisamment d'informations pour estimer intelligemment plus qu'un modèle constant. Étonnamment, c'est un problème commun que j'ai vu dans les données échantillonnées. L'utilisateur se demande pourquoi ils ne peuvent pas construire un bon modèle. Le problème est intégré aux données mêmes qu'ils ont échantillonnées.

Parfois, il s'agit simplement d'un choix de modèle. Cela peut être considéré comme "pas assez de données", mais de l'autre côté. Vous souhaitez estimer un modèle compliqué, mais n'avez pas fourni suffisamment de données pour le faire.

Dans tous les cas ci-dessus, la solution consiste à obtenir davantage de données, échantillonnées intelligemment depuis des lieux fournissant des informations sur le processus qui vous manque actuellement. La conception d'expériences est un bon point de départ.

Cependant, même de bonnes données sont parfois insuffisantes, du moins numériquement. (Pourquoi de mauvaises données arrivent-elles à de bonnes données?) Le problème ici peut être lié au modèle. Cela ne réside peut-être que dans un mauvais choix d’unités. Cela peut provenir de la programmation informatique effectuée pour résoudre le problème. (Ugh! Par où commencer?)

Tout d'abord, parlons des unités et de la mise à l'échelle. Supposons que j'essaie de résoudre un problème où une variable est BEAUCOUP plus grand qu'un ordre de grandeur. Par exemple, supposons que j'ai un problème qui implique ma taille et la taille de mes chaussures. Je vais mesurer ma taille en nanomètres. Donc, ma taille serait d'environ 1,78 milliard (1,78e9) nanomètres. Bien sûr, je choisirai de mesurer ma pointure en kilo-parsecs, soit 9,14 e-21 kilo-parsecs. Lorsque vous modélisez par régression, la régression linéaire concerne l’algèbre linéaire, qui implique des combinaisons linéaires de variables. Le problème, c’est que ces nombres diffèrent d’énormément d’ordre de grandeur (et même pas des mêmes unités.) nombre,

L'astuce consiste généralement à utiliser des unités communes, mais pour certains problèmes, même lorsque les variables varient d'un ordre de grandeur à l'autre. Le plus important est d’adapter vos chiffres à une magnitude similaire.

Ensuite, vous pouvez voir des problèmes avec les grands nombres et une petite variation dans ces nombres. Ainsi, supposons que vous essayiez de construire un modèle polynomial d'ordre relativement élevé avec des données où toutes vos entrées se situent dans l'intervalle [1,2]. La mise au carré, la mise en cubes, etc., des nombres de l'ordre de 1 ou 2 ne poseront aucun problème lorsque vous travaillerez en arithmétique en double précision. Sinon, ajoutez 1e12 à chaque numéro. En théorie, les mathématiques le permettront. Tout ce qu'il fait, c'est déplacer tout modèle polynomial construit sur l'axe des x. Il aurait exactement la même forme, mais se traduirait par 1e12 à droite. En pratique, l'algèbre linéaire échouera lamentablement à cause de problèmes d'insuffisance de rang. Vous n'avez rien fait d'autre que traduire les données, mais vous commencez tout à coup à voir apparaître des matrices singulières.

En général, le commentaire consiste à suggérer de "centrer et d’échelonner vos données". En pratique, cela signifie déplacer et mettre à l'échelle les données de manière à ce qu'elles aient une moyenne proche de zéro et un écart-type d'environ 1, ce qui améliorera considérablement le conditionnement de la plupart des modèles polynomiaux, réduisant ainsi le nombre de problèmes d'insuffisance de rang.

Il existe d’autres raisons pour la déficience de rang. Dans certains cas, il est construit directement dans le modèle. Par exemple, supposons que je fournisse la dérivée d'une fonction, puis-je déduire de manière unique la fonction elle-même? Bien sûr que non, car l’intégration implique une constante d’intégration, paramètre inconnu généralement déduit de la connaissance de la valeur de la fonction à un moment donné. En fait, cela se produit parfois aussi dans les problèmes d’estimation, où la singularité d’un système découle de la nature fondamentale du système à l’étude.

J’ai sûrement omis quelques-unes des nombreuses raisons de l’insuffisance de rang dans un système linéaire, et cela fait trop longtemps que je discute. J'espère que j'ai réussi à expliquer ceux que j'ai couverts en termes simples et un moyen d'atténuer le problème.

Pour définir le rang d’une matrice, vous pouvez vous référer à n’importe quel bon manuel d’algèbre linéaire ou consulter la page Wikipedia .

Une matrice est dite de rang plein si , et ses colonnes ne sont pas une combinaison linéaire de l'autre. Dans ce cas, la matrice est définie positive, ce qui implique qu'elle a un inverse .$n \times p$$X$$n \geq p$$p \times p$$X^TX$$(X^TX)^{-1}$

Si n’est pas un rang complet, l’une des colonnes est expliquée en détail par les autres, en ce sens qu’elle est une combinaison linéaire des autres. Un exemple trivial est quand une colonne est dupliquée. Cela peut également arriver si vous avez une variable 0-1 et qu'une colonne ne contient que 0 ou seulement 1. Dans ce cas, le rang de la matrice est inférieur à et n'a pas d'inverse.$X$$X$$n$$X^TX$

Puisque la solution de nombreux problèmes de régression (y compris la régression logistique) implique le calcul intermédiaire , il est alors impossible d'estimer les paramètres du modèle. Par curiosité, vous pouvez vérifier ici comment ce terme est impliqué dans la formule de la régression linéaire multiple.$(X^TX)^{-1}$

C'était tout pour le déficit de grade absolu. Mais parfois, le problème apparaît lorsque la matrice n’est "presque" pas classée dans son intégralité, comme détaillé par @woodchips. Ce problème est généralement appelé multicolinéarité . Cette question est assez courante, vous pouvez en savoir plus sur la façon de la traiter sur des posts liés ici et là .$X$

La réponse de user974 est fantastique du point de vue de la modélisation, tandis que celle de gui11aume est fantastique du point de vue mathématique. Je veux affiner l'ancienne réponse strictement du point de vue de la modélisation mixte: en particulier une généralisation perspective de modélisation mixte (GLMM). Comme vous pouvez le constater, vous avez référencé la fonction R mer_finalizequi se trouve dans le lme4package fantastique . Vous dites également que vous adaptez un modèle de régression logistique.

De nombreux problèmes surgissent avec de tels types d'algorithmes numériques. La question de la structure matricielle du modèle matrice d’effets fixes mérite certainement d’être examinée, comme l’a déjà évoqué l’utilisateur974. Mais ce qui est très facile à évaluer, simplement calculer la model.matrixde vos formula=et data=arguments dans un modèle et prendre son déterminant en utilisant la detfonction. Cependant, les effets aléatoires compliquent grandement l'interprétation, la routine d'estimation numérique et l'inférence sur les effets fixes (ce que vous appelez généralement des coefficients de régression dans un modèle de régression "normal").

Supposons que dans le cas le plus simple, vous n’êtes équipé que d’un modèle d’interception aléatoire. Ensuite, vous considérez fondamentalement qu'il existe des milliers de sources d'hétérogénéité non mesurées qui sont maintenues constantes dans des mesures répétées au sein de grappes. Vous estimez une "grande" interception, mais vous tenez compte de l'hétérogénéité en supposant que les intercepts spécifiques à un cluster ont une distribution moyenne moyenne. Les intercepts sont estimés de manière itérative et utilisés pour mettre à jour les effets du modèle jusqu'à ce que la convergence soit atteinte (la probabilité de log (ou une approximation de celle-ci) est maximisée). Le modèle mixte est très facile à envisager, mais mathématiquement, la probabilité est très complexe et présente des problèmes avec des singularités, des minimae locaux et des points limites (rapports de cotes = 0 ou infini). Les modèles mixtes n'ont pas de vraisemblance quadratique comme les GLM canoniques.

Malheureusement, Venerables et Ripley n’ont pas beaucoup investi dans le diagnostic des échecs de convergence comme le vôtre. Il est pratiquement impossible même de spéculer sur la myriade d’erreurs possibles menant à un tel message. Considérez alors les types de diagnostics que j'utilise ci-dessous:

1. Combien d'observations y a-t-il par groupe?
2. Quels sont les résultats d'un ajustement de modèle marginal utilisant GEE?
3. Quel est le CIC des clusters? L'hétérogénéité au sein du groupe est-elle proche de l'hétérogénéité entre les groupes?
4. Ajustez un estimateur en une étape et examinez les effets aléatoires estimés. Sont-ils à peu près normaux?
5. Ajuster le modèle mixte bayésien et regarder la distribution postérieure pour les effets fixes. Est-ce qu'ils semblent avoir une distribution approximativement normale?
6. Examinez un graphique des grappes montrant l’exposition ou le régresseur d’intérêt par rapport au résultat à l’aide d’un lisseur. Les tendances sont-elles cohérentes et claires ou existe-t-il de nombreuses manières possibles de les expliquer? (Par exemple, quel est le "risque" parmi les sujets non exposés, l'exposition semble-t-elle protectrice ou nuisible?)
7. Est-il possible de restreindre l'échantillon aux sujets ayant un nombre d'observations par personne suffisant (par exemple, n = 5 ou n = 10) pour estimer l'effet "idéal" du traitement?

Alternativement, vous pouvez envisager différentes approches de modélisation:

1. Y a-t-il assez peu de clusters ou de points de temps pour pouvoir utiliser un effet fixe (comme des indicateurs de groupe ou un effet de temps polynomial) pour modéliser l'hétérogénéité au niveau du cluster / autorégressive?
2. Un modèle marginal est-il approprié (utiliser un GEE pour améliorer l'efficacité de l'estimation de l'erreur standard, mais en utilisant uniquement des effets fixes)
3. Un modèle bayésien avec un préalable informatif sur les effets aléatoires pourrait-il améliorer l'estimation?