Wie passt man ein gemischtes Modell mit einer Antwortvariablen zwischen 0 und 1 an?

15

Ich versuche, lme4::glmer()ein binomiales verallgemeinertes gemischtes Modell (GLMM) mit abhängiger Variable anzupassen, die nicht binär ist, sondern eine kontinuierliche Variable zwischen Null und Eins. Man kann sich diese Variable als Wahrscheinlichkeit vorstellen; Tatsächlich ist es die Wahrscheinlichkeit, die von menschlichen Probanden angegeben wurde (in einem Experiment, das ich bei der Analyse unterstütze). Dh es ist kein "diskreter" Bruch, sondern eine stetige Variable.

Mein glmer()Anruf funktioniert nicht wie erwartet (siehe unten). Warum? Was kann ich tun?

Später bearbeiten: Meine Antwort unten ist allgemeiner als die Originalversion dieser Frage, daher habe ich die Frage geändert, um sie auch allgemeiner zu gestalten.


Mehr Details

Anscheinend ist es möglich, die logistische Regression nicht nur für binäres DV, sondern auch für kontinuierliches DV zwischen null und eins zu verwenden. In der Tat, wenn ich renne

glm(reportedProbability ~ a + b + c, myData, family="binomial")

Ich bekomme eine Warnmeldung

Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

aber eine sehr vernünftige Übereinstimmung (alle Faktoren sind kategorisch, so dass ich leicht überprüfen kann, ob die Modellvorhersagen nahe am subjektübergreifenden Mittelwert liegen und ob sie es sind).

Was ich aber eigentlich nutzen möchte ist

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

Es gibt mir die identische Warnung, gibt ein Modell zurück, aber dieses Modell ist eindeutig sehr abweichend. die schätzungen der fixen effekte sind sehr weit von glm()jenen und den themenübergreifenden mitteln entfernt. (Und ich muss glmerControl(optimizer="bobyqa")in den glmerAufruf einbeziehen, sonst konvergiert er überhaupt nicht.)

Amöbe sagt Reinstate Monica
quelle
1
Wie wäre es zuerst die Wahrscheinlichkeiten zu transformieren? Können Sie mit einer logit-Transformation etwas erreichen, das näher an der Normalverteilung liegt? Oder das ArcSin-Quadrat? Das wäre mir lieber, als glmer zu verwenden. In Ihrer Hack-Lösung können Sie auch versuchen, für jede Beobachtung einen Zufallseffekt hinzuzufügen, um die Unterdispersion aufgrund der von Ihnen gewählten Gewichte zu berücksichtigen.
Aaron - Wiedereinsetzung von Monica
Vielen Dank. Ja, ich kann die DV anmelden und dann das gemischte Gauß-Modell (lmer) verwenden, aber dies ist auch eine Art Hack, und ich habe gelesen, dass es nicht empfohlen wird. Ich werde für jede Beobachtung einen zufälligen Effekt ausprobieren! Im Moment versuche ich es mit einem Beta-Mischmodell. lme4 kann damit nicht umgehen, aber glmmadmb kann es. Wenn ich laufe glmmadmb(reportedProbability ~ a + b + c + (1 | subject), myData, family="beta"), erhalte ich eine korrekte Anpassung und angemessene Konfidenzintervalle, aber eine Konvergenz- Warnung ist fehlgeschlagen : Beta könnte für mich funktionieren, da ich keine Fälle mit DV = 0 oder DV = 1 habe.
Amöbe sagt Reinstate Monica
Ich weiß nicht für glmer, aber für glm kann dies helfen: stats.stackexchange.com/questions/164120/… :
1
@Aaron: Ich habe versucht + (1 | rowid), meinen glmer-Aufruf zu ergänzen, und dies ergibt stabile Schätzungen und stabile Konfidenzintervalle, unabhängig von meiner Gewichtswahl (ich habe 100 und 500 versucht). Ich habe auch versucht, lmer auf logit auszuführen (gemeldete Wahrscheinlichkeit) und ich bekomme fast genau das Gleiche. Beide Lösungen scheinen also gut zu funktionieren! Beta MM mit glmmadmb liefert ebenfalls sehr nahe beieinander liegende Ergebnisse, die jedoch aus irgendeinem Grund nicht vollständig konvergieren und ewig dauern. Stellen Sie eine Antwort auf diese Optionen und erläutern Sie ein wenig die Unterschiede und Vor- und Nachteile! (Vertrauensintervalle, die ich erwähne, sind alle Wald.)
Amöbe sagt Reinstate Monica
1
Und sie sind sich absolut sicher über ihren Wert wie 0,9, oder haben sie auch eine "Fehlerspanne"? Können Sie davon ausgehen, dass das von verschiedenen Probanden gemeldete Vertrauen gleich genau ist?

Antworten:

20

Es ist sinnvoll, mit einem einfacheren Fall ohne zufällige Effekte zu beginnen.

Es gibt vier Möglichkeiten, um mit einer kontinuierlichen Null-zu-Eins-Antwortvariablen umzugehen, die sich wie ein Bruchteil oder eine Wahrscheinlichkeit verhält ( dies ist unser kanonischster, am besten bewerteter und angesehener Thread zu diesem Thema, aber leider werden dort nicht alle vier Optionen behandelt):

  1. p=m/nnnN

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

    betareg(p ~ a+b+c, myData)
  3. Logit transformiert die Antwort und verwendet die lineare Regression. Dies wird normalerweise nicht empfohlen.

    lm(log(p/(1-p)) ~ a+b+c, myData)
  4. Passen Sie ein Binomialmodell an, aber berechnen Sie dann die Standardfehler unter Berücksichtigung der Überdispersion. Die Standardfehler können auf verschiedene Arten berechnet werden:

    • (a) skalierte Standardfehler über die Überdispersionsschätzung ( eins , zwei ). Dies wird als "Quasi-Binomial" -GLM bezeichnet.

    • (b) robuste Standardfehler über den Sandwich-Schätzer ( eins , zwei , drei , vier ). Dies wird in der Ökonometrie als "Fractional Logit" bezeichnet.


    Die Buchstaben (a) und (b) sind nicht identisch (siehe diesen Kommentar und die Abschnitte 3.4.1 und 3.4.2 in diesem Buch sowie diesen SO-Beitrag und auch diesen und diesen ), führen jedoch tendenziell zu ähnlichen Ergebnissen. Option (a) wird glmwie folgt umgesetzt:

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

Die gleichen vier Möglichkeiten stehen mit zufälligen Effekten zur Verfügung.

  1. Mit weightsArgument ( eins , zwei ):

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

    Gemäß dem zweiten obigen Link ist es möglicherweise eine gute Idee, eine Überdispersion zu modellieren (siehe dort und auch # 4 unten).

  2. Betamischmodell verwenden:

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

    oder

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

    Wenn die Antwortdaten exakte Nullen oder Einsen enthalten, kann man das auf Null / Eins aufgeblasene Betamodell in verwenden glmmTMB.

  3. Verwenden der logit-Transformation der Antwort:

    lmer(log(p/(1-p)) ~ a+b+c + (1|subject), myData)
  4. Berücksichtigung der Überdispersion im Binomialmodell. Dies verwendet einen anderen Trick: Hinzufügen eines zufälligen Effekts für jeden Datenpunkt:

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

    Aus irgendeinem Grund funktioniert dies nicht richtig, da sich glmer()Beschwerden über nicht ganzzahlige Werte pund Unsinnschätzungen ergeben. Eine Lösung, die ich mir ausgedacht habe, ist, eine gefälschte Konstante zu verwenden weights=kund sicherzustellen, dass diese p*kimmer eine Ganzzahl ist. Dies erfordert eine Rundung, paber wenn Sie keine ausreichend große Auswahl treffen , sollte dies keine große Rolle spielen. Die Ergebnisse scheinen nicht vom Wert von abzuhängen 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"))

    Späteres Update (Januar 2018): Dies ist möglicherweise ein ungültiger Ansatz. Siehe Diskussion hier . Ich muss das genauer untersuchen.


In meinem speziellen Fall ist Option 1 nicht verfügbar.

Option 2 ist sehr langsam und hat Probleme mit der Konvergenz: Die glmmadmbAusführung dauert fünf bis zehn Minuten (und beklagt sich immer noch, dass sie nicht konvergiert ist!), Während sie lmerin glmerSekundenbruchteilen funktioniert und einige Sekunden dauert. Update: Ich habe versucht, glmmTMBwie in den Kommentaren von @BenBolker vorgeschlagen, und es funktioniert fast so schnell wie glmerohne Konvergenzprobleme. Das ist es also, was ich verwenden werde.

Die Optionen 3 und 4 liefern sehr ähnliche Schätzungen und sehr ähnliche Wald-Konfidenzintervalle (erhalten mit confint). Ich bin allerdings kein großer Fan von # 3, weil es eine Art Betrug ist. Und # 4 fühlt sich etwas abgefahren an.

Vielen Dank an @Aaron, der mich in seinem Kommentar auf # 3 und # 4 hingewiesen hat.

Amöbe sagt Reinstate Monica
quelle
1
Schöne Antwort, gut erklärt und verbunden mit den No-Random-Effects-Modellen. Ich würde nicht sagen, dass # 3 (die Transformation) betrügt, diese Art von Transformationen sind in solchen Analysen sehr verbreitet. Ich würde stattdessen sagen, dass sowohl # 3 als auch # 4 Annahmen über die Beziehung zur Verteilung der Daten treffen und so auch über die Beziehung zwischen Mittelwert und Varianz, und nur weil # 4 auf der Skala der Daten modelliert wurde gesammelt am bedeutet nicht, dass diese Annahmen besser sein werden.
Aaron - Reinstate Monica
1
Bei # 3 wird davon ausgegangen, dass das Logit der Wahrscheinlichkeiten bei konstanter Varianz normal ist, während bei # 4 davon ausgegangen wird, dass die Varianz proportional zu p (1-p) ist. Nach Ihrer Beschreibung der Passung scheinen diese ähnlich genug zu sein, um nicht zu viel zu bedeuten. Und # 3 ist mit ziemlicher Sicherheit mehr Standard (abhängig von Ihrer Zielgruppe). Wenn die Diagnose angemessen ist, ist dies diejenige, die ich bevorzugen würde.
Aaron - Reinstate Monica
1
Eine andere Möglichkeit ist die Verwendung von glmmTMB . nachdem mit der Installation devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB"), Verwendung glmmTMB(p ~ a+b+c + (1|subject), myData, family=list(family="beta",link="logit"))sollte funktionieren ...
Ben Bolker
@BenBolker Danke! Gibt es einen Grund, glmmTMB gegenüber glmmADMB (für Betamodelle) vorzuziehen oder umgekehrt? Ist eines dieser Pakete neuer oder aktiver entwickelt? Abgesehen davon, darf ich fragen, welcher Ansatz unter den in dieser Antwort aufgelisteten - Gaußsches Glmm nach logit-Transformation, Beta-Glmm oder Binomial-Glmm mit (1 | rowid) -Term - ist Ihnen im Allgemeinen vorzuziehen?
Amöbe sagt Reinstate Monica
1
Ich bevorzuge das Beta-GLMM, wenn es machbar ist - es ist das statistische Modell, das Änderungen der Proportionen über Kovariaten / Gruppen hinweg messen soll . glmmTMBist schneller und stabiler als glmmADMBund unter (leicht) aktiverer Entwicklung, wenn auch nicht so ausgereift.
Ben Bolker