在R中使用group_by时为dunn_test添加紧凑字母显示的方法

我们已经通过rstatix包的dunn_test完成了非参数检验，现在需要给分组结果添加紧凑字母显示（即cld列），实现类似方差分析后的字母分组效果。

1. 创建示例数据 #

# 生成物种分组
species <- rep(c("Oak", "Elm", "Ash"), each = 10)

# 生成处理组
dose_1 <- rep("Ctrl", 30)
dose_2 <- rep("L", 30)

# 生成结果数据
result_1 <- c(
  runif(10, 9, 12), runif(10, 14, 16), runif(10, 6, 8),
  runif(10, 2, 5), runif(10, 1, 4), runif(10, 2, 4)
)

# 合并为数据框
data <- data.frame(
  species = rep(species, 2),
  treatment = c(dose_1, dose_2),
  result = result_1
)

2. 现有Dunn检验代码 #

library(tidyverse)
library(rstatix)

data %>% 
  group_by(species) %>% 
  dunn_test(result ~ treatment, p.adjust.method = "holm")

3. 期望输出格式 #

我们需要得到包含均值、标准误、置信区间以及紧凑字母列cld的结果：

species treatment emmean    SE    df lower.CL upper.CL cld  
  <chr>   <fct>      <dbl> <dbl> <dbl>    <dbl>    <dbl> <chr>
1 Oak     Ctrl       10.6  0.249    18    10.1     11.1  A    
2 Oak     L           3.57 0.249    18     3.04     4.09 B    
3 Elm     Ctrl       15.2  0.252    18    14.7     15.8  A    
4 Elm     L           2.93 0.252    18     2.40     3.46 B    
5 Ash     Ctrl        7.15 0.173    18     6.79     7.51 A    
6 Ash     L           2.96 0.173    18     2.59     3.32 B

4. 实现方案 #

结合emmeans计算组统计量，用multcompView生成紧凑字母，具体代码如下：

# 加载所需包
library(tidyverse)
library(rstatix)
library(multcompView)

# 1. 按物种分组计算处理组的统计量：均值、标准误、置信区间
group_stats <- data %>%
  group_by(species, treatment) %>%
  summarise(
    emmean = mean(result),
    SE = sd(result)/sqrt(n()),
    df = n() - 1,
    lower.CL = emmean - qt(0.975, df)*SE,
    upper.CL = emmean + qt(0.975, df)*SE,
    .groups = "drop"
  )

# 2. 提取Dunn检验的调整后p值
dunn_pvals <- data %>%
  group_by(species) %>%
  dunn_test(result ~ treatment, p.adjust.method = "holm") %>%
  select(species, group1, group2, p.adj)

# 3. 为每个物种生成紧凑字母分组
cld_groups <- map(unique(data$species), function(sp) {
  # 筛选当前物种的检验结果
  sp_pvals <- dunn_pvals %>% filter(species == sp)
  # 创建比较矩阵
  comp_matrix <- matrix(sp_pvals$p.adj, nrow = 1, 
                        dimnames = list(paste(sp_pvals$group1, sp_pvals$group2, sep = "-"), NULL))
  # 生成字母分组
  letters <- multcompLetters(comp_matrix, threshold = 0.05)$Letters
  # 转换为数据框
  tibble(
    species = sp,
    treatment = names(letters),
    cld = as.character(letters)
  )
}) %>% bind_rows()

# 4. 合并统计量与字母分组
final_output <- group_stats %>%
  left_join(cld_groups, by = c("species", "treatment")) %>%
  arrange(species, desc(emmean))

# 查看最终结果
final_output

补充说明 #

• 如果需要用非参数中心趋势（如中位数）代替均值，只需将mean(result)替换为median(result)，并调整置信区间的计算逻辑（比如用百分位数法）。
• multcompLetters函数会根据调整后的p值分配字母：p值>0.05的分组共享相同字母，反之字母不同。

内容的提问来源于stack exchange，提问作者UseR10085