0
$\begingroup$

I've conducted a GLMM on some zero-inflated pollinator count data using the glmmadmb() function. Here is the model, in which the counts are in the variable Total, fixed effects include mean temperature (Mean_temp), treatment (Trt), and an offset covariate which is number of floral units (FU), and the random effects include survey visit (Visit) and plot name (PlotCode):

fvs.glmm5 <- glmmadmb(Total ~ Mean_temp + Trt + (1|Visit)+ 
(1|PlotCode)+offset(log(FU)), data= fvs, family = "nbinom2",  
zeroInflation=TRUE)

Here is the summary of the results:

> summary(fvs.glmm5)

Call:
glmmadmb(formula = Total ~ Mean_temp + Trt + (1 | Visit) + (1 | 
PlotCode) + offset(log(FU)), data = fvs, family = "nbinom2", 
zeroInflation = TRUE)

AIC: 643.5 

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -5.4376     0.6622   -8.21  < 2e-16 ***
Mean_temp     0.0991     0.0244    4.07  4.7e-05 ***
TrtF          0.4249     0.2680    1.59   0.1129    
TrtM          0.3914     0.2612    1.50   0.1341    
TrtX          0.7958     0.2533    3.14   0.0017 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Number of observations: total=260, Visit=7, PlotCode=40 
Random effect variance(s):
Group=Visit
             Variance StdDev
(Intercept)   0.7657  0.875
Group=PlotCode
             Variance    StdDev
(Intercept) 1.126e-07 0.0003356

Negative binomial dispersion parameter: 2.5862 (std. err.: 0.79003)
Zero-inflation: 1.0247e-06  (std. err.:  0.00019615 )

Log-likelihood: -312.757 
Warning message:
In .local(x, sigma, ...) :
  'sigma' and 'rdig' arguments are present for compatibility only:   
ignored

The question I have is how to make sense of the results of the post-hoc tests. When I run the following code, I get results which suggest that treatment C is significantly different from treatment X, but that there are no other significant differences between treatments:

#post-hoc tests to see where differences lie
>fvs.ls = lsmeans(fvs.glmm5,
                     pairwise ~ Trt,
                     adjust="tukey",
                     type="response")         
>fvs.ls
$lsmeans
 Trt  response        SE df  asymp.LCL asymp.UCL
 C   0.5008787 0.4152155 NA 0.09865259  2.543061
 F   0.7660764 0.6596800 NA 0.14167383  4.142424
 M   0.7408063 0.6431391 NA 0.13512157  4.061483
 X   1.1100113 0.9612614 NA 0.20332621  6.059844

Confidence level used: 0.95 
Intervals are back-transformed from the log scale 

$contrasts
 contrast response.ratio        SE df z.ratio p.value
 C - F         0.6538234 0.1752378 NA  -1.585  0.3869
 C - M         0.6761265 0.1766245 NA  -1.498  0.4385
 C - X         0.4512375 0.1143120 NA  -3.141  0.0091
 F - M         1.0341117 0.3965509 NA   0.087  0.9998
 F - X         0.6901519 0.2453648 NA  -1.043  0.7240
 M - X         0.6673862 0.2421148 NA  -1.115  0.6805

P value adjustment: tukey method for comparing a family of 4 estimates 
Tests are performed on the log scale 

The treatment pairwise comparisons seem to make sense until you look at the values for the predicted least-square means and the confidence intervals around those values. If treatment C is truly significantly different from treatment X, then why are the confidence intervals around the estimated means for these two treatments so heavily overlapping?

Any help with making sense of these seemingly contradictory results would be much appreciated!

$\endgroup$
1
$\begingroup$

Because it is a mixed model, I.e., there is more than one source of variation. Apparently, treatment is a within-Visit and/or within-plotCode factor. Accordingly, the variance of the LSMs includes between-visits and between-plotcodes variations, while (all or some of) those variations cancel out when comparing them.

$\endgroup$
  • $\begingroup$ Thank you very much for your explanation. That does help to clarify the situation somewhat. Treatment is a within-plotCode factor, assuming I understand you properly. Each plot contains only one of the four treatments, but there are 10 replicates (plots) per treatment. Therefore, some of the variation in plotCodes is almost definitely due to variation in the treatments. Does this imply that I've set up my GLMM formula incorrectly? Is there any way that I can get confidence intervals for my treatment means that will reflect the significant difference between treatments C and X? $\endgroup$ – D.Hodgkiss Oct 11 '17 at 21:33
  • $\begingroup$ Nope that’s backwards. First, we are worried about RANDOM variations (due to visits an plots per your model). You are saying that treatments are between-plots. It’d be within-plot if each plot had every treatment, so that we could compare them on the same plot. But perhaps you’re saying that each visit has all treatments? That’s the one with the big variance. $\endgroup$ – rvl Oct 11 '17 at 23:16
  • $\begingroup$ Thanks for your explanation. Each visit has all treatments, yes, so I guess treatment is a within-visit factor. Does this mean that I should define the random effects differently when I run the GLMM formula? My goal is to be able to produce a graph with the ls means for each treatment and associated standard errors, but at the moment those SEs are so large that all treatments overlap. $\endgroup$ – D.Hodgkiss Oct 12 '17 at 7:04
  • $\begingroup$ It is WRONG to use side by side confidence intervals to do comparisons. Don’t do that, whether it looks good or not. It’s also wrong to change a valid statistical model into an invalid one in order to get the results you want. $\endgroup$ – rvl Oct 12 '17 at 13:38
  • $\begingroup$ You might try plot(lsmeans(...), comparisons = TRUE). This adds comparison arrows to the plot that are valid for use in judging comparisons. $\endgroup$ – rvl Oct 12 '17 at 13:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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