In R is there a way to loop multcompview functions (i.e. CLD) through lists of ANOVAs and Tukeys?

I eventually want to plot each variable with CLDs. But in the mean time I am having trouble just getting the CLDs. I have data that looks like this;

df <- data.frame(Treat = rep(LETTERS[1:4], 100, replace = TRUE),
                 A = rnorm(400),
                 B = rnorm(400),
                 C = rnorm(400),
                 D = rnorm(400),
                 E = rnorm(400),
                 F = rnorm(400),
                 G = rnorm(400),
                 H = rnorm(400))

Thanks to a very helpful answer from a previous question here, I then looped ANOVA's and Tukey's through each variable (and made data frames) like this;

## perform anova on dependent variables
aov_res <- apply(df[,2:ncol(df)], 2, function(x) aov(x ~ Treat, data = df))

# Apply Tukey's HSD test to the results of each ANOVA test
tukey_res <- sapply(aov_res, function(x) TukeyHSD(x, "Treat", ordered = TRUE))

# Convert the results of each ANOVA test into a tidy data frame using the broom package
aov_res_df <- do.call(rbind, lapply(aov_res, broom::tidy))

# Combine the results of the Tukey HSD tests into a single data frame
tukey_res_df <- as.data.frame(do.call(rbind, Map(cbind, Name = names(tukey_res), tukey_res)))

Now I just need to get the CLDs for each variable. I usually use multcompView::multcompLetters4() but I know there are alternatives. I also couldn't figure out lapply or mapply for working with two separate lists. Eventually I would love to geom_boxplot() the results, and I will be totally lost, but baby steps. Any help is greatly appreciated!

Here is another way using multcompView::multcompLetters4():

## perform anova on dependent variables
aov_res <- apply(df[,2:ncol(df)], 2, function(x) aov(x ~ Treat, data = df))

# Apply Tukey's HSD test to the results of each ANOVA test using lapply instead of sapply
tukey_res <- lapply(aov_res, function(x) TukeyHSD(x, "Treat", ordered = TRUE))

clds <- lapply(names(aov_res), function(x) multcompLetters4(aov_res[[x]], tukey_res[[x]]))

# tidy up 
names(clds) <- names(aov_res)
clds <- as.data.frame(clds)

To get the CLDs you can pass the 'aov_res' to first, the emmeans() function from emmeans package to obtain the marginal means with SEs and
confidence limits. Then this output would be used as a desired object for cld() function from mulicomp and multicompview
packages which add letters to compare the Treats with compact letter display.If you wish to do adjustment with Tukey method, simply remove adjust= option in the cld() function.

library(multcomp)
library(multcompView)
library(emmeans)
library(tidyverse)

df$Treat<-factor(df$Treat)

## perform anova on dependent variables
aov_res <- apply(df[,2:ncol(df)], 2, function(x) aov(x ~ Treat, data = df))

#Estimate marginal means (Least-squares means)
CLD <- lapply(aov_res ,FUN = function(i) emmeans::emmeans(i, specs = "Treat"))

# add your compact letters to the estimated means, note that "tukey" is only appropriate for one set of pairwise comparisons 
CLD_letters <- lapply(CLD
                      ,cld  #Set up a compact letter display of all pair-wise comparisons
                      ,adjust= 'sidak' # adjusting p-values using sidak method (other methods are fdr, BH,...)
                      ,Letters = letters # grouping with small letters
                      ,alpha = 0.05 # your significance level
                      ,reversed = T) # sort means if larger means are desired


# convert list of comparisons to dataframe including all informations on treatment means and for each dependent variable
CLD_letters_df <- bind_rows(CLD_letters, .id = "variable")

CLD_letters_df

Could you use the functions in multcomp:

library(multcomp)
library(broom)
library(dplyr)
df <- data.frame(Treat = rep(LETTERS[1:4], 100, replace = TRUE),
                 A = rnorm(400),
                 B = rnorm(400),
                 C = rnorm(400),
                 D = rnorm(400),
                 E = rnorm(400),
                 F = rnorm(400),
                 G = rnorm(400),
                 H = rnorm(400))

df$Treat <- as.ordered(df$Treat)

## perform anova on dependent variables
aov_res <- apply(df[,2:ncol(df)], 2, function(x) aov(x ~ Treat, data = df))

# Apply Tukey's HSD test to the results of each ANOVA test
tukey_res <- lapply(aov_res, function(a)summary(glht(a, linfct = mcp(Treat = "Tukey"))))

# Convert the results of each ANOVA test into a tidy data frame using the broom package
aov_res_df <- do.call(rbind, lapply(aov_res, broom::tidy))

# Combine the results of the Tukey HSD tests into a single data frame

tukey_res_df <- bind_rows(lapply(tukey_res, function(t)do.call(data.frame, t$test[-c(1,2)]) %>% 
                   as_tibble(rownames="comparison")), .id = "DV")

clds <- lapply(tukey_res, cld)

^{Created on 2023-02-23 with reprex v2.0.2}