GLMM para datos overdispersed

Question

Estoy analizando los datos de 3 experimentos de campo (granjas=3) para los cítricos de la flor de la enfermedad: la variable de respuesta es binomial porque la flor sólo puede estar enfermo o sano.

Tengo particular interés en la comparación de 5 fungicida equipos de pulverización (trt=5). No, no estoy interesado en el efecto de una granja específica, simplemente representan el total de fincas de la región donde quiero sugerir los mejores tratamientos.

Cada granja tiene 4 bloques (bk=4), incluyendo 2 de los árboles como submuestras (árbol=2) en la que evaluó 100 flores cada uno.

Este es un vistazo rápido de los datos:

dinc <- within(dinc, { tree_id <- as.factor(interaction(farm, trt, bk, tree)) })

farm      trt      bk   tree   dis   tot  tree_id
<fctr>   <fctr>  <fctr> <fctr><int> <int> <fctr>
iaras   Calendar   1      1     0   100  iaras.Calendar.1.1
iaras   Calendar   1      2     1   100  iaras.Calendar.1.2
iaras   Calendar   2      1     1   100  iaras.Calendar.2.1
iaras   Calendar   2      2     3   100  iaras.Calendar.2.2

El modelo que yo consideraba era:

resp <- with(df, cbind(dis, tot-dis)) 

m1 = glmer(resp ~ trt + (1|farm/bk) , family = binomial, data=df) 

He probado la sobredispersión con el overdisp_fun() de GLMM página

        chisq         ratio             p          logp 
 4.191645e+02  3.742540e+00  4.804126e-37 -8.362617e+01 

Como la proporción (residual dev/residual df) > 1, y el valor de p < 0.05, he considerado que añadir la observación de nivel de efectos aleatorios (enlace) para lidiar con la sobredispersión.

así que ahora se añade un efecto aleatorio para cada fila (tree_id) para el modelo, pero no estoy seguro de cómo incluirlo. Este es mi planteamiento:

m2 = glmer(resp ~ trt + (1|farm/bk) + (1|tree_id), family = binomial, data=df)

También me pregunto si farm debe ser un efecto fijo, ya que solo tiene 3 niveles de...

m3 = glmer(resp ~ trt * farm + (1|farm:bk) + (1|tree_id), family = binomial, data=df)

Realmente agradezco sus sugerencias acerca de mi modelo de especificaciones...

Answer 1

Sobredispersión Problema

Parece que estás modelando una variable de recuento como un binomio y creo que ese es el origen de la sobredispersión.

Usted podría modelo de todo como una distribución binomial, pero el total para cada observación es exactamente el mismo. Además, el recuento de plantas enfermas nunca alcanza el máximo de 100, así que no es realmente censurado la forma de un binomio sería.

EDIT: por Lo tanto, usted puede fácilmente este informe como una "tasa" de la enfermedad sobre el total de la muestra. De esta manera se podría analizar el 'conde' de la enfermedad o proporción (enfermedad / total) como un modelo binomial negativo.

EDIT2: Ya que parece que hay cierta incomodidad para el uso de una binomial negativa, aquí hay una lista de los últimos fitopatología de artículos (el mismo de la disciplina como OP) que el modelo de la enfermedad como una binomial negativa (Prager et al., 2014, Mori et al., 2008, Passey et al., 2017, Paiva de Almeida et al., 2016)

Un histograma de la y variable parece un cero inflado binomial negativa.

Tenga en cuenta la larga cola derecha que se suele ver con una negativa binomial o de Poisson.

Hay un par de maneras diferentes para controlar esto, pero hay una solución fácil:

m4<-glmer.nb(dis ~ trt + (1 | farm/bk),data = dinc)

summary(m4)
overdisp_fun(m4)

Tengo las siguientes sobredispersión resultados:

      chisq       ratio         rdf           p 
122.1655582   1.0811111 113.0000000   0.2617332 

Se ve bien, ¿verdad?

(EDIT: Ignorar tachado parte de abajo)

### Lado de la Cuestión: los Árboles son Observaciones Independientes

En primer lugar, se ve como cada uno de los dos árboles se debe a un efecto aleatorio.

Sin embargo, Árbol 1 en la granja 1 no es comparable a la de Árbol 1 en la granja 2. Por lo tanto, usted no desea modelar el efecto del Árbol como un efecto aleatorio. Imagínese si cada Árbol era una persona diferente. La adición de un efecto aleatorio para cada persona no importa a menos que usted tenía múltiples observaciones por persona.

De manera similar, incluyendo la granja de "bloque" en realidad no tiene un efecto en el modelo.

---

Modelos alternativos y reflexiones Finales

• Potencialmente podría retirar cero inflado binomial negativa
• A pesar de su dispersión no parece mala con la norma nb
• La MASA paquete es una forma alternativa de ejecutar un nb modelo
• Además, puede ejecutar esto como un Cuasi-Poisson
• Voy a incluir un poco de código a continuación, en caso de querer seguir este

---

 require("MASS")

 m5<-glmmPQL(dis ~ trt ,
             random = ~ 1 | farm/bk,
             family = negative.binomial(theta=9.86), 
             data = dinc)

 summary(m5)

 m6<-glmmPQL(dis ~ trt ,
             random = ~ 1 | farm/bk,
             family = quasipoisson(link='log'), 
             data = dinc)

 summary(m6)

La mejor de las suertes con tu modelo!

---

EDITAR En caso de que desee ejecutar esto como una "tasa", por favor pruebe este código:

dinc$dis_prob<- dinc$dis / dinc$tot 

m7<-glmmPQL(dis_prob ~ trt ,
             random = ~ 1 | farm/bk,
             family = quasipoisson(link='log'), 
             data = dinc)

summary(m7)

Answer 2

usted está en lo correcto que se han modelado de forma adecuada. Cada uno de sus flores es "anidada" debajo de un árbol y por lo tanto no son independientes el uno del otro. Su código es apropiado cuando usted ha permitido que la intersección para variar de los árboles.

También parece que han examinado la correlación intraclase (es decir, el overdisp_fun() que se utiliza).

Además, puesto que la finca cuenta con tres niveles, es apropiado sólo lo tratan como fijo (especialmente si usted no se preocupan realmente de la diferencia). En este caso, la prueba de la inclusión de los niveles fijos, y si no mejorar el ajuste, a continuación, usted puede deshacerse de ellos.

Asegúrese de que el examen de la AIC y BIC para ayudar con la construcción del modelo.