Binomial glmm con una variable categórica con éxitos completos

Estoy ejecutando un glmm con una variable de respuesta binomial y un predictor categórico. El efecto aleatorio viene dado por el diseño anidado utilizado para la recopilación de datos. Los datos se ven así:

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

El primer modelo que corro se ve así

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

Recibo dos advertencias que se ven así:

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?

El resumen del modelo muestra que uno de los tratamientos tiene un error estándar inusualmente grande, que puede ver aquí:

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 

Probé los diferentes optimizadores del control glmer y las funciones de otros paquetes, y obtuve un resultado similar. He ejecutado el modelo usando glm ignorando el efecto aleatorio, y el problema persiste. Mientras exploraba los datos me di cuenta de que el tratamiento con un alto estándar. El error solo tiene éxito en la variable de respuesta. Solo para verificar si eso podría estar causando el problema, agregué un punto de datos falso con un "fallo" para ese tratamiento y el modelo funciona sin problemas y proporciona un error estándar razonable. Se puede ver que aquí:

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 .

Me preguntaba si mi intuición es correcta acerca de la falta de fracasos para ese tratamiento que impide una buena estimación, y cómo puedo solucionar este problema.

¡Gracias por adelantado!

Tu intuición es exactamente la correcta. Este fenómeno se llama separación completa . Puede encontrar bastante (ahora que conoce su nombre) buscando en Google ... Se trata bastante a fondo aquí en un contexto general , y aquí en el contexto de los GLMM . La solución estándar a este problema es agregar un término pequeño que empuje los parámetros hacia cero; en contextos frecuentes, esto se llama un método penalizado o con corrección de sesgos . El algoritmo estándar se debe a Firth (1993, "Reducción de sesgo de las estimaciones de máxima verosimilitud" Biometrika 80, 27-38), y se implementa en el paquete logistfen CRAN. En contextos bayesianos, esto se enmarca como la adición de un débil antes de los parámetros de efectos fijos.

Que yo sepa, el algoritmo de Firth no se ha extendido a GLMM, pero puedes usar el truco bayesiano usando el paquete blme , que coloca una delgada capa bayesiana sobre la parte superior del lme4paquete. Aquí hay un ejemplo de la discusión GLMM vinculada anteriormente:

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

Las dos primeras líneas en este ejemplo son exactamente las mismas que usaríamos en el glmermodelo estándar ; el último especifica que el previo para los efectos fijos es una distribución normal multivariada con una matriz de varianza-covarianza diagonal. La matriz es 4x4 (porque tenemos 4 parámetros de efectos fijos en este ejemplo), y la varianza anterior de cada parámetro es 9 (correspondiente a una desviación estándar de 3, que es bastante débil, lo que significa que +/- 2SD es ( -6,6), que es un rango muy grande en la escala logit).

Los errores estándar muy grandes de los parámetros en su ejemplo son un ejemplo de un fenómeno estrechamente relacionado con la separación completa (ocurre cada vez que obtenemos valores de parámetros extremos en un modelo logístico) llamado efecto Hauck-Donner .

Dos referencias potencialmente más útiles (aún no las he investigado):

• Gelman A, Jakulin A, Pittau MG y Su TS (2008) Una distribución previa por defecto débilmente informativa para modelos logísticos y otros modelos de regresión. Anales de Estadística Aplicada , 2, 1360–383.
• José Cortiñas Abrahantes y Marc Aerts (2012) Una solución para la separación de datos binarios agrupados Modelización estadística 12 (1): 3–27 doi: 10.1177 / 1471082X1001200102

Una búsqueda más reciente de Google Scholar para "bglmer 'separación completa'" encuentra:

• Quiñones, AE y WT Wcislo. "Cuidado criado extendido de cría en la Abeja de sudor facultativamente eusocial Megalopta genalis ". Insectes Sociaux 62.3 (2015): 307–313.