I am having trouble with analysing some of my data. I'm trying to test the effect of a treatment with two levels (Treatment/Control) on the abundance of ants belonging to different dominance categories on plants (Response = Abundance, Predictors = Site, Treatment, Dominance).
There are two sites (SiteC/SiteT) each containing 3 sets of paired plots )((3*treatment + 3*control)*2 Sites). There are 3 categories of ant dominance (D, S or O). Each plot contains 9 plants. Sample provided below.
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
The initial model was: Average_Abundance ~ Treatment * Site * Dominance + (1|Pair)
The issues are that:
Sites responded differently to the treatment and there are interactions.
Very low replicates
Data is zero-inflated, overdispersed and are counts
Low number of levels in categorical variables (Site = 2, Treatment = 2 and Dominance = 3)
One solution suggested to me was to use simple statistics, specifically paired t-tests to analyse this data. However the non-parametric equivalent Wilcoxon signed-rank tests cannot return significance for n=3.
Some solutions I am trying to use:
Instead of using the average abundance per plot, I have used plant as the replicate and nested plot in plot pairs to account for spatial autocorrelation.
use a ZINB GLMM (glmmTMB) instead of glmer.nb (lme4) to address zero-inflation and overdispersion.
However! This brings me to my problem and my questions.
The full model below will not run when ziformula=~. (when the zero-inflated model is the same as the conditional model).
zinb_p1 <- glmmTMB(Abundance ~ Site * Treatment * Dominance + (1|Pair/Site_Plot), data=ant_stats, ziformula=~., family=nbinom1)
I get the following warning message
Warning message: In fitTMB(TMBStruc) : Model convergence problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')
If i change ziformula=~. to ziformula=~1 i get the following output, which I assume does not have a zero-inflated model fitted.
> 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
Is this model overparametised? Is that why it won't work?
If I want to continue to use plant as the replicate rather than plot but nesting plot in pair of plots causes working models to have convergence errors, what alternative is there to account for spatial autocorrelation?
Is there a different approach I could be using that would allow me to keep all of my fixed variables?