Multiple comparisons with gtsummary

Question

Since my question is similar to one that's been asked before, I'll steal the reprex (also below), for consistency's sake, from Summary Table (mean + std.error) with p-values for 2-way anova

I'm curious how to integrate a post hoc means comparisons (i.e. multcomp) and display letter groupings, like what the compact letter display function cld() would provide, directly in the gtsummary table.

Check out this table as an example of what I'm trying to achieve. But ideally, I'd like to use superscripts to denote letter groupings:

Wine grape example

library(gtsummary)
library(titanic)
library(tidyverse)
library(plotrix) #has a std.error function

packageVersion("gtsummary")
#> [1] '1.4.0'

# create smaller version of the dataset
df <- 
  titanic_train %>%
  select(Sex, Embarked, Age, Fare) %>%
  filter(Embarked != "") # deleting empty Embarked status

# first, write a little function to get the 2-way ANOVA p-values in a table
# function to get 2-way ANOVA p-values in tibble
twoway_p <- function(variable) {
  paste(variable, "~ Sex * Embarked") %>%
    as.formula() %>%
    aov(data = df) %>% 
    broom::tidy() %>%
    select(term, p.value) %>%
    filter(complete.cases(.)) %>%
    pivot_wider(names_from = term, values_from = p.value) %>%
    mutate(
      variable = .env$variable,
      row_type = "label"
    )
}

# add all results to a single table (will be merged with gtsummary table in next step)
twoway_results <-
  bind_rows(
    twoway_p("Age"),
    twoway_p("Fare")
  )
twoway_results
#> # A tibble: 2 x 5
#>            Sex Embarked `Sex:Embarked` variable row_type
#>          <dbl>    <dbl>          <dbl> <chr>    <chr>   
#> 1 0.00823      3.97e- 1         0.611  Age      label   
#> 2 0.0000000191 4.27e-16         0.0958 Fare     label


tbl <-
  # first build a stratified `tbl_summary()` table to get summary stats by two variables
  df %>%
  tbl_strata(
    strata =  Sex,
    .tbl_fun =
      ~.x %>%
      tbl_summary(
        by = Embarked,
        missing = "no",
        statistic = all_continuous() ~ "{mean} ({std.error})",
        digits = everything() ~ 1
      ) %>%
      modify_header(all_stat_cols() ~ "**{level}**")
  ) %>%
  # merge the 2way ANOVA results into tbl_summary table
  modify_table_body(
    ~.x %>%
      left_join(
        twoway_results,
        by = c("variable", "row_type")
      )
  ) %>%
  # by default the new columns are hidden, add a header to unhide them
  modify_header(list(
    Sex ~ "**Sex**", 
    Embarked ~ "**Embarked**", 
    `Sex:Embarked` ~ "**Sex * Embarked**"
  )) %>%
  # adding spanning header to analysis results
  modify_spanning_header(c(Sex, Embarked, `Sex:Embarked`) ~ "**Two-way ANOVA p-values**") %>%
  # format the p-values with a pvalue formatting function
  modify_fmt_fun(c(Sex, Embarked, `Sex:Embarked`) ~ style_pvalue) %>%
  # update the footnote to be nicer looking
  modify_footnote(all_stat_cols() ~ "Mean (SE)")