Glmm binomial avec une variable catégorique avec plein succès

11

Je lance un glmm avec une variable de réponse binomiale et un prédicteur catégorique. L'effet aléatoire est donné par le plan imbriqué utilisé pour la collecte de données. Les données ressemblent à ceci:

m.gen1$treatment
 [1] sucrose      control      protein      control      no_injection .....
Levels: no_injection control sucrose protein
m.gen1$emergence 
 [1]  1  0  0  1  0  1  1  1  1  1  1  0  0....
> m.gen1$nest
 [1] 1  1  1  2  2  3  3  3  3  4  4  4  .....
Levels: 1 2 3 4 5 6 8 10 11 13 15 16 17 18 20 22 24

Le premier modèle que je lance ressemble à ceci

m.glmm.em.<-glmer(emergence~treatment + (1|nest),family=binomial,data=m.gen1)

Je reçois deux avertissements qui ressemblent à ceci:

Warning messages:
1: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model failed to converge with max|grad| = 0.0240654 (tol = 0.001, component 4)
2: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv,  :
  Model is nearly unidentifiable: large eigenvalue ratio
 - Rescale variables?

Le résumé du modèle montre que l'un des traitements a une erreur standard anormalement grande, que vous pouvez voir ici:

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)  
(Intercept)         2.565      1.038   2.472   0.0134 *
treatmentcontrol   -1.718      1.246  -1.378   0.1681  
treatmentsucrose   16.863   2048.000   0.008   0.9934  
treatmentprotein   -1.718      1.246  -1.378   0.1681 

J'ai essayé les différents optimiseurs du contrôle glmer et des fonctions d'autres packages, et j'obtiens une sortie similaire. J'ai exécuté le modèle en utilisant glm en ignorant l'effet aléatoire et le problème persiste. En explorant les données, j'ai réalisé que le traitement avec un Std élevé. l'erreur ne réussit que dans la variable de réponse. Juste pour vérifier si cela pourrait être à l'origine du problème, j'ai ajouté un faux point de données avec un "échec" pour ce traitement et le modèle fonctionne correctement et donne une erreur standard raisonnable. Vous pouvez le voir ici:

Fixed effects:
                 Estimate Std. Error z value Pr(>|z|)  
(Intercept)        3.4090     1.6712   2.040   0.0414 *
treatmentcontrol  -1.8405     1.4290  -1.288   0.1978  
treatmentsucrose  -0.2582     1.6263  -0.159   0.8738  
treatmentprotein  -2.6530     1.5904  -1.668   0.0953 .

Je me demandais si mon intuition était bonne quant au manque d'échecs pour ce traitement empêchant une bonne estimation, et comment puis-je contourner ce problème.

Merci d'avance!

amibe dit réintégrer Monica
la source

Réponses:

15

Votre intuition est exacte. Ce phénomène est appelé séparation complète . Vous pouvez trouver beaucoup (maintenant que vous connaissez son nom) Googler autour ... C'est assez discuté ici dans un contexte général , et ici dans le contexte des GLMM . La solution standard à ce problème consiste à ajouter un petit terme qui repousse les paramètres vers zéro - dans les contextes fréquentistes, cela s'appelle une méthode pénalisée ou à correction de biais . L'algorithme standard est dû à Firth (1993, "Réduction du biais des estimations du maximum de vraisemblance" Biometrika 80, 27-38), et est implémenté dans le package logistfsur CRAN. Dans les contextes bayésiens, cela est défini comme l'ajout d'un faible avant les paramètres à effet fixe.

À ma connaissance, l'algorithme de Firth n'a pas été étendu aux GLMM, mais vous pouvez utiliser l'astuce bayésienne en utilisant le paquet blme , qui place une fine couche bayésienne sur le dessus du lme4paquet. Voici un exemple de la discussion GLMM liée ci-dessus:

cmod_blme_L2 <- bglmer(predation~ttt+(1|block),data=newdat,
                   family=binomial,
                   fixef.prior = normal(cov = diag(9,4)))

Les deux premières lignes de cet exemple sont exactement les mêmes que celles que nous utiliserions dans le glmermodèle standard ; la dernière spécifie que l'a priori pour les effets fixes est une distribution normale multivariée avec une matrice de variance-covariance diagonale. La matrice est 4x4 (car nous avons 4 paramètres à effet fixe dans cet exemple), et la variance antérieure de chaque paramètre est 9 (correspondant à un écart type de 3, ce qui est assez faible - cela signifie que +/- 2SD est ( -6,6), qui est une très large plage sur l'échelle logit).

Les très grandes erreurs standard des paramètres dans votre exemple sont un exemple d'un phénomène étroitement lié à la séparation complète (il se produit chaque fois que nous obtenons des valeurs de paramètres extrêmes dans un modèle logistique) appelé effet Hauck-Donner .

Deux autres références potentiellement utiles (je n'y ai pas encore creusé moi-même):

  • Gelman A, Jakulin A, Pittau MG et Su TS (2008) Une distribution a priori par défaut faiblement informative pour les modèles de régression logistique et autres. Annals of Applied Statistics , 2, 1360–1383.
  • José Cortiñas Abrahantes et Marc Aerts (2012) Une solution à la séparation pour les données binaires en cluster Modélisation statistique 12 (1): 3–27 doi: 10.1177 / 1471082X1001200102

Une recherche Google plus récente sur "bglmer" séparation complète "" révèle:

  • Quiñones, AE et WT Wcislo. "Cryptic Extended Brood Care in the Facultively Eusocial Sweat Bee Megalopta genalis ." Insectes Sociaux 62.3 (2015): 307–313.
Ben Bolker
la source
wow merci beaucoup !! Cela est parfaitement logique et le modèle fonctionne désormais correctement avec le bglmer. J'aurais juste une autre question, puis-je utiliser les méthodes comme dans lme4 pour évaluer les effets aléatoires et fixes, en d'autres termes pour comparer différents modèles?
2
Je dirais que oui, mais je ne sais pas s'il existe un soutien formel et / ou évalué par les pairs à mon avis ...
Ben Bolker
Merci! C'est exactement mon problème aussi. Un suivi rapide: contrairement à votre exemple, qui a un facteur à 4 niveaux, j'ai un design 2 x 2 où chaque facteur a 2 niveaux (donc le total est toujours à 4 niveaux). Puis-je également utiliser diag (9,4) pour mon modèle? Je ne connais pas bien les matrices, donc je voulais vérifier. De même, pour justifier cette solution dans mon article, dois-je citer Firth (1993) ou existe-t-il un article plus directement pertinent, qui a mis en œuvre votre solution en utilisant bglmer ()?
Sol
2
voir la réponse mise à jour.
Ben Bolker
2
Je pense que oui - le nombre total de paramètres d'effets fixes ne devrait avoir qu'une importance.
Ben Bolker