# Interpreting compact letter displays and Tukey multiple comparisons from glmm

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?

• As a general note, as far as I understand, the "Tukey" in mcp has nothing to do with Tukey HSD; it is simply code to tell the function to make all pairwise comparisons. You will see at the bottom of your output that the summary function uses a single step adjustment by default... As some general advice for your problem, 1) I would confirm that multcomp handles glmer obects correctly. I don't know where there is a list of models it supports. 2) That warning message in the output may be important. I don't know. (Cont...) – Sal Mangiafico Feb 14 '18 at 13:59
• 3) It's curious to me that for the significant contrasts, the estimate of differences is relatively small, but that the standard error is also small. Might there be something unbalanced in your data that's driving the smaller standard errors and therefore the results? 4) You might try the emmeans package and see how that handles your model and data. See ?emmeans::models ... As to how multcomp accounts for other effects in the model, the [documentation ](rdocumentation.org/packages/multcomp/versions/1.4-7/topics/glht) gives some details, but mostly refers you to other texts. – Sal Mangiafico Feb 14 '18 at 14:00
• I strongly urge you to not use CLDs. They highlight non-findings rather than findings. In significance testing, findings are those cases where the P value is small. But the groupings point to cases where the P value is large — where you CAN’T make a definitive conclusion. – Russ Lenth Jan 20 '20 at 14:51

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.