Comment adapter un modèle mixte avec une variable de réponse entre 0 et 1?

15

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' glmerappel, sinon il ne converge pas du tout.)

amibe dit réintégrer Monica
la source
1
Que diriez-vous de transformer les probabilités en premier? Pouvez-vous obtenir quelque chose qui est plus proche de la distribution normale avec, disons, une transformation logit? Ou l'arcsin-sqrt? Ce serait ma préférence plutôt que d'utiliser glmer. Ou dans votre solution de piratage, vous pouvez également essayer d'ajouter un effet aléatoire pour chaque observation pour tenir compte de la sous-dispersion en raison de votre choix de poids.
Aaron a quitté Stack Overflow le
Merci. Oui, je peux enregistrer le DV puis utiliser le modèle mixte gaussien (lmer), mais c'est aussi une sorte de hack, et j'ai lu que ce n'était pas recommandé. Je vais essayer un effet aléatoire pour chaque observation! En ce moment, j'essaie un modèle bêta mixte; lme4 ne peut pas le gérer, mais glmmadmb le peut. Lorsque je cours 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.
amibe dit Reinstate Monica
Je ne sais pas pour glmer mais pour glm cela peut aider: stats.stackexchange.com/questions/164120/… :
1
@Aaron: J'ai essayé d'ajouter + (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.)
Amoeba dit Reinstate Monica
1
Et ils sont absolument certains de leur valeur, comme 0,9, ou ont-ils également une `` marge d'erreur ''? Pouvez-vous supposer que la confiance rapportée par différents sujets est tout aussi précise?

Réponses:

21

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):

  1. p=m/nnnN

    glm(p ~ a+b+c, myData, family="binomial", weights=n)
  2. pp01

    betareg(p ~ a+b+c, myData)
  3. Logit transforme la réponse et utilise une régression linéaire. Ce n'est généralement pas conseillé.

    lm(log(p/(1-p)) ~ a+b+c, myData)
  4. 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 glmcomme suit:

    glm(p ~ a+b+c, myData, family="quasibinomial")

Les quatre mêmes méthodes sont disponibles avec des effets aléatoires.

  1. En utilisant l' weightsargument ( un , deux ):

    glmer(p ~ a+b+c + (1|subject), myData, family="binomial", weights=n)

    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).

  2. Utilisation du modèle mixte bêta:

    glmmadmb(p ~ a+b+c + (1|subject), myData, family="beta")

    ou

    glmmTMB(p ~ a+b+c + (1|subject), myData, 
            family=list(family="beta",link="logit"))

    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.

  3. En utilisant la transformation logit de la réponse:

    lmer(log(p/(1-p)) ~ a+b+c + (1|subject), myData)
  4. 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:

    myData$rowid = as.factor(1:nrow(myData))
    glmer(p ~ a+b+c + (1|subject) + (1|rowid), myData, family="binomial",
          glmerControl(optimizer="bobyqa"))

    Pour une raison quelconque, cela ne fonctionne pas correctement car se glmer()plaint de non-entier pet donne des estimations non- sens. Une solution que j'ai trouvée est d'utiliser une fausse constante weights=ket de s'assurer qu'elle p*kest toujours entière. Cela nécessite un arrondi, pmais en le sélectionnant k, il ne devrait pas avoir beaucoup d'importance. Les résultats ne semblent pas dépendre de la valeur de k.

    k = 100
    glmer(round(p*k)/k ~ a+b+c + (1|subject) + (1|rowid), myData, 
          family="binomial", weights=rowid*0+k, glmerControl(optimizer="bobyqa"))

    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: glmmadmbprend cinq à dix minutes pour s'exécuter (et se plaint toujours qu'elle n'a pas convergé!), Alors qu'elle lmerfonctionne en une fraction de seconde et glmerprend quelques secondes. Mise à jour: j'ai essayé glmmTMBcomme suggéré dans les commentaires de @BenBolker et cela fonctionne presque aussi vite que possible glmer, 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.

amibe dit réintégrer Monica
la source
1
Belle réponse, bien expliquée et liée aux modèles sans effets aléatoires. Je n'appellerais pas la triche # 3 (la transformation), ces types de transformations sont très courants dans des analyses comme celles-ci. Je dirais plutôt que # 3 et # 4 font des hypothèses sur la relation entre la distribution des données, et donc aussi sur la relation entre la moyenne et la variance, et simplement parce que # 4 modélise à l'échelle que les données a été recueilli ne signifie pas que ces hypothèses vont être meilleures.
Aaron a quitté Stack Overflow le
1
# 3 suppose que le logit des probabilités est normal avec une variance constante, tandis que # 4 suppose que la variance est proportionnelle à p (1-p). D'après votre description de l'ajustement, ceux-ci semblent être assez similaires pour ne pas trop d'importance. Et # 3 est presque certainement plus standard (en fonction de votre public), donc si les diagnostics sont raisonnables, c'est celui que je préférerais.
Aaron a quitté Stack Overflow le
1
une autre possibilité est d'utiliser glmmTMB ; après l'installation avec devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB"), l'utilisation glmmTMB(p ~ a+b+c + (1|subject), myData, family=list(family="beta",link="logit"))devrait fonctionner ...
Ben Bolker
@BenBolker Merci! Y a-t-il une raison de préférer glmmTMB à glmmADMB (pour les modèles bêta) ou vice versa? L'un de ces packages est-il plus récent ou développé plus activement? En dehors de cela, puis-je demander quelle approche parmi celles énumérées dans cette réponse - glmm gaussien après transformation logit, glmm bêta ou glmm binomial avec le terme (1 | rowid) - trouvez-vous généralement préférable?
Amoeba dit Reinstate Monica
1
Je préfère la version bêta GLMM si possible - c'est le modèle statistique qui vise à mesurer les changements de proportions entre les covariables / groupes. glmmTMBest plus rapide et plus stable que glmmADMB, et sous (légèrement) plus actif, bien qu'il ne soit pas aussi mature.
Ben Bolker