J'essaie d'utiliser lme4::glmer()
pour adapter un modèle mixte généralisé binomial (GLMM) avec une variable dépendante qui n'est pas binaire, mais une variable continue entre zéro et un. On peut considérer cette variable comme une probabilité; en fait , il est probable que rapporté par des sujets humains (dans une expérience que je l' aide à analyser). C'est-à-dire que ce n'est pas une fraction "discrète", mais une variable continue.
Mon glmer()
appel ne fonctionne pas comme prévu (voir ci-dessous). Pourquoi? Que puis-je faire?
Modification ultérieure: ma réponse ci-dessous est plus générale que la version originale de cette question, j'ai donc modifié la question pour qu'elle soit plus générale également.
Plus de détails
Apparemment, il est possible d'utiliser la régression logistique non seulement pour le DV binaire mais aussi pour le DV continu entre zéro et un. En effet, quand je cours
glm(reportedProbability ~ a + b + c, myData, family="binomial")
Je reçois un message d'avertissement
Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!
mais un ajustement très raisonnable (tous les facteurs sont catégoriques, donc je peux facilement vérifier si les prédictions du modèle sont proches des moyennes inter-sujets, et elles le sont).
Cependant, ce que je veux réellement utiliser est
glmer(reportedProbability ~ a + b + c + (1 | subject), myData, family="binomial")
Il me donne l'avertissement identique, renvoie un modèle, mais ce modèle est clairement très hors tension; les estimations des effets fixes sont très éloignées de glm()
celles et des moyennes transversales. (Et je dois inclure glmerControl(optimizer="bobyqa")
dans l' glmer
appel, sinon il ne converge pas du tout.)
glmmadmb(reportedProbability ~ a + b + c + (1 | subject), myData, family="beta")
, j'obtiens un ajustement correct et des intervalles de confiance raisonnables, mais un avertissement de convergence a échoué : - / J'essaie de comprendre comment augmenter le nombre d'itérations. La bêta pourrait fonctionner pour moi car je n'ai pas de cas DV = 0 ou DV = 1.+ (1 | rowid)
à mon appel glmer et cela donne des estimations stables et des intervalles de confiance stables, indépendamment de mon choix de poids (j'ai essayé 100 et 500). J'ai également essayé d'exécuter lmer sur logit (rapportéProbabilité) et j'obtiens presque exactement la même chose. Les deux solutions semblent donc bien fonctionner! Beta MM avec glmmadmb donne également des résultats très proches, mais pour une raison quelconque, il ne parvient pas à converger complètement et prend une éternité à s'exécuter. Pensez à publier une réponse répertoriant ces options et expliquant un peu les différences et les avantages / inconvénients! (Les intervalles de confiance que je mentionne sont tous Wald.)Réponses:
Il est logique de commencer par un cas plus simple sans effets aléatoires.
Il existe quatre façons de gérer la variable de réponse continue de zéro à un qui se comporte comme une fraction ou une probabilité ( c'est notre fil le plus canonique / le plus voté / le plus consulté sur ce sujet, mais malheureusement, les quatre options ne sont pas discutées ici):
n
Logit transforme la réponse et utilise une régression linéaire. Ce n'est généralement pas conseillé.
Ajuster un modèle binomial, puis calculer les erreurs standard en tenant compte de la sur-dispersion. Les erreurs standard peuvent être calculées de différentes manières:
(a) échelles standard échelonnées via l'estimation de surdispersion ( un , deux ). C'est ce qu'on appelle le GLM "quasi-binomial".
(b) des erreurs types robustes via l'estimateur sandwich ( un , deux , trois , quatre ). C'est ce qu'on appelle le "logit fractionnaire" en économétrie.
Les (a) et (b) ne sont pas identiques (voir ce commentaire , et les sections 3.4.1 et 3.4.2 dans ce livre , et ce message SO et aussi celui-ci et celui-ci ), mais ont tendance à donner des résultats similaires. L'option (a) est mise en œuvre
glm
comme suit:Les quatre mêmes méthodes sont disponibles avec des effets aléatoires.
En utilisant l'
weights
argument ( un , deux ):Selon le deuxième lien ci-dessus, ce pourrait être une bonne idée de modéliser la surdispersion, voir là (et aussi # 4 ci-dessous).
Utilisation du modèle mixte bêta:
ou
S'il y a des zéros ou des uns exacts dans les données de réponse, alors on peut utiliser le modèle bêta zéro / un gonflé dans
glmmTMB
.En utilisant la transformation logit de la réponse:
Prise en compte de la surdispersion dans le modèle binomial. Cela utilise une astuce différente: ajouter un effet aléatoire pour chaque point de données:
Pour une raison quelconque, cela ne fonctionne pas correctement car se
glmer()
plaint de non-entierp
et donne des estimations non- sens. Une solution que j'ai trouvée est d'utiliser une fausse constanteweights=k
et de s'assurer qu'ellep*k
est toujours entière. Cela nécessite un arrondi,p
mais en le sélectionnantk
, il ne devrait pas avoir beaucoup d'importance. Les résultats ne semblent pas dépendre de la valeur dek
.Mise à jour ultérieure (janvier 2018): il s'agit peut-être d'une approche non valide. Voir la discussion ici . Je dois enquêter davantage sur cela.
Dans mon cas spécifique, l'option # 1 n'est pas disponible.
L'option n ° 2 est très lente et a des problèmes de convergence:Mise à jour: j'ai essayéglmmadmb
prend cinq à dix minutes pour s'exécuter (et se plaint toujours qu'elle n'a pas convergé!), Alors qu'ellelmer
fonctionne en une fraction de seconde etglmer
prend quelques secondes.glmmTMB
comme suggéré dans les commentaires de @BenBolker et cela fonctionne presque aussi vite que possibleglmer
, sans aucun problème de convergence. Voilà donc ce que je vais utiliser.Les options 3 et 4 donnent des estimations très similaires et des intervalles de confiance de Wald très similaires (obtenus avec
confint
). Je ne suis pas un grand fan de # 3, car c'est une sorte de tricherie. Et # 4 se sent un peu hacky.Un grand merci à @Aaron qui m'a pointé vers # 3 et # 4 dans son commentaire.
la source
devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB")
, l'utilisationglmmTMB(p ~ a+b+c + (1|subject), myData, family=list(family="beta",link="logit"))
devrait fonctionner ...glmmTMB
est plus rapide et plus stable queglmmADMB
, et sous (légèrement) plus actif, bien qu'il ne soit pas aussi mature.