从 Tukey HSD 测试中获取 CDL

Posted 2023-02-18

【中文标题】从 Tukey HSD 测试中获取 CDL

【英文标题】：Obtain CLD from TukeyHSD test

【发布时间】：2022-01-08 08:45:26

【问题描述】：

我一直在努力从 TukeyHSD 测试的输出中制作自己的 CLD。

首先我做了一个双向方差分析：

aov2_arbuscular <- aov(arbuscular_count ~ block + pesticide*fertilizer, data = main_trial)

并进行了 TukeyHSD 测试作为事后测试：

tk_arbuscular <- TukeyHSD(aov2_arbuscular)

因为我无法使用 TukeyHSD 输出生成 CLD，所以我使用了 emmeans() 和 cld() 函数。

tk_arbuscular_model <- emmeans(aov2_arbuscular, 
                    pairwise ~ pesticide*fertilizer, 
                    adjust = "tukey")

tk_arbuscular_model_cld <- cld(tk_arbuscular_model$emmeans, 
                    alpha = .05,
                    Letters = letters)

我认为 TukeyHSD 和带有 adjust = "tukey" 的 emmeans 会产生相同的输出。他们这样做是为了获得最大的结果，不幸的是有一些例外。 我已经写了我的结果部分，不想再次调整所有的 p 值。因此，有人可以帮我生成带有 TukeyHSD 输出的 CLD，以便我可以将它们集成到 ggplot 中吗？

【问题讨论】：

为什么你认为TukeyHSD() 给出了不同的答案？我对 emmeans 所做的事情非常有信心，所以这让我认为你的 TukeyHSD 结果是错误的。

【参考方案1】：

您没有提供您的数据，因此我正在创建自己的代表，它也可以作为双因子块设计进行分析，即使用与您相似的模型。

如您所见，我无法重现您的问题 - 所有 p 值基本相同。我注意到您将adjust = "Tukey" 添加到emmeans() 语句中，但我通常将其添加到cld() 语句中——但这不应该是问题。

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

dataURL <- "https://raw.githubusercontent.com/SchmidtPaul/DSFAIR/master/data/Gomez%26Gomez1984.csv"
dat <- read_csv(dataURL) %>% 
  filter(G %in% c("A", "B") & N %in% c("N1", "N2")) %>% 
  mutate_at(vars(rep:N), as.factor)

aov <- aov(yield ~ G + N + G:N + rep, data = dat)


# get contrasts via 3 options ---------------------------------------------
option1 <- stats::TukeyHSD(aov) %>% 
  pluck("G:N") 

option2 <- emmeans::emmeans(aov, ~ G:N) %>% 
  emmeans::pairs(adjust = "Tukey") 

option3 <- emmeans::emmeans(aov, ~ G:N) %>% 
  multcomp::cld(adjust = "Tukey", details = TRUE)


# uniform format ----------------------------------------------------------
option1 <- option1 %>% 
  as_tibble(rownames = "contrast") %>% 
  transmute(contrast = contrast,
            estimate = diff,
            p.value = `p adj`)

option2 <- option2 %>% 
  as_tibble() %>% 
  dplyr::select(contrast, estimate, p.value)

option3 <- option3 %>% 
  pluck("comparisons") %>% 
  as_tibble() %>% 
  dplyr::select(contrast, estimate, p.value)


# compare -----------------------------------------------------------------
option1
#> # A tibble: 6 x 3
#>   contrast  estimate p.value
#>   <chr>        <dbl>   <dbl>
#> 1 B:N1-A:N1     53.3  0.999 
#> 2 A:N2-A:N1   1419.   0.0860
#> 3 B:N2-A:N1   1729.   0.0401
#> 4 A:N2-B:N1   1366    0.0984
#> 5 B:N2-B:N1   1676    0.0455
#> 6 B:N2-A:N2    310    0.910
option2
#> # A tibble: 6 x 3
#>   contrast    estimate p.value
#>   <chr>          <dbl>   <dbl>
#> 1 A N1 - B N1    -53.3  0.999 
#> 2 A N1 - A N2  -1419.   0.0860
#> 3 A N1 - B N2  -1729.   0.0401
#> 4 B N1 - A N2  -1366    0.0984
#> 5 B N1 - B N2  -1676    0.0455
#> 6 A N2 - B N2   -310.   0.910
option3
#> # A tibble: 6 x 3
#>   contrast    estimate p.value
#>   <chr>          <dbl>   <dbl>
#> 1 B N1 - A N1     53.3  0.999 
#> 2 A N2 - A N1   1419.   0.0860
#> 3 A N2 - B N1   1366    0.0984
#> 4 B N2 - A N1   1729.   0.0401
#> 5 B N2 - B N1   1676    0.0455
#> 6 B N2 - A N2    310.   0.910

tibble(
  o1_p = option1$p.value,
  o2_p = option2$p.value,
  o3_p = option3$p.value
) %>% cor()
#>           o1_p      o2_p      o3_p
#> o1_p 1.0000000 1.0000000 0.9967731
#> o2_p 1.0000000 1.0000000 0.9967731
#> o3_p 0.9967731 0.9967731 1.0000000

^{由reprex package (v2.0.1) 于 2021-12-02 创建}

这有帮助吗？如果没有，您能否尝试使用your data 创建a reproducible example，以便我们掌握所有信息？

另外，here 是我写的关于使用和解释紧凑型字母显示的一章。