¿Qué enfoque puedo utilizar para analizar datos de recuento inflados a cero, sobredispersados, con réplicas muy bajas y un efecto aleatorio anidado?

Question

Tengo problemas para analizar algunos de mis datos. Estoy intentando probar el efecto de un tratamiento con dos niveles (Tratamiento/Control) sobre la abundancia de hormigas pertenecientes a diferentes categorías de dominancia en las plantas (Respuesta = Abundancia, Predictores = Sitio, Tratamiento, Dominancia).

Hay dos sitios (SitioC/SitioT), cada uno de los cuales contiene 3 conjuntos de parcelas emparejadas ((3*tratamiento + 3*control)*2 Sitios). Hay 3 categorías de dominancia de hormigas (D, S u O). Cada parcela contiene 9 plantas. A continuación se ofrece una muestra.

 Site    Pair Plot Site_Plot Treatment   Check     Plant   Dominance Abundance Presence
1   SiteC SiteC_1   1C  SiteC_1C   Control Initial  Plant_73    Dominant         0        0
2   SiteC SiteC_1   1C  SiteC_1C   Control Initial  Plant_73 Subdominant         0        0
3   SiteC SiteC_1   1C  SiteC_1C   Control Initial  Plant_73 Subordinate         0        0
4   SiteC SiteC_1   1C  SiteC_1C   Control Initial  Plant_76    Dominant         0        0
5   SiteC SiteC_1   1C  SiteC_1C   Control Initial  Plant_76 Subdominant         0        0
6   SiteC SiteC_1   1C  SiteC_1C   Control Initial  Plant_76 Subordinate         1        1
308 SiteT SiteT_3   3T  SiteT_3T Treatment Initial  Plant_50    Dominant         0        0
309 SiteT SiteT_3   3T  SiteT_3T Treatment Initial  Plant_50 Subdominant         0        0
310 SiteT SiteT_3   3T  SiteT_3T Treatment Initial  Plant_50 Subordinate        10        1
311 SiteT SiteT_3   3T  SiteT_3T Treatment Initial  Plant_53    Dominant         0        0
312 SiteT SiteT_3   3T  SiteT_3T Treatment Initial  Plant_53 Subdominant         0        0
313 SiteT SiteT_3   3T  SiteT_3T Treatment Initial  Plant_53 Subordinate         1        1

El modelo inicial era: Cantidad_promedio ~ Tratamiento * Sitio * Dominancia + (1|Pareja)

Los problemas son que:

• Los sitios respondieron de manera diferente al tratamiento y hay interacciones.

• Réplicas muy bajas

• Los datos están inflados a cero, sobredispersados y son recuentos

• Bajo número de niveles en las variables categóricas (Sitio = 2, Tratamiento = 2 y Dominancia = 3)

Una de las soluciones que se me sugirieron fue utilizar estadísticas simples, concretamente pruebas t pareadas, para analizar estos datos. Sin embargo, las pruebas no paramétricas de rango con signo de Wilcoxon no son significativas para n=3.

Algunas soluciones que estoy tratando de utilizar:

• En lugar de utilizar la abundancia media por parcela, he utilizado la planta como réplica y he anidado la parcela en pares de parcelas para tener en cuenta la autocorrelación espacial.

• utilizar un GLMM ZINB (glmmTMB) en lugar de glmer.nb (lme4) para abordar la inflación cero y la sobredispersión.

Sin embargo Esto me lleva a mi problema y a mis preguntas.

El modelo completo a continuación no se ejecutará cuando ziformula=~. (cuando el modelo inflado a cero es el mismo que el modelo condicional).

zinb_p1 <- glmmTMB(Abundance ~ Site * Treatment * Dominance   + (1|Pair/Site_Plot),
                       data=ant_stats,
                       ziformula=~.,
                       family=nbinom1)

Me aparece el siguiente mensaje de advertencia

Warning message:
In fitTMB(TMBStruc) :
  Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')

Si cambio ziformula=~. a ziformula=~1 obtengo la siguiente salida, que supongo que no tiene un modelo inflado a cero ajustado.

 > summary(zinb_p3)
 Family: nbinom1  ( log )
Formula:          Abundance ~ Site * Treatment * Dominance + (1 | Pair/Site_Plot)
Zero inflation:             ~1
Data: ant_stats

     AIC      BIC   logLik deviance df.resid 
   722.2    782.2   -345.1    690.2      297 

Random effects:

Conditional model:
 Groups         Name        Variance  Std.Dev. 
 Site_Plot:Pair (Intercept) 5.902e-11 7.682e-06
 Pair           (Intercept) 5.120e-02 2.263e-01
Number of obs: 313, groups:  Site_Plot:Pair, 12; Pair, 6

Overdispersion parameter for nbinom1 family (): 4.74 

Conditional model:
                                                  Estimate Std. Error z value Pr(>|z|)   
(Intercept)                                     -5.453e-01  4.745e-01  -1.149  0.25045   
SiteTWP                                          1.120e+00  5.559e-01   2.015  0.04394 * 
TreatmentTreatment                              -1.175e-01  6.301e-01  -0.186  0.85209   
DominanceSubdominant                             1.455e-01  6.271e-01   0.232  0.81657   
DominanceSubordinate                             1.412e+00  5.107e-01   2.765  0.00569 **
SiteTWP:TreatmentTreatment                      -5.310e-01  7.825e-01  -0.679  0.49741   
SiteTWP:DominanceSubdominant                    -2.991e+01  6.858e+05   0.000  0.99997   
SiteTWP:DominanceSubordinate                    -2.251e+00  6.934e-01  -3.246  0.00117 **
TreatmentTreatment:DominanceSubdominant         -5.825e-01  9.608e-01  -0.606  0.54436   
TreatmentTreatment:DominanceSubordinate          1.213e-01  7.122e-01   0.170  0.86474   
SiteTWP:TreatmentTreatment:DominanceSubdominant -1.786e+00  3.371e+06   0.000  1.00000   
SiteTWP:TreatmentTreatment:DominanceSubordinate  1.612e+00  9.631e-01   1.674  0.09421 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Zero-inflation model:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -19.09    5788.09  -0.003    0.997

Mis preguntas:

1. ¿Este modelo está sobreparametrizado? ¿Es por eso que no funciona?

2. Si quiero seguir utilizando la planta como réplica en lugar de la parcela, pero la anidación de la parcela en un par de parcelas hace que los modelos de trabajo tengan errores de convergencia, ¿qué alternativa hay para tener en cuenta la autocorrelación espacial?

3. ¿Existe un enfoque diferente que pueda utilizar y que me permita mantener todas mis variables fijas?

Answer 1

Según la Introducción proporcionada por los creadores de glmmTMB, ziformula = ~ 1 significa que se ha aplicado un parámetro de inflación cero a todas las observaciones. ziformula = ~ 0 excluiría un ajuste por inflación cero.