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

Question

Actualización: Dado que ahora sé que mi problema se llama separación cuasi 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 pruebas y la response fue 1 ó 0. Además, manipulamos los materiales de manera que teníamos tres factores cruzados, p.validity (válido versus invá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, hay solo un pequeño número de 0s:

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

Una hipótesis es que existe un efecto de validity, sin embargo, el análisis preliminar sugiere que podría haber un efecto de counterexamples. Dado que tengo datos dependientes (cada participante trabajó en todas las pruebas), me gustaría usar un GLMM en los datos. Desafortunadamente, counterexamples separa cuasi-completamente 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]
## Efectos fijos:
##                                      Estimado Error estándar Valor z Pr(>|z|)
##   (Intercepto)                            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 de los parámetros son simplemente insensatos. Dado que mi objetivo final es evaluar si ciertos efectos son significativos, los errores estándar no son totalmente irrelevantes.

• ¿Cómo puedo lidiar con la separación cuasi 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 PRmodcomp del paquete pkrtest, pero este es otro paso no descrito aquí).

Enfoques utilizando otros paquetes también son válidos.

Answer 1

Tengo miedo de que haya un error tipográfico en tu título: no deberías intentar ajustar modelos mixtos, y mucho menos modelos mixtos no lineales, con solo 30 clústeres. A menos que creas que puedes ajustar una distribución normal a 30 puntos obstruidos por errores 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 verosimilitud, y es una forma de contracción. En términos bayesianos, los estimados resultantes son los modos posteriores del modelo con una priori de Jeffreys. 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.

Answer 2

Puedes utilizar un enfoque máximo a posteriori de Bayes con un prior débil en los efectos fijos para obtener aproximadamente el mismo efecto. En particular, el paquete blme para R (que es un envoltorio delgado alrededor del paquete lme4) hace esto, si especificas priors para los efectos fijos como en el ejemplo aquí (busca "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 ttt es un efecto fijo categórico con 4 niveles, por lo que el vector $\beta$ tendrá longitud 4. La matriz de varianza-covarianza prior especificada es $\Sigma = 9 I$, es decir, los parámetros de efecto fijo tienen priors independientes $N(\mu=0,\sigma^2=9)$ (o $\sigma$, es decir, desviación estándar, $=3$). Esto funciona bastante bien, aunque no es idéntico a la corrección de Firth (ya que Firth corresponde a un prior de Jeffreys, que no es exactamente lo mismo).

El ejemplo enlazado muestra que también puedes hacerlo con el paquete MCMCglmm, si quieres ir completamente bayesiano...