¿Cómo ajustar un modelo mixto con una respuesta variable entre 0 y 1?

Estoy tratando de usar lme4::glmer() para ajustar un modelo mixto binomial generalizado (GLMM) con una variable dependiente que no es binaria, sino una variable continua entre cero y uno. Uno puede pensar en esta variable como una probabilidad; de hecho, es la probabilidad según lo informado por sujetos humanos (en un experimento que ayudo a analizar). Es decir, no es una fracción "discreta", sino una variable continua.

Mi glmer()llamada no funciona como se esperaba (ver más abajo). ¿Por qué? ¿Que puedo hacer?

Edición posterior: mi respuesta a continuación es más general que la versión original de esta pregunta, por lo que modifiqué la pregunta para que sea más general también.

Más detalles

Aparentemente, es posible utilizar la regresión logística no solo para DV binario sino también para DV continuo 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 verificar fácilmente si las predicciones del modelo están cerca de las medias entre sujetos, y lo están).

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 apagado; Las estimaciones de los efectos fijos están muy lejos de glm()las de los medios transversales. (Y necesito incluir glmerControl(optimizer="bobyqa")en la glmerllamada, de lo contrario no converge en absoluto).

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

Hay cuatro formas de lidiar con la variable de respuesta continua de cero a uno que se comporta como una fracción o una probabilidad ( este es nuestro hilo más canónico / votado / visto sobre este tema, pero desafortunadamente no se discuten las cuatro opciones allí):

1. $p=m/n$$n$n$N$

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

2. $p$$p$$0$$1$

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

3. Logit transforma la respuesta y usa regresión lineal. Esto generalmente no se aconseja.

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

4. Ajuste un modelo binomial pero luego calcule los errores estándar teniendo en cuenta la dispersión excesiva. Los errores estándar se pueden calcular de varias maneras:

   • (a) errores estándar escalados a través de la estimación de sobredispersión ( uno , dos ). Esto se llama GLM "cuasi-binomial".

   • (b) errores estándar robustos a través del estimador sandwich ( uno , dos , tres , cuatro ). Esto se llama "logit fraccional" en econometría.

   
   Las letras (a) y (b) no son idénticas (vea este comentario , y las secciones 3.4.1 y 3.4.2 en este libro , y esta publicación SO y también esta y esta ), pero tienden a dar resultados similares. La opción (a) se implementa de la glmsiguiente manera:

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

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

1. Usando weightsargumento ( uno , dos ):

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

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

2. Usando el modelo beta mixto:

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

   o

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

   Si hay ceros exactos o unos en los datos de respuesta, entonces se puede usar el modelo beta cero / uno inflado en glmmTMB.

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

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

4. Contabilización de la sobredispersión en el modelo binomial. Esto usa un truco diferente: agregar 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, ya que se glmer()queja de no enteros py produce estimaciones sin sentido. Una solución que se me ocurrió es usar una constante falsa weights=ky asegurarme de que p*ksiempre sea un número entero. Esto requiere redondeo, ppero al seleccionar kque sea lo suficientemente grande, no debería importar mucho. Los resultados no parecen depender del valor de k.

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

   Actualización posterior (enero de 2018): este podría ser un enfoque no válido. Ver discusión aquí . Tengo que investigar esto más.

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

La opción n. ° 2 es muy lenta y tiene problemas para converger: glmmadmbtarda entre cinco y diez minutos en ejecutarse (¡y aún se queja de que no convergió!), Mientras que lmerfunciona en una fracción de segundo y glmertarda un par de segundos. Actualización: intenté glmmTMBcomo se sugiere en los comentarios de @BenBolker y funciona casi tan rápido como glmersin problemas de convergencia. Entonces esto es lo que usaré.

Las opciones 3 y 4 producen estimaciones muy similares e intervalos de confianza de Wald muy similares (obtenidos con confint). Sin embargo, no soy un gran fanático del # 3 porque es una especie de trampa. Y el # 4 se siente algo hacky.

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