¿Cómo lidiar con la separación casi completa en un GLMM logístico?

Actualización : dado que ahora sé que mi problema se llama separación casi completa , actualicé la pregunta para reflejar esto (gracias a Aaron).

Tengo un conjunto de datos de un experimento en el que 29 participantes humanos (factor code) trabajaron en un conjunto de ensayos y responsefue 1 o 0. Además, manipulamos los materiales para tener tres factores cruzados p.validity(válido versus no válido), type(afirmación versus negación) y counterexamples(pocos versus muchos):

d.binom <- read.table("http://pastebin.com/raw.php?i=0yDpEri8")
str(d.binom)
## 'data.frame':   464 obs. of  5 variables:
##      $ code           : Factor w/ 29 levels "A04C","A14G",..: 1 1 1 1 1 1 1 1 1 1 ...
##      $ response       : int  1 1 1 1 0 1 1 1 1 1 ...
##      $ counterexamples: Factor w/ 2 levels "few","many": 2 2 1 1 2 2 2 2 1 1 ...
##      $ type           : Factor w/ 2 levels "affirmation",..: 1 2 1 2 1 2 1 2 1 2 ...
##      $ p.validity     : Factor w/ 2 levels "invalid","valid": 1 1 2 2 1 1 2 2 1 1 ...

En general, solo hay un pequeño número de ceros:

mean(d.binom$response)
## [1] 0.9504

Una hipótesis es que hay un efecto de validity, sin embargo, el análisis preliminar sugiere que podría haber un efecto de counterexamples. Como tengo datos dependientes (cada participante trabajó en todos los ensayos) me gustaría usar un GLMM en los datos. Desafortunadamente, counterexamplessepare casi por completo los datos (al menos para un nivel):

with(d.binom, table(response, counterexamples))
##         counterexamples
## response few many
##        0   1   22
##        1 231  210

Esto también se refleja en el modelo:

require(lme4)
options(contrasts=c('contr.sum', 'contr.poly'))


m2 <- glmer(response ~ type * p.validity * counterexamples + (1|code), 
            data = d.binom, family = binomial)
summary(m2)
## [output truncated]
## Fixed effects:
##                                      Estimate Std. Error z value Pr(>|z|)
##   (Intercept)                            9.42     831.02    0.01     0.99
##   type1                                 -1.97     831.02    0.00     1.00
##   p.validity1                            1.78     831.02    0.00     1.00
##   counterexamples1                       7.02     831.02    0.01     0.99
##   type1:p.validity1                      1.97     831.02    0.00     1.00
##   type1:counterexamples1                -2.16     831.02    0.00     1.00
##   p.validity1:counterexamples1           2.35     831.02    0.00     1.00
##   type1:p.validity1:counterexamples1     2.16     831.02    0.00     1.00

Los errores estándar para los parámetros son simplemente una locura. Como mi objetivo final es evaluar si ciertos efectos son significativos o no, los errores estándar no son totalmente sin importancia.

• ¿Cómo puedo lidiar con la separación casi completa? Lo que quiero es obtener estimaciones a partir de las cuales pueda juzgar si un cierto efecto es significativo o no (por ejemplo, usando el PRmodcomppaquete pkrtest, pero este es otro paso que no se describe aquí).

Los enfoques que utilizan otros paquetes también están bien.

Me temo que hay un error tipográfico en su título: no debe intentar ajustar modelos mixtos, y mucho menos modelos mixtos no lineales, con solo 30 clústeres. No, a menos que crea que puede ajustar una distribución normal a 30 puntos obstruidos por error de medición, no linealidades y una separación casi completa (también conocida como predicción perfecta).

Lo que haría aquí es ejecutar esto como una regresión logística regular con la corrección de Firth :

library(logistf)
mf <- logistf(response ~ type * p.validity * counterexamples + as.factor(code),
      data=d.binom)

La corrección de Firth consiste en agregar una penalización a la probabilidad, y es una forma de contracción. En términos bayesianos, las estimaciones resultantes son los modos posteriores del modelo con un Jeffreys previo. En términos frecuentistas, la penalización es el determinante de la matriz de información correspondiente a una sola observación y, por lo tanto, desaparece asintóticamente.

Puede utilizar un enfoque Bayesiano de máximo a posteriori con un previo débil sobre los efectos fijos para obtener aproximadamente el mismo efecto. En particular, el paquete blme para R (que es una envoltura delgada alrededor del lme4paquete) hace esto, si especifica los previos para los efectos fijos como en el ejemplo aquí (busque "separación completa"):

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

Este ejemplo es de un experimento donde ttthay un efecto fijo categórico con 4 niveles, por lo que el$\beta$ el vector tendrá una longitud 4. La matriz de varianza-covarianza previa especificada es $\Sigma = 9 I$, es decir, los parámetros de efectos fijos tienen independiente $N(\mu=0,\sigma^2=9)$ (o $\sigma$, es decir, la devaluación estándar, $=3$) anteriores. Esto funciona bastante bien, aunque no es idéntico a la corrección de Firth (ya que Firth corresponde a un Jeffreys anterior , que no es lo mismo).

El ejemplo vinculado muestra que también puede hacerlo con el MCMCglmmpaquete, si desea ir completamente bayesiano ...