# Trouble modeling zero-inflated data. Estimates and standard errors are off with GLM, GLMM, and ZI models

I conducted a study looking at the attraction of different species of insects to 5 different chemical treatments (I have had other issues with this dataset explored here and here). This experiment was conducted over 5 experimental periods (i.e., 5 intervals of 2 weeks). Within each experimental period, the experiment was replicated 4 times. Thus, my dataset is 100 rows (5 treatments x 5 dates x 4 replicates within each date; one replicate was removed due to human error, leaving me with 95 rows) with ~50 columns detailing the capture of insects in my study.

For a reproducible example, one of my columns looks like:

       replicate treatment  date insect.species1
1        rep1   treatment1  date1       0
2         rep1  treatment2  date1       0
3         rep1  treatment3  date1       0
4         rep1  treatment4  date1       0
5         rep1  treatment5  date1       0
6         rep2  treatment1  date1       0
7         rep2  treatment2  date1       0
8         rep2  treatment3  date1       0
9         rep2  treatment4  date1       0
10        rep2  treatment5  date1       4
11        rep3  treatment1  date1       0
12        rep3  treatment2  date1       0
13        rep3  treatment3  date1       0
14        rep3  treatment4  date1       0
15        rep3  treatment5  date1       6
16        rep4  treatment1  date1       0
17        rep4  treatment2  date1       0
18        rep4  treatment3  date1       0
19        rep4  treatment4  date1       0
20        rep4  treatment5  date1       3
21        rep1  treatment1  date2       0
22        rep1  treatment2  date1       0
23        rep1  treatment3  date1       0
24        rep1  treatment4  date1       0
25        rep1  treatment5  date1       0
26        rep3  treatment1  date2       0
27        rep3  treatment2  date1       0
28        rep3  treatment3  date1       1
29        rep3  treatment4  date1       1
30        rep3  treatment5  date1       3
31        rep4  treatment1  date2       0
32        rep4  treatment2  date1       0
33        rep4  treatment3  date1       0
34        rep4  treatment4  date1       0
35        rep4  treatment5  date2       2
36        rep1  treatment1  date3       0
37        rep1  treatment2  date3       0
38        rep1  treatment3  date3       0
39        rep1  treatment4  date3       0
40        rep1  treatment5  date3       0
41        rep2  treatment1  date3       0
42        rep2  treatment2  date3       0
43        rep2  treatment3  date3       0
44        rep2  treatment4  date3       0
45        rep2  treatment5  date3       0
46        rep3  treatment1  date3       0
47        rep3  treatment2  date3       0
48        rep3  treatment3  date3       1
49        rep3  treatment4  date3       0
50        rep3  treatment5  date3       3
51        rep4  treatment1  date3       0
52        rep4  treatment2  date3       1
53        rep4  treatment3  date3       0
54        rep4  treatment4  date3       0
55        rep4  treatment5  date3       0
56        rep1  treatment1  date4       0
57        rep1  treatment2  date4       0
58        rep1  treatment3  date4       0
59        rep1  treatment4  date4       0
60        rep1  treatment5  date4       0
61        rep2  treatment1  date4       0
62        rep2  treatment2  date4       0
63        rep2  treatment3  date4       0
64        rep2  treatment4  date4       0
65        rep2  treatment5  date4       0
66        rep3  treatment1  date4       0
67        rep3  treatment2  date4       0
68        rep3  treatment3  date4       0
69        rep3  treatment4  date4       0
70        rep3  treatment5  date4       0
71        rep4  treatment1  date4       0
72        rep4  treatment2  date4       0
73        rep4  treatment3  date4       0
74        rep4  treatment4  date4       0
75        rep4  treatment5  date4       0
76        rep1  treatment1  date5       0
77        rep1  treatment2  date5       0
78        rep1  treatment3  date5       0
79        rep1  treatment4  date5       0
80        rep1  treatment5  date5       0
81        rep2  treatment1  date5       0
82        rep2  treatment2  date5       0
83        rep2  treatment3  date5       0
84        rep2  treatment4  date5       0
85        rep2  treatment5  date5       0
86        rep3  treatment1  date5       0
87        rep3  treatment2  date5       0
88        rep3  treatment3  date5       0
89        rep3  treatment4  date5       0
90        rep3  treatment5  date5       0
91        rep4  treatment1  date5       0
92        rep4  treatment2  date5       0
93        rep4  treatment3  date5       0
94        rep4  treatment4  date5       0
95        rep4  treatment5  date5       0


Aggregated, the data look as follows:

aggregate(insect.species1~treatment, mean, data=insectdata, na.rm=TRUE)

treatment insect.species1
1  treatment5 1.10526316
2  treatment3 0.10526316
3  treatment4 0.05263158
4  treatment2 0.05263158
5  treatment1 0.00000000


I have been attempting to analyze the effect of treatment on the number of insects captured with generalized linear and generalized linear mixed models. After a bit of model selection and fussing around with variables, the best model I have found is a negative binomial generalized linear model with the following form:

insectspecies1.nb = glm.nb(insect.species1 ~ treatment + date + replicate), data = insectdata)

summary(insectspecies1.nb)

glm.nb(formula = insect.species1 ~ treatment + date + replicate, data = insectdata, init.theta = 7230.591834, link = log)

Deviance Residuals:
Min        1Q    Median        3Q       Max
-1.36620  -0.01357  -0.00001   0.00000   2.07745

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)           -20.8970  8251.4182  -0.003  0.99798
treatment3             -2.3514     0.7400  -3.177  0.00149 **
treatment4             -3.0445     1.0236  -2.974  0.00294 **
treatment2             -3.0445     1.0236  -2.974  0.00294 **
treatment1            -22.6669 11061.2096  -0.002  0.99836
date1                   0.9555     0.5263   1.815  0.06946 .
date3                 -21.0477 10088.0573  -0.002  0.99834
date2                   0.5878     0.5990   0.981  0.32643
date4                 -21.0477 10088.0573  -0.002  0.99834
rep2                   20.8280  8251.4182   0.003  0.99799
rep3                   21.7444  8251.4182   0.003  0.99790
rep4                   20.8280  8251.4182   0.003  0.99799
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Negative Binomial(7230.592) family taken to be 1)

Null deviance: 121.88  on 94  degrees of freedom
Residual deviance:  19.61  on 83  degrees of freedom
AIC: 72.126

Number of Fisher Scoring iterations: 1

Theta:  7231
Std. Err.:  85985
Warning while fitting theta: iteration limit reached

2 x log-likelihood:  -46.126


One issue here is some of my treatments, dates, and replicates are missing in this output and it's not clear where they went. The main issue though, is that in spite of the above model being the "best" according to likelihood ratio tests and AIC, I am getting some odd estimates and standard errors. Removing date and replicate as factors in the model does not seem to help.

Based upon what I would expect from the data, treatment 5 should be significantly higher than all my other treatments. Performing pairwise comparisons gives me an answer that is close to what I'd anticipate, but is off because of the funky standard errors from treatment 1

cld(glht(h350.nb, mcp(treatment="Tukey")))

treatment5 treatment3 treatment4 treatment2 treatment1
"a"      "b"        "b"        "b"        "ab"


Because of the number of zeroes in my dataset, I recognize that my data are zero-inflated and I should attempt to use a zero-inflated model to see if I get a better fit...Well, that doesn't seem to work either:

insect.formula <- formula(insect.species1 ~ treatment + date + replicate)

Insect.species1.zi = zeroinfl(insect.formula, dist = "negbin",
link = "logit", data = insectdata)

Warning messages:
1: glm.fit: fitted rates numerically 0 occurred
2: glm.fit: fitted probabilities numerically 0 or 1 occurred
> summary(Zip1)

Call:
zeroinfl(formula = f1, data = sticky2016, dist = "negbin", link = "logit")

Pearson residuals:
Min         1Q     Median         3Q        Max
-1.011e+00 -4.293e-05 -1.104e-08 -4.178e-13  3.533e+00

Count model coefficients (negbin with log link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)           -19.9150 13082.6916  -0.002  0.99879
treatment3            -0.9029     0.8093  -1.116  0.26456
treatment4            -1.0856     1.1122  -0.976  0.32905
treatment2            -2.6902     1.0391  -2.589  0.00963 **
treatment1            -21.6671 13207.4652  -0.002  0.99869
date1                  1.1621     0.5500   2.113  0.03462 *
date3                 -20.0479 11188.0528  -0.002  0.99857
date2                  0.4060     0.6305   0.644  0.51969
date5                 -20.0479 11188.0528  -0.002  0.99857
rep2                   19.8510 13082.6916   0.002  0.99879
rep3                   20.5943 13082.6916   0.002  0.99874
rep4                   19.9374 13082.6916   0.002  0.99878
Log(theta)             28.7691     0.2615 110.016  < 2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept)             -11.53         NA      NA       NA
treatment3               52.62     548.92   0.096    0.924
treatment4               68.02     564.34   0.121    0.904
treatment2               22.03     544.62   0.040    0.968
treatment1               22.79 1876248.29   0.000    1.000
date1                    21.26     531.11   0.040    0.968
date3                    21.01 6706847.13   0.000    1.000
date2                   -14.74      67.61  -0.218    0.827
date5                    21.01 6706847.14   0.000    1.000
rep2                    -24.53         NA      NA       NA
rep3                    -49.49         NA      NA       NA
rep4                    -19.38         NA      NA       NA
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Theta = 3120901446336.07
Number of iterations in BFGS optimization: 264
Log-likelihood: -18.39 on 25 Df
Warning message:
In sqrt(diag(object\$vcov)) : NaNs produced


So I'm at a bit of a loss of what to do here. The AIC is smaller for the negative binomial model compared to the zero-inflated model.

I attempted to transform my data by adding a small number (0.01) to see if that influenced the estimates in my model. That model is the only one that seems to produce any sort of logical results with respect to pairwise comparisons and what I'd expect:

#Added 0.01 to each observation for insect.species1
insect.species1.transformed = lmer(insect.species1new ~ treatment + date + (1 | replicate), data = insectdata)

summary(insect.species1.transformed)

inear mixed model fit by REML ['lmerMod']
Formula: insect.species1new ~ treatment + date + (1 | replicate)
Data: insectdata

REML criterion at convergence: 232.1

Scaled residuals:
Min      1Q  Median      3Q     Max
-1.6719 -0.5126 -0.0326  0.2892  5.4240

Random effects:
Groups   Name        Variance Std.Dev.
replicate    (Intercept) 0.03767  0.1941
Residual             0.63196  0.7950
Number of obs: 95, groups:  block, 4

Fixed effects:
Estimate Std. Error t value
(Intercept)           1.1021     0.2600   4.238
treatment3           -1.0000     0.2579  -3.877
treatment4           -1.0526     0.2579  -4.081
treatment2           -1.0526     0.2579  -4.081
treatment1           -1.1053     0.2579  -4.285
date1                 0.4000     0.2514   1.591
date3                -0.2500     0.2514  -0.994
date2                 0.2121     0.2742   0.774
date4                -0.2500     0.2514  -0.994

cld(glht(insect.species1.transformed, mcp(treatment="Tukey")))

treatment5 treatment3 treatment4 treatment2 treatment1
"b"         "a"        "a"        "a"        "a"


As we'd expect though, this model has a high AIC (241.9) compared to our glm.nb model (72.1). So, this is my question/concern.

I'm not really sure where or how to move forward in figuring out what the best model is for this data?

I thought it would be some sort of zero-inflated model (this is my first time using them), but they seems to perform not really that much better than the other models. I don't really understand or know why I'm getting such large estimates for the standard error of treatments like treatment 1, where I captured no insects whatsoever. I presume this is some function of the math and how zeroes interact with one another during the model building process. But, I obviously can't say, keep my negative binomial model and just say, "ignore the ab there", because if that's happening, there is obviously something wrong with my whole model correct? I have tried modeling with fewer factors (i.e., removing date and replicate, they aren't really factors of interest), but the AIC on those models are higher and likelihood ratio tests suggest they are significantly worse than my full models.

I'd appreciate any input or insight on how to move forward and to better understand what is going on here. I also hope Cross-Validated was the correct site (opposed to say, stackoverflow) for this question. If I'm missing any necessary info, please let me know and I'd be happy to fill anything in. Thanks much.

• You have lots of zeroes in your data! Within each date, do you apply all 5 treatments such that you get 4 replicates per treatment within that date? What was the rationale between getting 4 replicates? In particular, would you be able to define your outcome as the average of whatever you are measuring for each replicate (if focused on one species)? How many insect species do you have and what exactly are you measuring for each of them? Would it be possible to analyze the entire "community of species" rather than one species at a time? – Isabella Ghement Mar 26 at 3:17
• Thanks for the comment. The number of zeroes in my data is typical of this sort of trapping experiment. When we conduct these sorts of experiments we often find a moderate to strong response to one or more of our treatments (lots of our target insects), and no response whatsoever to the rest of the treatments. Within each replicate, all 5 treatments are deployed in the field. The purpose of these replicates is that our insects are often very density dependent and if the traps aren't put in the "right" area, we catch none of them. 1/n – Todd Johnson Mar 26 at 19:47
• In the sense of the type of zeroes, those are not true zeroes, because the insects we are targeting weren't in that area, rather than the zeroes being reflective of a response to our treatments. It's typical by many researchers in my field to actually drop replicates and dates where none of our target insects are captured in any treatments within a replicate or date, because it's viewed as they either weren't in that location, they were in the location, but abiotic effects such as temperature caused them to not be responsive to treatments, or 2/n – Todd Johnson Mar 26 at 19:50
• they normally would be in the location, but they are not present in the environment during that period of time (we would refer to this as phenology). I have attempted to analyze the data in this fashion, but by dropping large numbers of replicates and treatments often leaves a very small dataset to work with. The Friedman test is usually used to analyze data like this, and with small datasets, it seems we do not have the power to detect actual responses to treatments, in spite of what the means and our plots of our data show. 3/n – Todd Johnson Mar 26 at 19:53
• With respect to your actual questions, I have tried averaging in various ways, but I'm not sure if I have tried averaging over replicate before my analysis. I will give that a try and report back. As for the number of insects I am looking at, it's about 50 species currently. It would be possible to look at the community response, but those data would have the same structure as these data. Because of how I sampled as well, looking at my data from a community perspective would not be appropriate. 4/4 – Todd Johnson Mar 26 at 19:56

Plotting the data first:

I would agree with comments that say that mixed models are overkill for this data set: in the first four treatments you have a total of four non-zero observations (all ones). You could fit a Poisson or negative binomial to a model of treatment alone: zero-inflation does not seem necessary here; these data are perfectly consistent with Poisson or NB distributions with very low means. However, here you'll run into complete separation/Hauck-Donner problems. So let's use bias-reduced GLM, from the brglm2 package:

summary(m1 <- glm(insect.species~treatment,
data=dd,
family=poisson,
method="brglmFit"))
## ...
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)           -3.638      1.414  -2.572  0.01011 *
## treatmenttreatment2    1.099      1.633   0.673  0.50110
## treatmenttreatment3    1.609      1.549   1.039  0.29886
## treatmenttreatment4    1.099      1.633   0.673  0.50110
## treatmenttreatment5    3.761      1.431   2.629  0.00856 **


This confirms statistically what we think the picture shows; there are no clear differences among treatments except for the last treatment.

Some alternative/follow-up tests or analyses:

You could refit the model with sum-to-zero contrasts:

summary(m2 <- update(m1, contrasts=list(treatment=contr.sum)))


Plot coefficients:

dotwhisker::dwplot(m1)+geom_vline(xintercept=0,lty=2)


Compute and plot expected marginal means:

plot(emmeans::emmeans(m1,~treatment))

• I appreciate you taking the time to look at my data and suggest some paths forward. I was unaware of the brglm2 package and it seems to report some more realistic estimates and standard errors. – Todd Johnson Mar 30 at 20:08