I am conducting Tukey HSDs with compact letter displays on a series of glmms. One example is how the presence/absence of ticks on an animal varies based on several parameters (using lme4 in R):
model <- glmer(ticks.pa ~ Month.factor + Perturb3*Sex + Age2+(1|CurrentChip),
family=binomial, data=temp)
I am then using Tukey HSDs to look post-hoc at differences between specific levels of Month and Perturb3 (both categorical; using multcomp package). For example:
summary(glht(model,mcp(Month.factor="Tukey")))
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glmer(formula = ticks.pa ~ Month.factor + Perturb3 * Sex + Age2 + (1 | CurrentChip), data = temp, family = binomial)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
2 - 1 == 0 16.0812 66.0132 0.244 1.00000
3 - 1 == 0 15.2714 66.0136 0.231 1.00000
5 - 1 == 0 14.6938 66.0132 0.223 1.00000
6 - 1 == 0 14.2554 66.0132 0.216 1.00000
7 - 1 == 0 15.4118 66.0253 0.233 1.00000
8 - 1 == 0 -16.1893 48.4404 -0.334 1.00000
9 - 1 == 0 -18.2058 48.4348 -0.376 1.00000
10 - 1 == 0 -17.4837 48.4350 -0.361 1.00000
11 - 1 == 0 -18.0147 48.4351 -0.372 1.00000
12 - 1 == 0 -1.7331 0.4929 -3.516 0.00728 **
3 - 2 == 0 -0.8098 0.5412 -1.496 0.83237
5 - 2 == 0 -1.3874 0.4859 -2.855 0.06405 .
6 - 2 == 0 -1.8258 0.4832 -3.779 0.00264 **
7 - 2 == 0 -0.6694 1.3943 -0.480 0.99998
8 - 2 == 0 -32.2705 77.9376 -0.414 0.99999
9 - 2 == 0 -34.2869 77.9334 -0.440 0.99999
10 - 2 == 0 -33.5649 77.9337 -0.431 0.99999
11 - 2 == 0 -34.0959 77.9335 -0.437 0.99999
12 - 2 == 0 -17.8143 66.0149 -0.270 1.00000
5 - 3 == 0 -0.5776 0.4841 -1.193 0.95602
6 - 3 == 0 -1.0160 0.4773 -2.129 0.37218
7 - 3 == 0 0.1405 1.3983 0.100 1.00000
8 - 3 == 0 -31.4607 77.9378 -0.404 1.00000
9 - 3 == 0 -33.4771 77.9337 -0.430 0.99999
10 - 3 == 0 -32.7550 77.9340 -0.420 0.99999
11 - 3 == 0 -33.2860 77.9338 -0.427 0.99999
12 - 3 == 0 -17.0044 66.0153 -0.258 1.00000
6 - 5 == 0 -0.4384 0.3786 -1.158 0.96403
7 - 5 == 0 0.7180 1.3258 0.542 0.99993
8 - 5 == 0 -30.8831 77.9375 -0.396 1.00000
9 - 5 == 0 -32.8995 77.9334 -0.422 0.99999
10 - 5 == 0 -32.1775 77.9336 -0.413 0.99999
11 - 5 == 0 -32.7085 77.9335 -0.420 0.99999
12 - 5 == 0 -16.4269 66.0148 -0.249 1.00000
7 - 6 == 0 1.1564 1.3529 0.855 0.99633
8 - 6 == 0 -30.4447 77.9375 -0.391 1.00000
9 - 6 == 0 -32.4611 77.9334 -0.417 0.99999
10 - 6 == 0 -31.7391 77.9337 -0.407 1.00000
11 - 6 == 0 -32.2701 77.9335 -0.414 0.99999
12 - 6 == 0 -15.9885 66.0149 -0.242 1.00000
8 - 7 == 0 -31.6011 77.9478 -0.405 1.00000
9 - 7 == 0 -33.6176 77.9437 -0.431 0.99999
10 - 7 == 0 -32.8955 77.9440 -0.422 0.99999
11 - 7 == 0 -33.4265 77.9438 -0.429 0.99999
12 - 7 == 0 -17.1449 66.0270 -0.260 1.00000
9 - 8 == 0 -2.0164 0.8282 -2.435 0.19535
10 - 8 == 0 -1.2944 0.8088 -1.600 0.76750
11 - 8 == 0 -1.8254 0.9091 -2.008 0.45936
12 - 8 == 0 14.4562 48.4431 0.298 1.00000
10 - 9 == 0 0.7221 0.3530 2.046 0.43147
11 - 9 == 0 0.1911 0.4397 0.435 0.99999
12 - 9 == 0 16.4727 48.4374 0.340 1.00000
11 - 10 == 0 -0.5310 0.4876 -1.089 0.97654
12 - 10 == 0 15.7506 48.4376 0.325 1.00000
12 - 11 == 0 16.2816 48.4377 0.336 1.00000
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)
Warning message: In RET$pfunction("adjusted", ...) : Completion with error > abseps
I understand the general idea of the Tukey HSD, but I am a bit confused as to interpretation of results with respect to the glmm model. How does it account for the other factors in the model (e.g., Sex and Perturb3 here, as well as CurrentChip) in comparing the means for levels of Month? If it is looking strictly at the model coefficients, then the results would only apply to the reference levels of the other variables in the model (sex, age, perturb)? The results are different if I simply use Month.factor + (1|CurrentChip) in the model, so it clearly is not ignoring the other parameters in the model. Is there a better way to assess the differences between months, while accounting for variation in the other variables?
In this instance, only two comparisons are significant: Dec vs Jan (12-1) and Feb vs Jun (6-2). I am then also conducting a compact letter display to visually help discern significantly different means in my plots. Here:
cld(glht(model, linfct=mcp(Month.factor = "Tukey")))
1 2 3 5 6 7 8 9 10 11 12
"cd" "ad" "ac" "ac" "bc" "ac" "ac" "ac" "ac" "ac" "ab"
For some parasites, the output makes sense to me, but some are like this one for ticks. Based on the original Tukey results, we know that there are differences between 12-1 and 6-2. This makes sense in that neither of those pairs share letters (12 is ab vs 1 is cd; 6 is bc vs 2 is ad). However, I don't really understand how there can be non-consecutive letters for a level. (I am assuming this is a statistical interpretation issue, and not some type of coding issue in R…please correct me if I'm wrong and I can redirect this portion of the question to stackoverflow) The mean for 12 is below that of 1, which makes sense that a and b are lower than c and d. 6 is in between, with bc. Still ok. But how is 2 ad? How can it overlap with the higher and lower groups, but not the middle ones?
Please let me know if any more info is needed, and thanks for any insights!
Best Answer
Answer:
You are correct, the letters indicating homogenous subsets appear to look messed up but, in this instance, it is an accurate portrayal of differences detected by pairwise comparisons.
In most cases, compact letter displays (clds) are generated by ordering least squared or marginal means in ascending or descending order and then drawing a series of unbroken lines that unify means for which there are no pairwise differences detected (Gramm et al. 2007). In this case, lines can be converted to letters and vice versa. This is the most common outcome of pairwise assessments. Here is a figure to illustrate:
However, in some instances and especially for comparisons where the standard errors of differences are not homogeneous, drawing these unbroken lines is not possible (Piepho 2000, Gramm et al. 2007). In these cases, letters represent broken lines and result in non-consecutive letters to be generated. Here is a figure to illustrate:
Likely, in part, due to this confusing scenario of non-sequential cld's, the creator of the emmeans package has developed an alternative to compact-letter displays described as “pairwise P-value plot”, which displays all the P-values in pairwise comparisons (see function emmeans::pwpp()). It produces a figure that has the estimated marginal means of the groups ranked from largest to smallest on the y-axis and the p-value of the Tukey post-hoc comparison on the x-axis, allowing for a display of pairwise comparisons that range from most different (towards the left side of the figure) to most similar. Here is an example:
References:
Piepho, H.P., 2000. Multiple treatment comparisons in linear models when the standard error of a difference is not constant. Biometrical Journal 42 (7), 823–835.
Gramm et al. 2007. Algorithms for compact letter displays: comparison and evaluation. Computational Statistics & Data Analysis 52:725-736.
Lenth, R. 2019. Comparisions and contrasts in emmeans (V.1.4.3.1) vignette. URL https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html. Accessed Jan 2020.