Binomial glmm mit einer kategorialen Variablen mit vollem Erfolg

11

Ich führe ein glmm mit einer binomialen Antwortvariablen und einem kategorialen Prädiktor aus. Der zufällige Effekt ergibt sich aus dem verschachtelten Design, das für die Datenerfassung verwendet wird. Die Daten sehen folgendermaßen aus:

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

Das erste Modell, das ich laufe, sieht so aus

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

Ich erhalte zwei Warnungen, die so aussehen:

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?

Die Modellzusammenfassung zeigt, dass eine der Behandlungen einen ungewöhnlich großen Standardfehler aufweist, den Sie hier sehen können:

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 

Ich habe die verschiedenen Optimierer von glmer control und Funktionen von anderen Paketen ausprobiert und erhalte eine ähnliche Ausgabe. Ich habe das Modell mit glm ausgeführt und den zufälligen Effekt ignoriert, und das Problem besteht weiterhin. Während ich die Daten untersuchte, stellte ich fest, dass die Behandlung mit einem hohen Std. Fehler hat nur Erfolge in der Antwortvariablen. Um zu überprüfen, ob dies das Problem verursachen könnte, habe ich einen gefälschten Datenpunkt mit einem "Fehler" für diese Behandlung hinzugefügt, und das Modell läuft reibungslos und gibt einen angemessenen Standardfehler aus. Das können Sie hier sehen:

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 .

Ich habe mich gefragt, ob meine Intuition in Bezug auf das Fehlen von Fehlern bei dieser Behandlung, die eine gute Einschätzung verhindern, richtig ist, und wie ich dieses Problem umgehen kann.

Danke im Voraus!

Amöbe sagt Reinstate Monica
quelle

Antworten:

15

Ihre Intuition ist genau richtig. Dieses Phänomen wird als vollständige Trennung bezeichnet . Sie können ziemlich viel finden (jetzt, wo Sie den Namen kennen). Googeln ... Es wird hier in einem allgemeinen Kontext und hier im Kontext von GLMMs ziemlich gründlich diskutiert . Die Standardlösung für dieses Problem besteht darin, einen kleinen Term hinzuzufügen, der die Parameter auf Null zurückschiebt. In häufig auftretenden Kontexten wird dies als bestrafte oder vorurteilskorrigierte Methode bezeichnet. Der Standardalgorithmus basiert auf Firth (1993, "Bias Reduction of Maximum Likelihood Estimations ", Biometrika 80, 27-38) und ist im logistf-Paket implementiertauf CRAN. In Bayes'schen Kontexten wird dies als Hinzufügen einer Schwachstelle vor den Parametern mit festem Effekt dargestellt.

Meines Wissens wurde der Algorithmus von Firth nicht auf GLMMs erweitert, aber Sie können den Bayes'schen Trick verwenden, indem Sie das blme- Paket verwenden, das eine dünne Bayes'sche Schicht über das Paket legt lme4. Hier ist ein Beispiel aus der oben verlinkten GLMM-Diskussion:

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

Die ersten beiden Zeilen in diesem Beispiel sind genau die gleichen, die wir im Standardmodell glmerverwenden würden. Der letzte gibt an, dass der Prior für die festen Effekte eine multivariate Normalverteilung mit einer diagonalen Varianz-Kovarianz-Matrix ist. Die Matrix ist 4x4 (weil wir in diesem Beispiel 4 Parameter mit festem Effekt haben), und die vorherige Varianz jedes Parameters ist 9 (entsprechend einer Standardabweichung von 3, was ziemlich schwach ist - das heißt, +/- 2SD ist ( -6,6), was ein sehr großer Bereich auf der Logit-Skala ist).

Die sehr großen Standardfehler der Parameter in Ihrem Beispiel sind ein Beispiel für ein Phänomen, das eng mit der vollständigen Trennung zusammenhängt (es tritt immer dann auf, wenn in einem logistischen Modell extreme Parameterwerte auftreten), der als Hauck-Donner-Effekt bezeichnet wird .

Zwei weitere potenziell nützliche Referenzen (ich habe mich selbst noch nicht damit befasst):

  • Gelman A, Jakulin A, Pittau MG und Su TS (2008) Eine schwach informative Standard-Vorverteilung für logistische und andere Regressionsmodelle. Annals of Applied Statistics , 2, 1360–383.
  • José Cortiñas Abrahantes und Marc Aerts (2012) Eine Lösung zur Trennung für gruppierte Binärdaten Statistical Modeling 12 (1): 3–27 doi: 10.1177 / 1471082X1001200102

Eine neuere Google-Gelehrten-Suche nach "bglmer 'vollständige Trennung'" findet:

  • Quiñones, AE und WT Wcislo. "Cryptic Extended Brood Care in der fakultativ eusozialen Schweißbiene Megalopta genalis ." Insectes Sociaux 62.3 (2015): 307–313.
Ben Bolker
quelle
wow vielen dank !! Das macht durchaus Sinn und das Modell läuft jetzt reibungslos mit dem bglmer. Ich hätte nur noch eine Frage: Kann ich die Methoden wie in lme4 verwenden, um die zufälligen und festen Effekte zu bewerten, mit anderen Worten, um verschiedene Modelle zu vergleichen?
2
Ich würde es sagen, aber ich weiß nicht, ob es formelle und / oder von Experten geprüfte Unterstützung für meine Meinung gibt ...
Ben Bolker
Vielen Dank! Das ist auch genau mein Problem. Ein kurzes Follow-up: Im Gegensatz zu Ihrem Beispiel, das einen Faktor mit 4 Ebenen hat, habe ich ein 2 x 2-Design, bei dem jeder Faktor 2 Ebenen hat (die Summe beträgt also immer noch 4 Ebenen). Kann ich für mein Modell auch diag (9,4) verwenden? Ich bin mit Matrizen nicht vertraut, deshalb wollte ich es noch einmal überprüfen. Um diese Lösung in meinem Artikel zu rechtfertigen, sollte ich Firth (1993) zitieren, oder gibt es einen direkteren Artikel, der Ihre Lösung mit bglmer () implementiert hat?
Sol
2
siehe aktualisierte Antwort.
Ben Bolker
2
Ich denke schon - es sollte nur eine Rolle spielen, wie viele feste Effektparameter es insgesamt gibt.
Ben Bolker