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

Question

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!

score 1 · Answer 1 · answered Feb 24 '23 at 01:18

1

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)

answered Feb 24 '23 at 01:18

asaei

491
3
5

Thanks for a great answer. On my real data I am getting an error that clds cannot be coerced into a dataframe because it is type mulcompview. Even though class(clds) = list. Any idea why the error? – Dustin Feb 24 '23 at 03:58
I am not able to reproduce it. Would you be able to share part of your data using ```dput()```. Please see [here](https://meta.stackoverflow.com/questions/315885/what-is-the-correct-way-to-share-r-data-on-stack-overflow) for more details on how to share data. – asaei Feb 24 '23 at 07:25

S-SHAAF · Accepted Answer · 2023-03-05T00:18:06.060

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

Wow, this is a very elegant solution that gives me exactly what I am looking for. However, I am curious about adjust = 'tukey' giving me an error that this is "only appropriate for one set of pairwise comparisons." To me it looks like this IS only being applied to one set of pairwise comparisons (e.g. Treat A separate from Treat B, etc.). Am I wrong in my understanding? But if not, then why the error? — Dustin, Feb 24 '23 at 03:07
The p-values of the pairwise comparisons are adjusted with the Tukey-method, while the Sidak adjustment is applied to the confidence intervals of the means (columns lower.CL and upper.CL) and It is not a problem with your statistical analysis . If you wish to do adjustment with Tukey method, simply remove adjust= option in the cld() function. — S-SHAAF, Feb 24 '23 at 12:58
Ah, ok, as long as the pairwise comparisons are still Tukey we are in business. Thank you! — Dustin, Feb 24 '23 at 22:12

score 0 · Answer 3 · answered Feb 24 '23 at 00:50

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}

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

3 Answers3

Linked