Dwumianowy glmm ze zmienną kategorialną z pełnymi sukcesami

Korzystam z glmm ze zmienną dwumianową i predyktorem jakościowym. Efekt losowy wynika z zagnieżdżonego projektu stosowanego do gromadzenia danych. Dane wyglądają następująco:

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

Pierwszy model, który uruchamiam, wygląda tak

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

Otrzymuję dwa ostrzeżenia, które wyglądają tak:

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?

Podsumowanie modelu pokazuje, że jeden z zabiegów ma niezwykle duży błąd standardowy, który można zobaczyć tutaj:

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 

Wypróbowałem różne optymalizatory z kontroli glittera i funkcji z innych pakietów i otrzymałem podobny wynik. Uruchomiłem model przy użyciu glm, ignorując efekt losowy, a problem nadal występuje. Podczas eksploracji danych zdałem sobie sprawę, że leczenie z wysokim Std. błąd ma tylko sukcesy w zmiennej odpowiedzi. Aby sprawdzić, czy może to powodować problem, dodałem fałszywy punkt danych z „niepowodzeniem” dla tego leczenia, a model działa płynnie i daje rozsądny błąd standardowy. Możesz to zobaczyć tutaj:

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 .

Zastanawiałem się, czy moja intuicja ma rację co do braku niepowodzeń tego leczenia uniemożliwiającego dobre oszacowanie i jak mogę obejść ten problem.

Z góry dziękuję!

Twoja intuicja jest dokładnie właściwa. Zjawisko to nazywa się całkowitą separacją . Możesz znaleźć całkiem sporo (teraz, gdy znasz jego nazwę) Googling wokoło ... Jest dość dokładnie omawiany tutaj w ogólnym kontekście , a tutaj w kontekście GLMM . Standardowym rozwiązaniem tego problemu jest dodanie małego terminu, który przesuwa parametry z powrotem do zera - w kontekstach częstych nazywa się to metodą karaną lub skorygowaną o błąd systematyczny. Standardowy algorytm wynika z Firtha (1993, „Redukcja błędu maksymalnego prawdopodobieństwa szacunków” Biometrika 80, 27-38) i jest zaimplementowany w pakiecie logistfna CRAN. W kontekście bayesowskim jest to sformułowane jako dodanie osłabienia przed parametrami o stałym efekcie.

O ile wiem, algorytm Firtha nie został rozszerzony na GLMM, ale możesz użyć sztuczki bayesowskiej, używając pakietu blme , który nakłada cienką warstwę bayesowską na wierzch lme4pakietu. Oto przykład z wyżej wymienionej dyskusji GLMM:

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

Pierwsze dwie linie w tym przykładzie są dokładnie takie same, jak w standardowym glmermodelu; ostatni określa, że ​​wcześniej dla ustalonych efektów jest wielowymiarowy rozkład normalny z diagonalną macierzą wariancji-kowariancji. Matryca to 4x4 (ponieważ w tym przykładzie mamy 4 parametry o stałym efekcie), a wcześniejsza wariancja każdego parametru wynosi 9 (co odpowiada standardowemu odchyleniu 3, które jest dość słabe - oznacza to, że +/- 2SD wynosi ( -6,6), co stanowi bardzo duży zakres w skali logit).

Bardzo duże standardowe błędy parametrów w twoim przykładzie są przykładem zjawiska ściśle związanego z całkowitą separacją (występuje za każdym razem, gdy otrzymujemy ekstremalne wartości parametrów w modelu logistycznym) o nazwie efekt Haucka-Donnera .

Dwa kolejne potencjalnie przydatne odniesienia (sam jeszcze ich nie kopałem):

• Gelman A, Jakulin A, Pittau MG i Su TS (2008) Słabo informacyjna domyślna wcześniejsza dystrybucja dla modeli logistycznych i innych modeli regresji. Annals of Applied Statistics , 2, 1360–383.
• José Cortiñas Abrahantes i Marc Aerts (2012) Rozwiązanie separacji dla klastrowych danych binarnych Modelowanie statystyczne 12 (1): 3–27 doi: 10.1177 / 1471082X1001200102

W nowszym wyszukiwaniu Google dla słowa „bglmer” zupełna separacja ”znaleziono:

• Quiñones, AE i WT Wcislo. „Cryptic Extended Brood Care in Facective Eusocial Sweat Bee Megalopta genalis .” Insectes Sociaux 62.3 (2015): 307–313.