0
$\begingroup$

I have two datasets from different years. They are structured almost identically. One is substantially larger than the other and has more time points. When I prepare a negative binomial generalized linear mixed model for one using glmmTMB and then conduct a multiple comparison (Dunnett's test) using emmeans for each 'level' of time (although time is a continuous predictor, not a factor), the contrast output correctly shows treatment vs control at each 'level' of time. When I repeat the exact procedure using the same code for the second dataset, the output contrast is not at each 'level' of time, but instead one single output is generated at the average time. Why is this happening...?

Code for the dataset that works:


## Create the model
model2.08nymph <- glmmTMB(count ~ treatment * time + (1 + time|block), family=nbinom2, data = master08.nymph)
summary(model2.08nymph)

 Family: nbinom2  ( log )
Formula:          count ~ treatment * time + (1 + time | block)
Data: master08.nymph

     AIC      BIC   logLik deviance df.resid 
  1553.3   1606.9   -764.7   1529.3      628 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev. Corr  
 block  (Intercept) 1.967315 1.40261        
        time        0.003557 0.05964  -0.82 
Number of obs: 640, groups:  block, 16

Overdispersion parameter for nbinom2 family (): 1.18 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)   
(Intercept)             1.30927    1.00266   1.306  0.19162   
treatmentflint          1.97536    1.50028   1.317  0.18795   
treatmentpristine       8.85113    3.63800   2.433  0.01498 * 
treatmentsylgard        0.24207    1.71022   0.142  0.88744   
time                   -0.02576    0.04116  -0.626  0.53140   
treatmentflint:time    -0.11199    0.06174  -1.814  0.06970 . 
treatmentpristine:time -0.49252    0.16129  -3.054  0.00226 **
treatmentsylgard:time  -0.06850    0.06989  -0.980  0.32702   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## Check distribution of simulated residuals
resid.sim.08nymph <- simulateResiduals(model2.08nymph)
plot(resid.sim.08nymph) # perfect

Diagnostic plot 1


## Explore significant interactions
emmeans(model2.08nymph, trt.vs.ctrl ~ treatment | time)

$emmeans
time = 22:
 treatment emmean    SE  df lower.CL upper.CL
 control    0.743 0.431 628   -0.104   1.5892
 flint      0.254 0.437 628   -0.605   1.1129
 pristine  -1.242 0.493 628   -2.210  -0.2730
 sylgard   -0.522 0.467 628   -1.440   0.3949

time = 29:
 treatment emmean    SE  df lower.CL upper.CL
 control    0.562 0.519 628   -0.456   1.5808
 flint     -0.710 0.536 628   -1.763   0.3422
 pristine  -4.870 1.170 628   -7.167  -2.5725
 sylgard   -1.182 0.571 628   -2.303  -0.0615

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
time = 22:
 contrast           estimate    SE  df t.ratio p.value
 flint - control      -0.488 0.614 628  -0.795  0.7412
 pristine - control   -1.984 0.655 628  -3.029  0.0073
 sylgard - control    -1.265 0.636 628  -1.990  0.1218

time = 29:
 contrast           estimate    SE  df t.ratio p.value
 flint - control      -1.272 0.745 628  -1.708  0.2158
 pristine - control   -5.432 1.279 628  -4.247  0.0001
 sylgard - control    -1.744 0.770 628  -2.265  0.0644

Results are given on the log (not the response) scale. 
P value adjustment: dunnettx method for 3 tests 

Code for the dataset that is not working



## Create the model
model2.12nymph <- glmmTMB(count ~ treatment * time + (1 + time|block), family=nbinom2, data = master12.nymph)
summary(model2.12nymph)

 Family: nbinom2  ( log )
Formula:          count ~ treatment * time + (1 + time | block)
Data: master12.nymph

     AIC      BIC   logLik deviance df.resid 
  6788.7   6863.7  -3382.3   6764.7     3828 

Random effects:

Conditional model:
 Groups Name        Variance Std.Dev. Corr  
 block  (Intercept) 0.135402 0.36797        
        time        0.005443 0.07377  -0.83 
Number of obs: 3840, groups:  block, 16

Overdispersion parameter for nbinom2 family (): 0.291 

Conditional model:
                       Estimate Std. Error z value Pr(>|z|)    
(Intercept)            -1.70055    0.23203  -7.329 2.32e-13 ***
treatmentflint         -2.32149    0.37643  -6.167 6.95e-10 ***
treatmentpristine      -3.52118    0.49348  -7.135 9.65e-13 ***
treatmentsylgard       -0.03619    0.32829  -0.110 0.912221    
time                    0.26782    0.04241   6.314 2.71e-10 ***
treatmentflint:time     0.21137    0.06215   3.401 0.000672 ***
treatmentpristine:time  0.18701    0.06988   2.676 0.007451 ** 
treatmentsylgard:time   0.01009    0.05999   0.168 0.866394    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

## Check distribution of simulated residuals
resid.sim.12nymph <- simulateResiduals(model2.12nymph)
plot(resid.sim.12nymph) # outliers present

Diagnostic plot 2

## Explore significant interactions
emmeans(model2.12nymph, trt.vs.ctrl ~ treatment | time)

$emmeans
time = 5:
 treatment emmean    SE   df lower.CL upper.CL
 control   -0.361 0.131 3828   -0.618  -0.1054
 flint     -1.626 0.164 3828   -1.947  -1.3049
 pristine  -2.948 0.234 3828   -3.406  -2.4891
 sylgard   -0.347 0.131 3828   -0.604  -0.0908

Results are given on the log (not the response) scale. 
Confidence level used: 0.95 

$contrasts
time = 5:
 contrast           estimate    SE   df t.ratio p.value
 flint - control     -1.2646 0.209 3828  -6.042  <.0001
 pristine - control  -2.5861 0.268 3828  -9.662  <.0001
 sylgard - control    0.0143 0.185 3828   0.077  0.9977

Results are given on the log (not the response) scale. 
P value adjustment: dunnettx method for 3 tests 

The only major difference between these two analyses that I can see is that the second has one more troublesome residuals. Would this lead to the model being created differently or stop emmeans from properly identifying the 'levels' of time? Maybe there's another approach to make these comparisons? (glht from the package multComp requires time to be a factor)

Any and all help appreciated :)

$\endgroup$
2
  • $\begingroup$ Time is inputted as continuous (numeric) in both models $\endgroup$ Mar 14, 2022 at 23:00
  • $\begingroup$ So honestly think time has a linear effect? If so, you post hoc comparisons should make use of that fact, e.g. estimating and comparing slopes. $\endgroup$
    – Russ Lenth
    Mar 14, 2022 at 23:23

1 Answer 1

2
$\begingroup$

I think you have time as a continuous predictor in these models, rather than as a factor. When there are only 2 different times, that makes no difference in terms of predictions, and emmeans by default treats it like a factor.

But when there are more than 2 distinct times, it does make a difference. And emmeans will use the average time. My guess is the second model is wrong and you should refit it with factor(time) instead of time.

Addendum

If it is indeed deemed appropriate to model times as a linear effect, it is easy to obtain the results you want for a selected set of times...

To use all the different values of time:

emmeans(model2.12nymph, trt.vs.ctrl ~ treatment | time,
    cov.reduce = FALSE)

To use a selected set of times, say 3, 5, and 7:

emmeans(model2.12nymph, trt.vs.ctrl ~ treatment | time,
    at = list(time = c(3, 5, 7))

But another style of post hoc comparisons may be more appropriate -- estimating and comparing the slopes of the fitted trend lines:

emtrends(model2.12nymph, trt.vs.ctrl ~ treatment,
    var = "time")

Consider this last option, because the changes with time are linear, so the contrasts at specified times have systematic relationships with one another.

$\endgroup$
3
  • $\begingroup$ Is it not logical to use time as a continuous predictor and to also investigate the interaction of treatment at every level of time? I expect that in many cases the main effect of time won't be meaningful, for example because data was collected long after treatments wore off.. but it feels dirty turning continuous variables into factors $\endgroup$ Mar 15, 2022 at 1:01
  • $\begingroup$ Dirty? I don't think so, when only a small number of times are part of the design. It's fine to model time as a continuous variable, but I suggest diagnostics to make sure the linear fit is adequate. One such diagnostic would be a lack-of-fir test, which (ironically?) is a comparison with the model with time as a factor. I will add to my answer. $\endgroup$
    – Russ Lenth
    Mar 15, 2022 at 16:02
  • $\begingroup$ For some of my datasets, counts were recorded for populations that came and went, so parabolic over time. In these cases, I won't include time as a continuous predictor but as a factor instead. Does it still make sense to do post hoc tests with emmeans (treatment vs control) on these models? If I get main effects (or are they simple effects?) as well as interaction effects of specific timepoints and specific treatments, the interactions give me the same resolution that Dunnett's would, correct? $\endgroup$ Mar 18, 2022 at 1:35

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.