r - 来自成对比较 p 值的紧凑字母显示 (CLD)

Question

我一直在努力根据 p 值的成对比较表来制作我自己的 CLD。我知道 multcomp 是可能的，但我想生成我自己的 DIY 函数，它可以适应不同的事后输出。当然，有两个具有挑战性的方面:组生成背后的逻辑和编程实现。

我的“群组生成背后的逻辑”是这样的:

1. 按平均顺序处理。
2. 将均值最高的治疗分配给“a”组
3. 从那里开始，循环所有治疗。将每种治疗i与所有现有组中的治疗进行比较。
4. 如果治疗 i 与现有组的所有治疗没有区别，则将其分配给该组。
5. 如果治疗 i 与所有现有组中的至少一个元素不同，则创建一个新组，其中治疗 i 和所有与 i 没有显着差异的治疗已分配。

我使用循环是因为我发现它更容易看到我在做什么。

如果有人能指出逻辑或实现中的任何问题，那就太好了。如果有人能发现任何错误或提供有关如何使其工作的任何提示，我们将不胜感激。

我分别上传数据、每个处理的平均值，以及所有组之间的成对比较(用 agricolae::HSD.test 生成)和 p 值(成对比较文件包括重复数据，因为有是两列中所有可能的组合，这意味着所有处理都显示在两列中)

#LOADING THE DATA
contrasts <- read.table("https://raw.githubusercontent.com/paracon/cld_data/main/contrasts.csv",
                          header=TRUE)

means <- read.table("https://raw.githubusercontent.com/paracon/cld_data/main/means.csv",
                    header=TRUE)

#################################################################################

# WE ORDER TREATMENTS BY MEAN VALUE
means <- means[rev(order(means$mean)),]

# WE CREATE A DATAFRAME WHERE ALL GROUPS WILL BE ADDED AND 
# WHERE THE TREATMENT WITH HIGHEST MEAN IS ASSIGNED GROUP "a"
existing_groups <- rbind(data.frame(trts=character(),
                                    groups=character()), 
                         data.frame(trts=as.character(means$treat[1]),
                                    groups=as.character(letters[1]))
                         )


# WE LOOP ALONG ALL TREATMENTS (FROM means FILE) AFTER THE ONE WITH HIGHEST MEAN
for (i in 2:length(means$treat)){
  
  # WE SUBSET FOR ALL THE CONTRASTS FOR TREATMENT i
  contrasts_i <- contrasts[as.character(contrasts$col1)==as.character(gsub(" ", "", means[i,]$treat, fixed = TRUE)),]

  # WE CREATE AN EMPTY DATAFRAME WHERE WE WILL ADD ALL TREATMENTS NOT DIFFERENT FROM TREATMENT i
  same_as_checked <-  data.frame(trts=as.character(),
                                 groups=as.character())
    
  # WE LOOP ALONG ALL ALREADY EXISTING GROUPS
  for (g in unique(existing_groups$groups)){
    
    # WE SUBSET FOR ALL THE TREATMENTS IN GROUP g
    existing_groups_g <- existing_groups[existing_groups$groups==g,]
    
    # WE LOOP ALONG ALL THE TREATMENTS IN GROUP g
    for (j in 1:length(existing_groups_g$trts)){

      existing_groups_j <- existing_groups_g[j,]
      existing_groups_j$trts <- as.character(gsub(" ", "", existing_groups_j$trts, fixed = TRUE))
      # WE CHECK PAIRWISE COMPARISON BETWEEN TREATMENT j IN THE GROUP AND contrasts_i$col2
      # AND ALL ELEMENTS OF THAT GROUP
      
      try(if(contrasts_i[contrasts_i$col2==existing_groups_j$trts,]$p_val>=0.05){
        
        same_as_checked <-  rbind(same_as_checked,
                                  data.frame(trts=as.character(existing_groups_j$trts),
                                             groups=NA))
      },silent=TRUE)
    }
  }  
  
  print(means[i,]$treat)
  print(same_as_checked$trts)
  
  # same_as_checked SHOULD INCLUDE ALL THE TREATMENTS WHICH ARE NOT DIFFERENT FROM TREATMENT i
  
  
  group_with_no_differences_exists <- "no"
  
  # WE LOOP AGAIN ALONG ALL ALREADY EXISTING GROUPS
  # NOW TO COMPARE THE TREATMENTS IN same_as_checked WITH ALL EXISTING GROUPS
  for (g in unique(existing_groups$groups)){
  
    # WE SUBSET FOR ALL THE TREATMENTS IN GROUP g
    existing_groups_g <- existing_groups[existing_groups$groups==g,]
    
    # WE CHECK IF GROUP g IS IDENTICAL TO same_as_checked
    try(group_with_no_differences_exists <- ifelse(isTRUE(all.equal(unique(same_as_checked[order(same_as_checked$trts),]$trts),
                        unique(existing_groups_g[order(existing_groups_g$groups),]$trts))),
       "yes","no"), silent=TRUE)
    
    # IF GROUP IS IDENTICAL, WE WILL ADD TREATMENT i TO THIS GROUP               
    if (group_with_no_differences_exists=="yes"){
      new_groups <- data.frame(trts=as.character(unique(contrasts_i$col1)),
                        groups=as.character(letters[which(letters==existing_groups_j$groups)])
                        )
    }

    # IF GROUP IS DIFFERENT, WE CREATE A NEW GROUP, WITH TREATMENT i AND ALL THE TREATMENTS IN same_as_checked

    #same_as_checked IS NOT EMPTY:
    if (group_with_no_differences_exists=="no" & nrow(same_as_checked)!=0){
      new_groups <- rbind(data.frame(trts=same_as_checked$trts,
                                     groups=rep(x=as.character(letters[which(letters==existing_groups_j$groups) + 1]),
                                                times=length(same_as_checked$groups))
                                     ),
                          data.frame(trts=as.character(unique(contrasts_i$col1)),
                                     groups=as.character(letters[which(letters==existing_groups_j$groups) + 1])
                                     )
                           )
    }
    #same_as_checked IS EMPTY, TREATMENT i 'IS ALONE' IN THE NEW GROUP:
    if (group_with_no_differences_exists=="no" & nrow(same_as_checked)==0){
      new_groups <- data.frame(trts=as.character(unique(contrasts_i$col1)),
                               groups=as.character(letters[which(letters==existing_groups_j$groups) + 1])
                               )
                               
    }
  }

  existing_groups <- rbind(existing_groups, new_groups)

}

#################################################################################

library(tidyverse)

unique(existing_groups) %>%
  group_by(trts) %>%
  dplyr::summarise(
    groups = paste(as.character(groups), collapse="")
  )

这可行，但会生成冗余组，必须将其删除。要删除这些组，我应用了不太漂亮的代码:

#COMPARING EACH GROUP WITH THE NEXT GROUP (BY ALPHABETICAL ORDER),
#IF ALL TREATMENTS IN A GROUP ARE INCLUDED IN THE NEXT GROUP, THE FORMER IS REMOVED (CONVERTED TO NA)
for (g in unique(as.character(existing_groups$groups)))try({
  existing_trts_g   <- as.character(unique(existing_groups[existing_groups$groups==g,]$trts))
  existing_trts_g_1 <- as.character(unique(existing_groups[existing_groups$groups==as.character(letters[which(letters==g) + 1]),]$trts))
  if (existing_trts_g_1[!(is.na(existing_trts_g_1))] %contain% existing_trts_g[!(is.na(existing_trts_g))]){
    existing_groups[!(is.na(existing_groups$groups)) & existing_groups$groups==g,]$groups <- NA
  }
  existing_groups <- existing_groups[!(is.na(existing_groups$groups)),]
},silent=TRUE)


for (g in unique(as.character(existing_groups$groups))){
  nth <- match(g,unique(as.character(existing_groups$groups)))
  existing_groups[existing_groups$groups==g,]$groups <- letters[nth]
}

预期的组(根据 agricolae::HSD.test)是:

一个:B100；ab: KT50, T10;abc: T35, KT10;公元前:T50、T100、KT35；c: KT100

您也可以在这里看到它们:https://raw.githubusercontent.com/paracon/cld_data/main/output.csv

Answer 1

也许这有帮助。

查找更多详情 here .

mod <- lm(Sepal.Width ~ Species, data = iris)

mod_means_contr <- emmeans::emmeans(object = mod,
                                    pairwise ~ "Species",
                                    adjust = "tukey")

mod_means <- multcomp::cld(object = mod_means_contr$emmeans,
                           Letters = letters)

### Bonus plot
library(ggplot2)

ggplot(data = mod_means,
       aes(x = Species, y = emmean)) +
  geom_point() +
  geom_errorbar(aes(ymin = lower.CL, 
                    ymax = upper.CL), 
                width = 0.2) +
  geom_text(aes(label = gsub(" ", "", .group)),
            position = position_nudge(x = 0.2)) +
  labs(caption = "Means followed by a common letter are\nnot significantly different according to the Tukey-test")

^{由 reprex package 创建于 2021-06-03 (v2.0.0)}