4
$\begingroup$

A little bit about my data:

I have four treatment groups: control, early, late, both. For each group, I counted nymphs and eggs on leaves on five different dates. The design is randomized complete block design. A treatment was applied at different times for each treatment group (vertical red dashed lines), which is suspected to reduce egg and nymph counts. For the early group, the treatment was applied on an early date. For the late group, the treatment was applied slightly later than the early group. For the both group, the early and late treatments were both applied. I am analyzing nymphs and eggs separately. I am interested in evaluating the relative effect of treatment timing for each morphology.

Here's what the data looks like: enter image description here

For my analysis, I created the following models:

A negative-binomial hurdle model

model1 <- hurdle(count ~ treatment*date+block, data = masterdata.egg, dist = "negbin", zero.dist = "binomial")

... which, compared to Poisson, fit pretty well: enter image description here

I also tried creating multiple generalized linear mixed models, using different methods:

# Hurdle zero-inflated negative binomial
model2 <- glmmTMB(count ~ treatment*date + (1|block),
               zi=~treatment*date,
               family=truncated_nbinom1, data=masterdata.egg)

# Zero-inflated negative binomial
model3 <- glmmTMB(count ~ treatment*date + (1|block),
                  zi=~treatment*date,
                  family=nbinom1, data=masterdata.egg)

# Negative binomial
model4 <- glmmTMB(count ~ treatment*date + (1|block),
                  family=nbinom1, data=masterdata.egg)

# Hurdle zero-inflated Poisson
model5 <- glmmTMB(count ~ treatment*date + (1|block),
                  zi=~treatment*date,
                  family=truncated_poisson, data=masterdata.egg)

# Zero-inflated Poisson
model6 <- glmmTMB(count ~ treatment*date + (1|block),
                  zi=~treatment*date,
                  family=poisson, data=masterdata.egg)

# Poisson
model7 <- glmmTMB(count ~ treatment*date + (1|block),
                  family=poisson, data=masterdata.egg)

I then compared everything by AIC:

> AIC(model1, model2, model3, model4, model5, model6, model7)
       df      AIC
model1 55 13294.08
model2 50 13431.39
model3 50 13405.21
model4 26 13603.39
model5 49 36352.85
model6 49 36352.71
model7 25 52937.56

... and concluded that the first hurdle model is the most parsimonious!

So then I call the summary...

> summary(model1)

Call:
hurdle(formula = count ~ treatment * date + block, data = masterdata.egg, 
    dist = "negbin", zero.dist = "binomial")

Pearson residuals:
    Min      1Q  Median      3Q     Max 
-1.0051 -0.6738 -0.3888  0.3214  9.7273 

Count model coefficients (truncated negbin with log link):
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    3.24041    0.12546  25.829  < 2e-16 ***
treatmentearly                 0.35026    0.16040   2.184 0.028989 *  
treatmentlate                 -0.08835    0.16089  -0.549 0.582904    
treatmentboth                  0.02973    0.16387   0.181 0.856030    
date2013-05-30                 0.50692    0.16837   3.011 0.002606 ** 
date2013-06-11                 0.50022    0.16564   3.020 0.002529 ** 
date2013-06-24                 0.14098    0.16017   0.880 0.378755    
date2013-07-09                -1.06334    0.19508  -5.451 5.02e-08 ***
date2013-07-22                -0.68781    0.17798  -3.865 0.000111 ***
block2                         0.17429    0.07950   2.192 0.028347 *  
block3                         0.26878    0.08143   3.301 0.000965 ***
block4                         0.34784    0.07954   4.373 1.22e-05 ***
treatmentearly:date2013-05-30 -1.05361    0.24080  -4.375 1.21e-05 ***
treatmentlate:date2013-05-30  -0.24402    0.24119  -1.012 0.311667    
treatmentboth:date2013-05-30  -1.05727    0.25195  -4.196 2.71e-05 ***
treatmentearly:date2013-06-11 -1.44106    0.23749  -6.068 1.30e-09 ***
treatmentlate:date2013-06-11  -0.70219    0.23927  -2.935 0.003339 ** 
treatmentboth:date2013-06-11  -1.14215    0.24599  -4.643 3.43e-06 ***
treatmentearly:date2013-06-24 -1.04128    0.23017  -4.524 6.07e-06 ***
treatmentlate:date2013-06-24  -0.73985    0.23587  -3.137 0.001708 ** 
treatmentboth:date2013-06-24  -1.35066    0.24364  -5.544 2.96e-08 ***
treatmentearly:date2013-07-09 -0.66616    0.26680  -2.497 0.012530 *  
treatmentlate:date2013-07-09  -0.55470    0.27691  -2.003 0.045162 *  
treatmentboth:date2013-07-09  -0.53726    0.29771  -1.805 0.071133 .  
treatmentearly:date2013-07-22  0.06485    0.24165   0.268 0.788437    
treatmentlate:date2013-07-22   0.63994    0.25004   2.559 0.010486 *  
treatmentboth:date2013-07-22   0.45122    0.25189   1.791 0.073234 .  
Log(theta)                     0.12591    0.04787   2.630 0.008534 ** 
Zero hurdle model coefficients (binomial with logit link):
                              Estimate Std. Error z value Pr(>|z|)    
(Intercept)                    1.63126    0.38289   4.260 2.04e-05 ***
treatmentearly                 0.16880    0.54765   0.308 0.757914    
treatmentlate                  0.16880    0.54765   0.308 0.757914    
treatmentboth                  0.12048    0.54865   0.220 0.826190    
date2013-05-30                -1.10252    0.45820  -2.406 0.016118 *  
date2013-06-11                -0.80725    0.47055  -1.716 0.086246 .  
date2013-06-24                 0.16880    0.54765   0.308 0.757914    
date2013-07-09                -2.23588    0.43964  -5.086 3.66e-07 ***
date2013-07-22                -1.59664    0.44521  -3.586 0.000335 ***
block2                         0.72200    0.14711   4.908 9.20e-07 ***
block3                         0.72018    0.14715   4.894 9.87e-07 ***
block4                         0.96257    0.15175   6.343 2.25e-10 ***
treatmentearly:date2013-05-30 -0.54751    0.65364  -0.838 0.402234    
treatmentlate:date2013-05-30  -0.60573    0.65281  -0.928 0.353468    
treatmentboth:date2013-05-30  -1.04664    0.64948  -1.612 0.107068    
treatmentearly:date2013-06-11 -0.66068    0.66546  -0.993 0.320802    
treatmentlate:date2013-06-11  -0.84278    0.66240  -1.272 0.203260    
treatmentboth:date2013-06-11  -1.07630    0.65997  -1.631 0.102924    
treatmentearly:date2013-06-24 -1.14485    0.73388  -1.560 0.118763    
treatmentlate:date2013-06-24  -1.69881    0.72107  -2.356 0.018475 *  
treatmentboth:date2013-06-24  -1.99763    0.71757  -2.784 0.005371 ** 
treatmentearly:date2013-07-09  0.30400    0.63713   0.477 0.633264    
treatmentlate:date2013-07-09   0.02745    0.63693   0.043 0.965620    
treatmentboth:date2013-07-09  -0.70185    0.63892  -1.098 0.271989    
treatmentearly:date2013-07-22  0.77513    0.66746   1.161 0.245516    
treatmentlate:date2013-07-22  -0.11161    0.64367  -0.173 0.862339    
treatmentboth:date2013-07-22  -0.12048    0.64376  -0.187 0.851547    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Theta: count = 1.1342
Number of iterations in BFGS optimization: 36 
Log-likelihood: -6592 on 55 Df

How do I interpret this?

I would assume that

Count model coefficients (truncated negbin with log link):
                              Estimate Std. Error z value Pr(>|z|)    
date2013-06-11                 0.50022    0.16564   3.020 0.002529 ** 

could be interpreted that there are significant differences between my treatments on this date.

But what does something like

Count model coefficients (truncated negbin with log link):
                              Estimate Std. Error z value Pr(>|z|)  
treatmentearly:date2013-06-11 -1.44106    0.23749  -6.068 1.30e-09 ***

tell me? The early treatment, on 2013-06-11, is significant in what way? What is it compared to?

I've never dealt with such a complex summary readout before, so I don't even know where to start here.

$\endgroup$
1
  • $\begingroup$ "I counted nymphs and eggs on leaves on five different dates" You are sure you have 5 dates ? What does levels(masterdata.egg$date) produce ? Also, do you have specific interesy in the date variable, or is it just a replication to improve precision ? $\endgroup$ Commented Aug 10, 2020 at 7:19

1 Answer 1

5
$\begingroup$

There are a lot of estimates in this model, mainly because you have quite a few dates and you are coding the date variable as categorical. If you coded it as numeric there would be much less output. However, looking at the main effects for the date variable there doesn't appear to be a linear trend or any systematic trend, so unless date is an actual variable you have inerest in the estimates for (which doesn't appear to be your research question) then you might consider fitting random intercepts for it instead. 5 is rather few to treat it as random, but it will make the model much more easily intepretable. You could look at model with and without random intercepts for date and compare them. Hopefully they will give you similar insight.

It is common for models with large amounts of output to be a bit intimidating, however there are simple rules that apply to all models, so it's just a question of applying them methodically.

First, the intercept is the estimate of the "outcome" when all other variables are zero, or if they are categorical, then when they are at their reference level.

All of the main effects for categorical variables, such as date013-06-11 that you mention, or treatmentearly are interpreted as contrasts with the reference level for that variable - the reference level for treatment seems to be control. For date it is unclear since you say there are 5 dates so there should only be 4 estimates. I suspect you actually have 6 dates.

However, when a variable is involved in an interaction, the main effect is conditional on the other variable that it is interacted with being at it's reference level. So date2013-06-11 is the estimated difference between the "outcome" at date 2013-06-11 and whatever the reference level for date is, when treatment is at it's reference level - ie, in the control group.

The interactions then tell you the difference between the "outcome" at the relevant date and the reference level date, for the other treatment group. For example treatmentearly:date2013-06-11 is the estimated difference between the "outcome" at 2013-06-11 and whateve the reference date is, for the early treatment group compared to the control group.

I put "outcome" in quotes above because, since this is a hurdle model, each estimate has to be interpreted in terms of which part it belongs to (the 0/not0 or the negative bimomial part). For the 0/not zero part (the 2nd section of output) the estimates are on the log-odds scale. These should be exponentiated and will tell you (for categorical) variable the odds ratios. For the negative binomial part (1st part of the output), it is on the log scale so exponentiating this will tell you the expected change in counts.

$\endgroup$
6
  • $\begingroup$ Thank you Robert for this excellent response. You are correct that I do in fact have six sample dates, not five. To address some of your points, I coded the date as character as using them as class "date" caused convergence issues in the glmmTMB and produced NaNs in the hurdle model. Perhaps these models don't work well with continuous date inputs? $\endgroup$ Commented Aug 10, 2020 at 16:11
  • $\begingroup$ You asked in a previous comment "Also, do you have specific interested in the date variable, or is it just a replication to improve precision?" -- I think date is very important in my model since for my different treatments, the measurement dates directly relate to the treatment application dates. But with that being said, using the first of the six dates as a baseline doesn't seem like the most intuitive method of analysis for my data. Using control as baseline, at every sample date, would make more sense to me $\endgroup$ Commented Aug 10, 2020 at 16:17
  • $\begingroup$ No problem, you are welcome. Rather than using type "date", try it with type numeric with the first date as 0, and the others the number of days after 0. $\endgroup$ Commented Aug 10, 2020 at 16:19
  • $\begingroup$ I think the default in R is that dates are treated as the number of days after 1 Jan 1970, so those will be big numbers which might be the reason you had convergence problems. $\endgroup$ Commented Aug 10, 2020 at 16:26
  • $\begingroup$ Your suggestion of working with day 0 and days after 0 worked great. The output is much simpler. I have a couple more questions if you don't mind... in my new output, I get treatmentearly:date (p = 0.009) but for treatmentlate:date and treatmentboth:date I get non-significance ( p ~= 0.25). Do these p values represent the overall trend of each treatment across the dates, relative to to control:date? $\endgroup$ Commented Aug 10, 2020 at 16:41

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.