¿Cómo encajar el GLMM binomial con una respuesta continua entre 0 y 1 que no sea una fracción?

Question

Estoy tratando de usar lme4::glmer() para ajustar un GLMM binomial con una variable dependiente que no es binaria, sino una variable continua entre cero y uno. Se puede pensar en esta variable como una probabilidad; de hecho es probabilidad como la reportada por los sujetos humanos (en un experimento que ayudo a analizar). El glmer() produce un modelo que está claramente fuera de lugar, y muy lejos del que tengo con glm() así que algo va mal. ¿Por qué? ¿Qué puedo hacer?

---

Más detalles

Aparentemente es posible utilizar la regresión logística no sólo para la DV binaria sino también para la DV continua entre cero y uno. De hecho, cuando corro

glm(reportedProbability ~ a + b + c, myData, family="binomial")

Recibo un mensaje de advertencia

Warning message:
In eval(expr, envir, enclos) : non-integer #successes in a binomial glm!

pero un ajuste muy razonable (todos los factores son categóricos, por lo que puedo comprobar fácilmente si las predicciones de los modelos se acercan a los medios de los sujetos, y así es).

Sin embargo, lo que realmente quiero usar es

glmer(reportedProbability ~ a + b + c + (1 | subject), myData, family="binomial")

Me da la misma advertencia, devuelve un modelo, pero este modelo está claramente muy fuera de lugar; las estimaciones de los efectos fijos están muy lejos de la glm() y de los medios de la materia. (Y necesito incluir glmerControl(optimizer="bobyqa") en el glmer llamada, de lo contrario no converge en absoluto.)

Answer 1

Tiene sentido empezar con un caso más simple sin efectos aleatorios.

Hay cuatro maneras de tratar con la variable de respuesta continua de cero a uno que se comporta como una fracción o una probabilidad:

1. Si es una fracción $p=m/n$ y todos $n$ se conocen, entonces se puede usar la regresión logística estándar, alias binomio GLM. Una forma de codificarlo en R es (asumiendo que n es un vector de $n$ para cada punto de datos):

   glm(p ~ a+b+c, myData, family="binomial", weights=n)

2. Si $p$ no es una fracción, entonces se puede usar la regresión beta. Esto sólo funcionará si $p$ nunca es igual a $0$ o $1$ . Si es así, entonces son posibles modelos beta más complicados de cero/uno-inflado, pero esto se vuelve más complicado.

   betareg(p ~ a+b+c, myData)

3. Logit transformar la respuesta y utilizar la regresión lineal. Esto normalmente no se aconseja.

   lm(log(p/(1-p)) ~ a+b+c, myData)

4. Ajustar un modelo binomial pero luego calcular los errores estándar teniendo en cuenta la sobredispersión. Los errores estándar pueden ser calculados de varias maneras:

   • a) Errores estándar escalados mediante la estimación de la sobredispersión ( uno , dos ). Esto se llama "pseudobinomio" en la comunidad R;

   • b) Errores estándar robustos mediante el estimador de sándwiches ( uno , dos , tres , cuatro ). Esto se llama "logit fraccionario" en la econometría;

   • c) Tal vez algunos otros enfoques sólidos.

   Los apartados a) y b) no son idénticos (véase este comentario y las secciones 3.4.1 y 3.4.2 de este libro y este puesto de SO y también este y este ), pero tienden a dar resultados similares. La opción a) se aplica en glm de la siguiente manera:

   glm(p ~ a+b+c, myData, family="pseudobinomial")

---

Las mismas cuatro formas están disponibles con efectos aleatorios.

1. Usando weights argumento ( uno , dos ):

   glmer(p ~ a+b+c + (1|subject), myData, family="binomial", weights=n)

   De acuerdo con el segundo enlace de arriba, podría ser una buena idea modelar la sobredispersión, ver allí (y también el #4 abajo).

2. Usando un modelo mixto beta:

   glmmadmb(p ~ a+b+c + (1|subject), myData, family="beta")

   o

   devtools::install_github("glmmTMB/glmmTMB",sub="glmmTMB")
   glmmTMB(p ~ a+b+c + (1|subject), myData, 
           family=list(family="beta",link="logit"))

3. Usando la transformación de logit de la respuesta:

   lmer(log(p/(1-p)) ~ a+b+c + (1|subject), myData)

4. Teniendo en cuenta la sobredispersión en el modelo del binomio. Esto utiliza un truco diferente: añadir un efecto aleatorio para cada punto de datos:

   myData$rowid = as.factor(1:nrow(myData))
   glmer(p ~ a+b+c + (1|subject) + (1|rowid), myData, family="binomial",
         glmerControl(optimizer="bobyqa"))

   Por alguna razón esto no funciona correctamente como glmer() se queja de los no enteros p y produce estimaciones sin sentido. Una solución que se me ocurrió es usar una constante falsa weights y asegurarse de que p*n siempre es un número entero. Esto requiere redondear p pero seleccionando n que es lo suficientemente grande no debería importar mucho. Los resultados no parecen depender del valor de n .

   n = 100
   glmer(round(p*n)/n ~ a+b+c + (1|subject) + (1|rowid), myData, 
         family="binomial", weights=rowid*0+n, glmerControl(optimizer="bobyqa"))

---

En mi caso específico la opción #1 no está disponible.

La opción 2 es muy lenta y tiene problemas de convergencia: glmmadmb tarda cinco-diez minutos en funcionar (¡y todavía se queja de que no convergió!), mientras que lmer funciona en una fracción de segundo y glmer toma un par de segundos. [ Actualizar: Lo intenté. glmmTMB como se sugiere en los comentarios de @BenBolker y funciona casi tan rápido como glmer sin problemas de convergencia. Podría terminar usando esta opción].

Las opciones 3 y 4 arrojan estimaciones muy similares e intervalos de confianza de Wald muy parecidos (obtenidos con confint ). No soy un gran fan del número 3 porque es una especie de trampa. Así que probablemente usaré el número 4.

Enorme gracias a @Aaron que me señaló hacia el #3 y #4 en su comentario.