Lme4::glmer vs. command melogit de Stata

Question

Récemment, j'ai essayé d'ajuster de nombreux modèles d'effets aléatoires à des ensembles de données relativement importants. Disons environ 50 000 personnes (ou plus) observées jusqu'à 25 points temporels. Avec un échantillon aussi important, nous incluons de nombreux prédicteurs pour lesquels nous ajustons - peut-être une cinquantaine d'effets fixes. Je modélise le modèle sur un résultat binaire en utilisant lme4::glmer dans R, avec des intercepteurs aléatoires pour chaque sujet. Je ne peux pas entrer dans les détails sur les données, mais le format de base de la commande glmer que j'ai utilisée était :

fit <-  glmer(outcome ~ treatment + study_quarter + dd_quarter + (1|id),
              family = "binomial", data = dat)

où study_quarter et dd_quarter sont des facteurs avec environ 20 niveaux chacun.

Lorsque j'essaie d'ajuster ce modèle dans R, cela prend environ 12 à 15 heures et renvoie une erreur indiquant qu'il n'a pas convergé. J'ai fait beaucoup de dépannage (par exemple, en suivant ces directives), sans amélioration. Et la convergence n'est même pas proche à la fin (le gradient maximal est d'environ 5 à 10, alors que le critère de convergence est de 0,001 je pense).

J'ai alors essayé d'ajuster le modèle dans Stata, en utilisant la commande melogit. Le modèle a été ajusté en moins de 2 minutes, sans problèmes de convergence. La commande Stata correspondante est

melogit outcome treatment i.study_quarter i.dd_quarter || id:

Qu'est-ce qui se passe ? Est-ce que Stata dispose simplement d'un meilleur algorithme d'ajustement, ou mieux optimisé pour les modèles et les ensembles de données volumineux ? C'est vraiment surprenant à quel point les temps d'exécution étaient différents.

Answer 1

Le ajustement de glmer sera probablement beaucoup plus rapide avec l'argument optionnel nAGQ=0L. Vous avez de nombreux paramètres à effets fixes (20 niveaux pour chacun de study_quarter et dd_quarter générant un total de 28 contrastes) et la méthode d'optimisation par défaut (correspondant à nAGQ=1L) place tous ces coefficients dans l'appel d'optimisation non linéaire général. Avec nAGQ=0L, ces coefficients sont optimisés dans l'algorithme bien plus rapide de moindres carrés itérativement pondérés et pénalisés (PIRLS). La valeur par défaut fournit généralement une meilleure estimation dans le sens où le déviance à l'estimation est plus faible, mais la différence est généralement très faible et la différence de temps est énorme.

J'ai rédigé les différences entre ces algorithmes dans un Jupyter notebook nAGQ.ipynb. Cette explication utilise le package MixedModels pour Julia au lieu de lme4 mais les méthodes sont similaires. (Je suis l'un des auteurs de lme4 et l'auteur de MixedModels.)

Si vous avez l'intention de réaliser un grand nombre de ajustements GLMM, je vous conseillerais de le faire en Julia avec MixedModels. C'est souvent beaucoup plus rapide que R, même avec tout le code compliqué dans lme4.

Answer 2

Êtes-vous sûr que Stata lit l'intégralité du fichier?

http://www.stata.com/manuals13/rlimits.pdf

Je demande car il me semble que vous avez 50k personnes observées 25 fois (1,250,000 lignes); selon la version de Stata que vous utilisez, vous pourriez obtenir des résultats tronqués.

MODIFICATION Puisqu'il ne s'agit pas d'un problème de longueur de fichier, avez-vous essayé les autres packages pour les effets mixtes comme nlme? Je pense que le modèle d'effets mixtes non linéaires traiterait vos données un peu plus rapidement.

MODIFICATION Cette ressource peut être plus utile que tout ce qui concerne les différents modèles: https://stats.stackexchange.com/questions/173813/r-mixed-models-lme-lmer-or-both-which-one-is-relevant-for-my-data-and-why