在 R 中是否有一种方法可以通过方差分析和 Tukeys 列表循环 multcompview 函数（即 CLD）？

Question

我最终想用 CLD 绘制每个变量。但与此同时，我在获取 CLD 时遇到了麻烦。我有这样的数据；

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))

感谢上一个问题 here 的非常有用的答案，然后我通过每个变量（并制作数据框）循环方差分析和 Tukey 像这样；

## 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)))

现在我只需要获取每个变量的 CLD。我通常使用

multcompView::multcompLetters4()

，但我知道还有其他选择。我也想不出

lapply

或

mapply

来处理两个单独的列表。最终我很想

geom_boxplot()

结果，我会完全迷失，但步步为营。非常感谢任何帮助！

Answer 1

你能使用

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)

^{创建于 2023-02-23 与 reprex v2.0.2}

Answer 2

这是使用

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)

Answer 3

要获得 CLD，您可以首先将 'aov_res' 传递给 emmeans 包中的 emmeans() 函数，以获得 SE 的边际均值和 置信限度。然后此输出将用作 mulicomp 和 multicompview 的 cld() 函数的所需对象 添加字母以比较带有紧凑字母显示的零食的包装。

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 
CLD_letters <- lapply(CLD
                      ,cld  #Set up a compact letter display of all pair-wise comparisons
                      ,adjust= 'fdr' # adjusting p-values using false dicovery rate (other methods are tukey, 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