1
$\begingroup$

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.

Improve this question
$\endgroup$
10
  • $\begingroup$ 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? $\endgroup$
    – Isabella Ghement
    Commented Mar 26, 2019 at 3:17
  • $\begingroup$ 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 $\endgroup$
    – Todd D. Johnson
    Commented Mar 26, 2019 at 19:47
  • $\begingroup$ 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 $\endgroup$
    – Todd D. Johnson
    Commented Mar 26, 2019 at 19:50
  • $\begingroup$ 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 $\endgroup$
    – Todd D. Johnson
    Commented Mar 26, 2019 at 19:53
  • $\begingroup$ 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 $\endgroup$
    – Todd D. Johnson
    Commented Mar 26, 2019 at 19:56

1 Answer 1

Reset to default
2
$\begingroup$

Plotting the data first:

enter image description here

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))
Improve this answer
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$
    – Todd D. Johnson
    Commented Mar 30, 2019 at 20:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.