Any other options besides the traditional CLD bar graph?

Question

I am looking for an alternative approach to plotting results from pairwise comparisons than traditional bar plots. If possible, I would like to create a plot like the one shown below [1], but for a model that includes an interaction effect. R code for the plot below is online [2]. Is there a way to revise or add onto this code to include an interaction effect?

Example of my data set (too large to include in its entirety but I can send upon request) and the model used:

aq <- tibble::tribble(
  ~trt, ~season, ~site,         ~sp,
  "herbicide", "early",    1L, 0.120494496,
  "herbicide", "early",    1L,  0.04057757,
  "herbicide", "early",    1L, 0.060556802,
  "herbicide", "early",    1L, 0.050567186,
  "herbicide", "early",    1L, 0.110504881,
  "herbicide", "early",    1L, 0.090525649,
  "herbicide", "early",    1L, 0.100515265,
  "herbicide", "early",    1L, 0.030587954,
  "herbicide", "early",    1L, 0.080536033,
  "herbicide", "early",    1L, 0.010608723,
  "herbicide", "early",    1L, 0.080536033,
  "herbicide", "early",    1L,  0.04057757,
  "herbicide",   "mid",    1L, 0.050567186,
  "herbicide",   "mid",    1L, 0.050567186,
  "herbicide",   "mid",    1L,  0.04057757,
  "herbicide",   "mid",    1L,  0.04057757,
  "herbicide",   "mid",    1L, 0.140473728,
  "herbicide",   "mid",    1L, 0.030587954,
  "herbicide",   "mid",    1L, 0.150463344,
  "herbicide",   "mid",    1L, 0.020598339,
  "herbicide",   "mid",    1L, 0.120494496,
  "herbicide",   "mid",    1L,  0.04057757,
  "herbicide",   "mid",    1L, 0.050567186,
  "herbicide",  "late",    1L, 0.090525649,
  "herbicide",  "late",    1L, 0.070546417,
  "herbicide",  "late",    1L, 0.150463344,
  "herbicide",  "late",    1L, 0.070546417,
  "herbicide",  "late",    1L, 0.220390654,
  "herbicide",  "late",    1L, 0.120494496,
  "herbicide",  "late",    1L, 0.150463344,
  "herbicide",  "late",    1L, 0.130484112,
  "herbicide",  "late",    1L, 0.090525649,
  "herbicide",  "late",    1L, 0.020598339,
  "herbicide",  "late",    1L, 0.170442575,
  "herbicide",  "late",    1L, 0.050567186,
  "herbicide", "early",    1L, 0.010608723,
  "herbicide", "early",    1L, 0.060556802,
  "herbicide", "early",    1L, 0.000619107,
  "herbicide", "early",    1L, 0.050567186,
  "herbicide", "early",    1L, 0.030587954,
  "herbicide", "early",    1L, 0.010608723,
  "herbicide", "early",    1L, 0.000619107,
  "herbicide", "early",    1L, 0.000619107,
  "herbicide", "early",    1L, 0.020598339,
  "herbicide", "early",    1L, 0.000619107,
  "herbicide", "early",    1L, 0.030587954,
  "herbicide", "early",    1L, 0.010608723,
  "herbicide",   "mid",    1L,  0.04057757,
  "herbicide",   "mid",    1L, 0.050567186,
  "herbicide",   "mid",    1L, 0.010608723,
  "herbicide",   "mid",    1L, 0.010608723,
  "herbicide",   "mid",    1L,  0.04057757,
  "herbicide",   "mid",    1L, 0.010608723,
  "herbicide",   "mid",    1L, 0.050567186,
  "herbicide",   "mid",    1L, 0.010608723,
  "herbicide",   "mid",    1L, 0.010608723,
  "herbicide",   "mid",    1L, 0.070546417,
  "herbicide",   "mid",    1L, 0.020598339,
  "herbicide",   "mid",    1L, 0.060556802,
  "herbicide",  "late",    1L, 0.030587954,
  "herbicide",  "late",    1L, 0.030587954,
  "herbicide",  "late",    1L, 0.070546417,
  "herbicide",  "late",    1L,  0.04057757,
  "herbicide",  "late",    1L, 0.010608723,
  "herbicide",  "late",    1L, 0.080536033,
  "herbicide",  "late",    1L, 0.000619107,
  "herbicide",  "late",    1L, 0.010608723,
  "herbicide",  "late",    1L, 0.010608723,
  "herbicide",  "late",    1L, 0.030587954,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L,  0.04057757,
  "mow", "early",    1L,  0.04057757,
  "mow", "early",    1L, 0.080536033,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L, 0.020598339,
  "mow", "early",    1L, 0.060556802,
  "mow", "early",    1L, 0.000619107,
  "mow", "early",    1L,  0.04057757,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L, 0.020598339,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L, 0.070546417,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L,  0.04057757,
  "mow",   "mid",    1L,  0.04057757,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L, 0.030587954,
  "mow",   "mid",    1L, 0.010608723,
  "mow",   "mid",    1L, 0.010608723,
  "mow",  "late",    1L,  0.04057757,
  "mow",  "late",    1L, 0.020598339,
  "mow",  "late",    1L,  0.04057757,
  "mow",  "late",    1L, 0.020598339,
  "mow",  "late",    1L, 0.020598339,
  "mow",  "late",    1L, 0.020598339,
  "mow",  "late",    1L, 0.030587954,
  "mow",  "late",    1L, 0.030587954,
  "mow",  "late",    1L, 0.020598339,
  "mow",  "late",    1L, 0.000619107,
  "mow",  "late",    1L, 0.030587954,
  "mow",  "late",    1L, 0.030587954,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L, 0.010608723,
  "mow", "early",    1L, 0.100515265,
  "mow", "early",    1L, 0.110504881,
  "mow", "early",    1L,  0.04057757,
  "mow", "early",    1L, 0.030587954,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L,  0.04057757,
  "mow", "early",    1L, 0.050567186,
  "mow", "early",    1L, 0.010608723,
  "mow", "early",    1L, 0.010608723,
  "mow", "early",    1L, 0.000619107,
  "mow",   "mid",    1L, 0.060556802,
  "mow",   "mid",    1L, 0.010608723,
  "mow",   "mid",    1L, 0.000619107,
  "mow",   "mid",    1L, 0.030587954,
  "mow",   "mid",    1L, 0.060556802,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L, 0.050567186,
  "mow",   "mid",    1L,  0.04057757,
  "mow",   "mid",    1L, 0.020598339,
  "mow",   "mid",    1L,  0.04057757,
  "mow",   "mid",    1L, 0.030587954,
  "mow",   "mid",    1L, 0.030587954,
  "mow",  "late",    1L, 0.050567186,
  "mow",  "late",    1L, 0.050567186,
  "mow",  "late",    1L, 0.010608723,
  "mow",  "late",    1L, 0.030587954,
  "mow",  "late",    1L, 0.010608723,
  "mow",  "late",    1L, 0.010608723,
  "mow",  "late",    1L, 0.060556802,
  "mow",  "late",    1L, 0.020598339,
  "mow",  "late",    1L, 0.050567186,
  "mow",  "late",    1L,  0.04057757,
  "mow",  "late",    1L, 0.010608723,
  "mow",  "late",    1L, 0.070546417,
  "herbicide", "early",    2L,  0.04057757,
  "herbicide", "early",    2L, 0.450151817,
  "herbicide", "early",    2L, 0.000619107,
  "herbicide", "early",    2L, 0.500099896,
  "herbicide", "early",    2L, 0.010608723,
  "herbicide", "early",    2L, 0.190421807,
  "herbicide", "early",    2L, 0.180432191,
  "herbicide", "early",    2L, 0.130484112,
  "herbicide", "early",    2L, 0.020598339,
  "herbicide", "early",    2L, 0.360245275,
  "herbicide", "early",    2L, 0.010608723,
  "herbicide", "early",    2L, 0.030587954,
  "herbicide",   "mid",    2L, 0.050567186,
  "herbicide",   "mid",    2L, 0.370234891,
  "herbicide",   "mid",    2L, 0.010608723,
  "herbicide",   "mid",    2L, 0.250359502,
  "herbicide",   "mid",    2L, 0.050567186,
  "herbicide",   "mid",    2L, 0.080536033,
  "herbicide",   "mid",    2L,  0.04057757,
  "herbicide",   "mid",    2L, 0.050567186,
  "herbicide",   "mid",    2L, 0.050567186,
  "herbicide",   "mid",    2L,  0.16045296,
  "herbicide",   "mid",    2L, 0.000619107,
  "herbicide",   "mid",    2L, 0.000619107,
  "herbicide",  "late",    2L, 0.050567186,
  "herbicide",  "late",    2L, 0.540058359,
  "herbicide",  "late",    2L,  0.04057757,
  "herbicide",  "late",    2L, 0.260349117,
  "herbicide",  "late",    2L, 0.070546417,
  "herbicide",  "late",    2L, 0.120494496,
  "herbicide",  "late",    2L, 0.030587954,
  "herbicide",  "late",    2L, 0.070546417,
  "herbicide",  "late",    2L, 0.020598339,
  "herbicide",  "late",    2L, 0.120494496,
  "herbicide",  "late",    2L,  0.04057757,
  "herbicide",  "late",    2L, 0.000619107,
  "herbicide", "early",    2L, 0.010608723,
  "herbicide", "early",    2L, 0.050567186,
  "herbicide", "early",    2L, 0.010608723,
  "herbicide", "early",    2L, 0.010608723,
  "herbicide", "early",    2L, 0.060556802,
  "herbicide", "early",    2L,  0.04057757,
  "herbicide", "early",    2L, 0.210401038,
  "herbicide", "early",    2L, 0.060556802,
  "herbicide", "early",    2L, 0.100515265,
  "herbicide", "early",    2L, 0.090525649,
  "herbicide", "early",    2L, 0.010608723,
  "herbicide", "early",    2L, 0.000619107,
  "herbicide",   "mid",    2L, 0.060556802,
  "herbicide",   "mid",    2L, 0.020598339,
  "herbicide",   "mid",    2L, 0.030587954,
  "herbicide",   "mid",    2L, 0.010608723,
  "herbicide",   "mid",    2L, 0.000619107,
  "herbicide",   "mid",    2L, 0.010608723,
  "herbicide",   "mid",    2L, 0.030587954,
  "herbicide",   "mid",    2L, 0.070546417,
  "herbicide",   "mid",    2L, 0.020598339,
   
     

    library(tidyverse)
    library(betareg)
    library(emmeans)
    library(lmtest)
    library(multcomp)
    library(lme4)
    library(car)
    library(glmmTMB)
    
    trt_key <- c(ctrl = "Control", mow = "FallMow", herbicide = "SpotSpray")
    aq$trt <- recode(aq$trt, !!!trt_key)  
    aq$trt <- factor(aq$trt, levels = c("Control", "FallMow", "SpotSpray"))
    season_key <- c(early = "Early", mid = "Mid", late = "Late")
    aq$season <- recode(aq$season, !!!season_key)
    aq$season <- factor(aq$season, levels=c("Early","Mid","Late"))
    
    glm.soil <- glmmTMB(sp ~ trt + season + trt*season + (1 | site), data = aq, 
                     family = list(family = "beta", link = "logit"), dispformula = ~trt)
    
    #Interaction
    lsm <- emmeans(glm.soil, pairwise ~ trt:season, type="response", adjust = "tukey")
    lsmtab <- cld(lsm, Letter=letters, sort = F)
    colnames(lsmtab)[1] <- "Treatment"
    colnames(lsmtab)[2] <- "Season"
    colnames(lsmtab)[8] <- "letter"
    
    df <- as.data.frame(lsmtab)
    print(df)

This is my first post, so I apologize in advance if I've overlooked any posting protocols. Thanks!


  [1]: https://i.stack.imgur.com/GJ8VA.png
  [2]: https://schmidtpaul.github.io/DSFAIR/compactletterdisplay.html

Answer 1

I'm the author of the plot/code you linked. You are not the first one asking how to create an analogous plot when interactions are present. I am suggesting two options below using your data.

^{(Note that in the following reprex I deleted the part of the code with the data, because my post was reaching the character limit.)}

# packages  ---------------------------------------------------------------
library(emmeans)
library(glmmTMB)
library(car)
library(multcomp)
library(multcompView)
library(scales)
library(tidyverse)

# format ------------------------------------------------------------------
trt_key <- c(Control = "ctrl", FallMow = "mow", SpotSpray = "herbicide")
season_key <- c(Early = "early", Mid = "mid", Late = "late")

aq <- aq %>%
  mutate(
    trt    = trt    %>% fct_recode(!!!trt_key)    %>% fct_relevel(names(trt_key)),
    season = season %>% fct_recode(!!!season_key) %>% fct_relevel(names(season_key))
  )

# model -------------------------------------------------------------------
glm.soil <-
  glmmTMB(
    sp ~ trt + season + trt:season + (1 | site),
    data = aq,
    family = list(family = "beta", link = "logit"),
    dispformula = ~ trt
  )

Anova(glm.soil) # interaction is significant!
#> Analysis of Deviance Table (Type II Wald chisquare tests)
#> 
#> Response: sp
#>             Chisq Df Pr(>Chisq)    
#> trt        48.422  2  3.057e-11 ***
#> season     18.888  2  7.916e-05 ***
#> trt:season 16.980  4   0.001951 ** 
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here the ANOVA tells us that there are significant interaction effects between treatment and season. Simply put, this means that treatments behave differently depending on the season. In the extreme case this could mean that the same treatment could have the best performance in one season, but the worst performance in another season. Therefore it would be misleading to simply estimate one mean per treatment across seasons (or one mean per season across treatments). Instead, we should look at the means of all treatment-season combinations. In this scenario, there are 9 season-treatment combinations and we can estimate them using emmeans() via either ~ trt:season or ~ trt|season. These are the two options I was talking about above. Again, both estimate the same means for the same 9 combinations, but what is different is which of these means are compared to each other. None of the two appraoches is more correct than the other. Instead, you as the analyst must decide which approach is more informative for what you are trying to investigate.

# emmean comparisons ------------------------------------------------------
emm1 <- emmeans(glm.soil, ~ trt:season, type = "response")
emm2 <- emmeans(glm.soil, ~ trt|season, type = "response")

cld1 <- cld(emm1, Letter=letters, adjust = "tukey")
cld2 <- cld(emm2, Letter=letters, adjust = "tukey")

~ trt:season Comparing all combinations to all other combinations

Here, all combinations are compared to all other combinations. In order to plot this, I create a new column trt_season that represents each combination (and I also sort its factor levels according to their estimated mean) and put it on the x-axis. Note that I also filled the boxes and colored the dots according to their treatment, but this is optional.

cld1
#>  trt       season response      SE  df lower.CL upper.CL .group
#>  Control   Early    0.0560 0.00540 950   0.0428   0.0730  abc  
#>  FallMow   Early    0.0316 0.00248 950   0.0254   0.0393     d 
#>  SpotSpray Early    0.0498 0.00402 950   0.0397   0.0622  abc  
#>  Control   Mid      0.0654 0.00610 950   0.0504   0.0845  a    
#>  FallMow   Mid      0.0427 0.00305 950   0.0350   0.0520   b   
#>  SpotSpray Mid      0.0609 0.00473 950   0.0490   0.0754  a c  
#>  Control   Late     0.0442 0.00464 950   0.0330   0.0590   bcd 
#>  FallMow   Late     0.0437 0.00314 950   0.0358   0.0533   b   
#>  SpotSpray Late     0.0630 0.00482 950   0.0508   0.0777  a c  
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 9 estimates 
#> Intervals are back-transformed from the logit scale 
#> P value adjustment: tukey method for comparing a family of 9 estimates 
#> Tests are performed on the log odds ratio scale 
#> significance level used: alpha = 0.05 
#> NOTE: Compact letter displays can be misleading
#>       because they show NON-findings rather than findings.
#>       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

cld1_df <- cld1 %>% as.data.frame()

cld1_df <- cld1_df %>% 
  mutate(trt_season = paste0(season, "-", trt)) %>% 
  mutate(trt_season = fct_reorder(trt_season, response))

aq <- aq %>% 
  mutate(trt_season = paste0(season, "-", trt)) %>% 
  mutate(trt_season = fct_relevel(trt_season, levels(cld1_df$trt_season)))

ggplot() +
  # y-axis
  scale_y_continuous(
    name = "sp",
    limits = c(0, NA),
    breaks = pretty_breaks(),
    expand = expansion(mult = c(0, 0.1))
  ) +
  # x-axis
  scale_x_discrete(name = "Treatment-Season combination") +
  # general layout
  theme_bw() +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    vjust = 1
  ),
  legend.position = "top") +
  # black data points
  geom_point(
    data = aq,
    aes(y = sp, x = trt_season, color = trt),
    #shape = 16,
    alpha = 0.5,
    position = position_nudge(x = -0.2)
  ) +
  # black boxplot
  geom_boxplot(
    data = aq,
    aes(y = sp, x = trt_season, fill = trt),
    width = 0.05,
    outlier.shape = NA,
    position = position_nudge(x = -0.1)
  ) +
  # red mean value
  geom_point(
    data = cld1_df,
    aes(y = response, x = trt_season),
    size = 2,
    color = "red"
  ) +
  # red mean errorbar
  geom_errorbar(
    data = cld1_df,
    aes(ymin = lower.CL, ymax = upper.CL, x = trt_season),
    width = 0.05,
    color = "red"
  ) +
  # red letters
  geom_text(
    data = cld1_df,
    aes(
      y = response,
      x = trt_season,
      label = str_trim(.group)
    ),
    position = position_nudge(x = 0.1),
    hjust = 0,
    color = "red"
  ) +
  labs(
    caption = str_wrap("Black dots represent raw data. Red dots and error bars represent backtransformed estimated marginal means ± 95% confidence interval per group. Means not sharing any letter are significantly different by the Tukey-test at the 5% level of significance.", width = 100)
  )

~ trt|season Comparing all trt combinations only within each season

Here, fewer mean comparisons are made, i.e. only 3 comparisons (Control vs. FallMow, FallMow vs. SpotSpray, Control vs. SpotSpray) for each of the seasons. This means that e.g. Early-Control is never compared to Mid-Control. Moreover, the letters from the compact letter display are also created separately for each season. This means that e.g. the a assigned to the Early-Control mean has nothing to do with the a assigned to the Mid-Control mean. This is crucial and I made sure to explicitly state this in the plot's caption. I also used facettes which separate the results for the seasons visually. That being said, presenting the results in this way may actually be more suited for your investigation. (Obviously you could use colors per treatment or season here as well)

cld2
#> season = Early:
#>  trt       response      SE  df lower.CL upper.CL .group
#>  Control     0.0560 0.00540 950   0.0444   0.0704  a    
#>  FallMow     0.0316 0.00248 950   0.0262   0.0381   b   
#>  SpotSpray   0.0498 0.00402 950   0.0410   0.0603  a    
#> 
#> season = Mid:
#>  trt       response      SE  df lower.CL upper.CL .group
#>  Control     0.0654 0.00610 950   0.0522   0.0816  a    
#>  FallMow     0.0427 0.00305 950   0.0360   0.0506   b   
#>  SpotSpray   0.0609 0.00473 950   0.0505   0.0733  a    
#> 
#> season = Late:
#>  trt       response      SE  df lower.CL upper.CL .group
#>  Control     0.0442 0.00464 950   0.0344   0.0567  a    
#>  FallMow     0.0437 0.00314 950   0.0368   0.0519  a    
#>  SpotSpray   0.0630 0.00482 950   0.0524   0.0755   b   
#> 
#> Confidence level used: 0.95 
#> Conf-level adjustment: sidak method for 3 estimates 
#> Intervals are back-transformed from the logit scale 
#> P value adjustment: tukey method for comparing a family of 3 estimates 
#> Tests are performed on the log odds ratio scale 
#> significance level used: alpha = 0.05 
#> NOTE: Compact letter displays can be misleading
#>       because they show NON-findings rather than findings.
#>       Consider using 'pairs()', 'pwpp()', or 'pwpm()' instead.

cld2_df <- cld2 %>% as.data.frame()

ggplot() +
  facet_grid(cols = vars(season), labeller = label_both) +
  # y-axis
  scale_y_continuous(
    name = "sp",
    limits = c(0, NA),
    breaks = pretty_breaks(),
    expand = expansion(mult = c(0, 0.1))
  ) +
  # x-axis
  scale_x_discrete(name = "Treatment") +
  # general layout
  theme_bw() +
  # black data points
  geom_point(
    data = aq,
    aes(y = sp, x = trt),
    shape = 16,
    alpha = 0.5,
    position = position_nudge(x = -0.2)
  ) +
  # black boxplot
  geom_boxplot(
    data = aq,
    aes(y = sp, x = trt),
    width = 0.05,
    outlier.shape = NA,
    position = position_nudge(x = -0.1)
  ) +
  # red mean value
  geom_point(
    data = cld2_df,
    aes(y = response, x = trt),
    size = 2,
    color = "red"
  ) +
  # red mean errorbar
  geom_errorbar(
    data = cld2_df,
    aes(ymin = lower.CL, ymax = upper.CL, x = trt),
    width = 0.05,
    color = "red"
  ) +
  # red letters
  geom_text(
    data = cld2_df,
    aes(
      y = response,
      x = trt,
      label = str_trim(.group)
    ),
    position = position_nudge(x = 0.1),
    hjust = 0,
    color = "red"
  ) +
  labs(
    caption = str_wrap("Black dots represent raw data. Red dots and error bars represent backtransformed estimated marginal means ± 95% confidence interval per group. For each season separately, means not sharing any letter are significantly different by the Tukey-test at the 5% level of significance.", width = 100)
  )

So these are the two options I think about when I have a significant two-way interaction and want to compare means. Note that you could also switch treatment and season to ~ season|trt and also plot it the other way around.

Further reading

• emmeans documentation
• chapter on Compact Letter Display (CLD) - you will also find a brief discussion on the weakness of CLDs that Russ Lenth mentioned in his comment
• related stackoverflow post #1
• related stackoverflow post #2

Bonus: Raincloud plot

This has nothing to do with your question, but I'd like to point out that you have much more data points than I have in the the plot you linked. Because of this, the original plotting approach could potentially be improved because geom_point() leads to many dots being plotted on top of each other so that we have no way seeing how much data there really is. Therefore, you could simply replace geom_point() with geom_jitter() or even go further and create these raincloud plots mentioned in this blogpost. I've created one (without the emmeans part) below that is analogous to the first plot above.

ggplot() +
  # y-axis
  scale_y_continuous(
    name = "sp",
    limits = c(0, NA),
    breaks = pretty_breaks(),
    expand = expansion(mult = c(0, 0.1))
  ) +
  # x-axis
  scale_x_discrete(name = "Treatment-Season combination") +
  # general layout
  theme_bw() +
  theme(axis.text.x = element_text(
    angle = 45,
    hjust = 1,
    vjust = 1
  ),
  legend.position = "top") +
  # add half-violin from {ggdist} package
  ggdist::stat_halfeye(
    data = aq,
    aes(y = sp, x = trt_season, fill = trt),
    adjust = .5, 
    width = .5, 
    .width = 0, 
    justification = -.2, 
    point_colour = NA, 
    show.legend = FALSE
  ) + 
  # boxplot
  geom_boxplot(
    data = aq,
    aes(y = sp, x = trt_season, fill = trt),
    width = 0.1,
    outlier.shape = NA
  ) +
  # add justified jitter from the {gghalves} package
  gghalves::geom_half_point(
    data = aq,
    aes(y = sp, x = trt_season, color = trt),
    side = "l", 
    range_scale = .4, 
    alpha = .2
  )

^{Created on 2022-01-26 by the reprex package (v2.0.1)}

Answer 2

IMO, almost anything is better than a CLD. They display non-findings rather than findings.

I suggest presenting the simple comparisons in tabular form

lsm = emmeans(glm.soil, ~ trt:season, type = "response")
pairs(lsm, by = "trt")
pairs(lsm, by = "season")

If you really want the diagonal comparisons, you could consider pwpp() or pwpm() which can display quite a few comparisons compactly. See the documentation in emmeans