Ajuste del modelo de regresión binomial en R: fórmula correcta, prueba de significación y sobredispersión

Question

Estoy utilizando modelos lineales generalizados para probar el efecto de varios predictores en algunos datos binomiales. Mi respuesta es un vector binomial de éxitos y fracasos.

Quiero comprobar si mi predictor categórico de interés (P1) es un predictor significativo de mi respuesta. Sin embargo, hay otra variable categórica (P2) que quiero tener en cuenta. Por lo tanto, la incluí como efecto aleatorio en mi modelo, que no parecía tener mucho impacto sobre si P1 era significativo o no:

mod1<-glmer(Response~P1+(1|P2),family=binomial)
mod0<-glmer(Response~1+(1|P2),family=binomial)
anova(mod0,mod1,test="Chisq")

Sin embargo, los revisores me han pedido ahora más información sobre el efecto de P2: creen que será importante y quieren saber si es significativo o no. Si me limito a comprobar el efecto de P2 como efecto fijo por sí solo, sería significativo, pero creo que esto se debe únicamente a que existe una correlación con P1, no a que sea importante. Por lo tanto, estaba pensando que podría probar si P2 es un predictor significativo dentro de cada nivel de P1, y trató de hacer esto así:

mod2<-glmer(Response~P2+(1|P1),family=binomial)
mod3<-glmer(Response~1+(1|P1),family=binomial)
anova(mod2,mod3,test="Chisq")

Sin embargo, me preocupa que las fórmulas que estoy utilizando sean incorrectas, ya que he leído esta guía sobre las fórmulas de modelos mixtos en R:

"Random effects are specified as e|g, where e is an effect and g is a grouping factor"

Entonces, ¿debería estar haciendo esto?

mod4<-glmer(Response~P2|P1,family=binomial)

Y si es así, ¿cómo compruebo la significación?

Como pregunta adicional, mis datos están muy sobredispersados, pero noto que no puedo usar quasibinomial en lme4. He leído sobre la inclusión de un efecto aleatorio a nivel de observación. Utilizando mi ejemplo anterior, ¿es ésta la forma correcta de hacerlo?

mod5<-glmer(Response~P2+(1|P1)+(1|obs),family=binomial)
mod6<-glmer(Response~1+(1|P1)+(1|obs),family=binomial)
anova(mod5,mod6,test="Chisq")

Esto no parece tener ningún efecto obvio en los resultados, así que no estaba seguro de si lo estaba haciendo bien.

Muchas gracias

Answer 1

Aquí se plantean diversas cuestiones:

• ojalá P2 tiene suficientes niveles (por ejemplo, más de 5 o 6) como para que sea razonable tratarlo como un efecto aleatorio; de lo contrario, es bastante probable que la variabilidad asociada a P2 se pondrá a cero (basta con comprobar VarCorr(mod1) ).

• si no, probablemente tenga sentido tratar P2 como fijo, es decir, utilizar

  mod2 <- glm(Response~P1+P2,family=binomial) drop1(mod2,test="Chisq")

o

mod1 <- glm(Response~P1,family=binomial)
anova(mod2,mod1,test="Chisq")

Usted podría utilice Response~P2|P1 pero que en realidad busca un interacción - diferencias en el efecto de P2 a diferentes niveles de P1 - que no es lo que dijiste que querías. (Además, el efecto aleatorio está centrado en cero, por lo que probablemente debería añadir un efecto principal (fijo) de P2 al modelo).

También podría comprobar si existe un efecto significativo de P2 considerado como un efecto aleatorio comparando glm(Response~P1,...) a glmer(Response~P1+(1|P2),...) pero hay algunas cuestiones sobre las pruebas de cociente de verosimilitud de efectos aleatorios (por ejemplo, véase http://glmm.wikidot.com/faq ); la prueba ingenua será conservadora.

Tenga en cuenta que las "variables molestas" (variables que no le interesan directamente) y los "efectos aleatorios" son no sinónimos; véase, por ejemplo esta respuesta . Yo no sugeriría cambiar la representación de una variable de fija a aleatoria o viceversa dependiendo de las inferencias concretas que intentes hacer.

Su enfoque para tratar la sobredispersión parece razonable, aunque véase este artículo de PeerJ para tener algo de perspectiva... Si la varianza a nivel de observación se estimara como cero (véanse de nuevo los resultados de VarCorr(model) ), eso explicaría por qué los resultados no cambiaron...

Answer 2

Hubo una cosa en tu descripción de lo que querías hacer que captó mi interés, escribiste:

Por lo tanto, estaba pensando que podría probar si P2 es un significativo dentro de cada nivel de P1,

Esa idea suena mucho a que P1 y P2 son realmente efectos fijos. Tal vez estoy sobreinterpretando su elección de términos.

De todos modos, si los niveles de P1 y P2 son interesantes en sí mismos, y no sólo le interesa controlar la molestia de P2, entonces le sugiero, como ya hizo Ben Bolker, que los trate a ambos como efectos fijos. A partir del glmm-faq :

2. Los efectos son fijos si son interesantes en sí mismos o aleatorios si hay interés en la población subyacente. Searle, Casella y McCulloch [(1992), Sección 1.4] exploran esta distinción en profundidad. 3. "Cuando una muestra agota la población, la variable correspondiente es fija; cuando la muestra es una parte pequeña (es decir, insignificante) de la población, la variable correspondiente es aleatoria" [Green y Tukey (1960)].

Si opta por esta vía, tal vez también resulte esclarecedor incorporar un término de interacción en su modelo.

fm1 <- glm(Response ~ P1 * P2, family = binomial)

O con alguna familia que explique la sobredispersión.