Two-Way Comparison Figures with ggcompare()

A common publication figure compares an outcome across two categorical factors at once – for example a score across gender and education level – and annotates it with the interaction from a two-way ANOVA plus pairwise comparisons within each group. ggcompare() assembles this figure in a single call.

The one-call two-way figure

Map the first factor to x and the second factor to color (or fill). That second mapping switches ggcompare() into its two-way mode.

library(ggpubr)
library(rstatix)
#> 
#> Attaching package: 'rstatix'
#> The following object is masked from 'package:ggcorrplot':
#> 
#>     cor_pmat
#> The following object is masked from 'package:stats':
#> 
#>     filter

data("jobsatisfaction", package = "datarium")

ggcompare(
  jobsatisfaction,
  x = "gender", y = "score", color = "education_level"
)

The figure carries:

  • dodged boxes for the two factors, with jittered points and group means;
  • pairwise comparison brackets of the education levels within each gender;
  • a subtitle naming the two-way ANOVA interaction; and
  • a caption describing the pairwise test.

In two-way mode the pairwise method defaults to "emmeans_test" (simple comparisons from the pooled model), p.adjust.method to "bonferroni", and hide.ns to TRUE. Each of these is a default only – pass the argument to override it.

It matches the by-hand analysis

The annotation reproduces the numbers you would compute by hand. The subtitle is the interaction row of the two-way ANOVA:

jobsatisfaction %>% anova_test(score ~ gender * education_level)
#> ANOVA Table (type II tests)
#> 
#>                   Effect DFn DFd       F        p p<.05   ges
#> 1                 gender   1  52   0.745 3.92e-01       0.014
#> 2        education_level   2  52 187.892 1.60e-24     * 0.878
#> 3 gender:education_level   2  52   7.338 2.00e-03     * 0.220

and the brackets are the pooled-model emmeans comparisons within each gender:

jobsatisfaction %>%
  group_by(gender) %>%
  emmeans_test(score ~ education_level, p.adjust.method = "bonferroni") %>%
  select(gender, group1, group2, p.adj, p.adj.signif)
#> # A tibble: 6 × 5
#>   gender group1  group2        p.adj p.adj.signif
#>   <fct>  <chr>   <chr>         <dbl> <chr>       
#> 1 male   school  college    1.01e- 2 *           
#> 2 male   school  university 2.06e-20 ****        
#> 3 male   college university 2.53e-16 ****        
#> 4 female school  college    1.49e- 2 *           
#> 5 female school  university 1.82e-14 ****        
#> 6 female college university 5.52e-10 ****

Simple main effects on the plot

When the interaction is significant, interpret it by breaking it into simple main effects – the effect of the focal factor within each level of the other. Group by the moderator and run a one-way anova_test(), passing the full two-way model as error so the test borrows the pooled error term from the whole design:

model <- lm(score ~ gender * education_level, data = jobsatisfaction)

sme <- jobsatisfaction %>%
  group_by(gender) %>%
  anova_test(score ~ education_level, error = model)
sme
#> # A tibble: 2 × 8
#>   gender Effect            DFn   DFd     F        p `p<.05`   ges
#> * <fct>  <chr>           <dbl> <dbl> <dbl>    <dbl> <chr>   <dbl>
#> 1 male   education_level     2    52 132.  3.92e-21 *       0.836
#> 2 female education_level     2    52  62.8 1.35e-14 *       0.707

The pooled error is what keeps the denominator degrees of freedom at the full-model residual df – here F(2, 52) for both genders. Dropping error = model would instead use each subset’s own residuals, giving a different (smaller) df and a different F, so keep the error argument when you report a simple main effect.

Layer those per-group F values onto the ggcompare() figure with a geom_text(), placing each label above its cluster’s brackets:

p <- ggcompare(
  jobsatisfaction,
  x = "gender", y = "score", color = "education_level"
)

# position each label just above that group's tallest bracket
pwc <- jobsatisfaction %>%
  group_by(gender) %>%
  emmeans_test(score ~ education_level, p.adjust.method = "bonferroni") %>%
  add_xy_position(x = "gender")
tops <- tapply(pwc$y.position, pwc$gender, max)

sme_lab <- transform(
  sme,
  x = as.integer(factor(gender, levels = levels(jobsatisfaction$gender))),
  y = tops[as.character(gender)] + 0.4,
  label = sprintf(
    "F(%g,%g) = %.3g, %s", DFn, DFd, F,
    ifelse(p < 0.0001, "p < 0.0001", paste0("p = ", signif(p, 2)))
  )
)

p +
  geom_text(
    data = sme_lab, inherit.aes = FALSE,
    aes(x = x, y = y, label = label), size = 3, fontface = "italic"
  ) +
  scale_y_continuous(expand = expansion(mult = c(0.05, 0.28)))

The figure now reports itself completely: the interaction on top, each gender’s simple main effect of education above its boxes, and which education levels differ within each gender marked by the brackets.

Customizing

Every element is overridable. For instance, compare the genders within each education level instead (swap the comparison direction with pwc.group.by), switch the base geometry, and show non-significant brackets:

ggcompare(
  jobsatisfaction,
  x = "gender", y = "score", color = "education_level",
  base = "violin",
  pwc.group.by = "legend.var",
  hide.ns = FALSE
)

You can also compare each group to a reference level (ref.group), choose a different pairwise method/p.adjust.method, or pass extra styling through pwc.args, add.params and ... (forwarded to the base builder). See ?ggcompare for the full argument list.

Naming an effect on any plot

The subtitle machinery is also available directly through add_test_label(). Given a plot with x and a second factor, group.by computes the two-way ANOVA and effect chooses which effect to report – the interaction (default), a named term, or "all" for a multi-line label of both main effects and the interaction:

p <- ggboxplot(jobsatisfaction, x = "gender", y = "score",
  color = "education_level")

add_test_label(p, group.by = "education_level", effect = "all", type = "text")