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)