如何执行Kruskal-Wallis检验并为数据框添加紧凑字母显示

核心思路 #

Kruskal-Wallis是ANOVA的非参数替代方案，适配非正态分布数据；事后检验用Dunn检验，再通过multcompView包生成紧凑字母显示（CLD），最后结合中位数（非参数场景下替代均值）用ggplot完成绘图。

所需R包 #

先安装并加载必要工具包：

install.packages(c("FSA", "multcompView", "ggplot2", "dplyr"))
library(FSA)
library(multcompView)
library(ggplot2)
library(dplyr)

示例代码（单分组场景） #

1. 模拟测试数据

set.seed(123)
data <- data.frame(
  group = rep(c("A", "B", "C", "D"), each = 20),
  value = c(rnorm(20, 5), rnorm(20, 7), rnorm(20, 5), rnorm(20, 9))
)

2. 执行Kruskal-Wallis整体检验

kw_test <- kruskal.test(value ~ group, data = data)
# 提取并格式化整体p值
kw_pval <- round(kw_test$p.value, 3)

3. 执行Dunn事后检验（Bonferroni校正）

dunn_result <- dunnTest(value ~ group, data = data, method = "bonferroni")
# 转换结果为数据框方便后续处理
dunn_df <- as.data.frame(dunn_result$res)

4. 生成紧凑字母显示（CLD）

先整理检验对比对和校正后p值，再生成字母标签：

comp_pairs <- dunn_df %>%
  select(Comparison, P.adj) %>%
  mutate(
    group1 = gsub(" - .*", "", Comparison),
    group2 = gsub(".* - ", "", Comparison)
  ) %>%
  select(group1, group2, P.adj)

# 基于校正p值生成字母分组
cld_result <- cldList(P.adj ~ group1 + group2, data = comp_pairs, threshold = 0.05)

5. 合并中位数统计量与字母标签

summary_stats <- data %>%
  group_by(group) %>%
  summarise(
    median_val = median(value),
    .groups = "drop"
  ) %>%
  left_join(cld_result, by = c("group" = "Group"))

6. 绘图（含箱线图、中位数标签、紧凑字母及整体p值）

非参数场景用四分位距（IQR）作为误差线更合理：

ggplot(data, aes(x = group, y = value)) +
  geom_boxplot(fill = "lightblue", alpha = 0.7) +
  # 添加中位数上方的紧凑字母
  geom_text(data = summary_stats, aes(x = group, y = median_val + 1, label = Letter), 
            size = 5, fontface = "bold") +
  # 添加整体检验p值标注
  annotate("text", x = 2.5, y = max(data$value) + 1, 
           label = paste("Kruskal-Wallis p =", kw_pval), size = 4) +
  theme_minimal() +
  labs(x = "分组", y = "数值")

分物种批量处理场景 #

如果需要按物种分组执行检验并绘图，用dplyr的嵌套分组实现批量处理：

# 模拟带物种维度的测试数据
set.seed(456)
data_sp <- data.frame(
  species = rep(c("Sp1", "Sp2"), each = 80),
  group = rep(c("A", "B", "C", "D"), each = 20, times = 2),
  value = c(rnorm(20, 5), rnorm(20, 7), rnorm(20, 5), rnorm(20, 9),
            rnorm(20, 6), rnorm(20, 6), rnorm(20, 8), rnorm(20, 8))
)

# 批量处理每个物种的检验与统计量
processed_data <- data_sp %>%
  group_by(species) %>%
  nest() %>%
  mutate(
    # 批量执行Kruskal-Wallis检验
    kw_test = map(data, ~kruskal.test(value ~ group, data = .x)),
    kw_pval = map_dbl(kw_test, ~round(.x$p.value, 3)),
    # 批量执行Dunn事后检验
    dunn_result = map(data, ~dunnTest(value ~ group, data = .x, method = "bonferroni")),
    dunn_df = map(dunn_result, ~as.data.frame(.x$res)),
    # 批量生成紧凑字母
    cld = map(dunn_df, function(df) {
      comp_pairs <- df %>%
        select(Comparison, P.adj) %>%
        mutate(
          group1 = gsub(" - .*", "", Comparison),
          group2 = gsub(".* - ", "", Comparison)
        ) %>%
        select(group1, group2, P.adj)
      cldList(P.adj ~ group1 + group2, data = comp_pairs, threshold = 0.05)
    }),
    # 合并中位数与字母标签
    summary_stats = map2(data, cld, function(dat, cld_dat) {
      dat %>%
        group_by(group) %>%
        summarise(median_val = median(value), .groups = "drop") %>%
        left_join(cld_dat, by = c("group" = "Group"))
    })
  ) %>%
  select(species, summary_stats, kw_pval) %>%
  unnest(summary_stats)

# 分物种绘图
ggplot(data_sp, aes(x = group, y = value)) +
  geom_boxplot(fill = "lightgreen", alpha = 0.7) +
  geom_text(data = processed_data, aes(x = group, y = median_val + 1, label = Letter),
            size = 5, fontface = "bold") +
  annotate("text", x = 2.5, y = max(data_sp$value) + 1.5,
           label = paste("Kruskal-Wallis p =", kw_pval), size = 4) +
  facet_wrap(~species, scales = "free_y") +
  theme_minimal() +
  labs(x = "分组", y = "数值")

关键说明 #

• 非参数检验场景下，用中位数替代均值更符合统计逻辑，误差线优先选四分位距（IQR）或百分位数范围
• Dunn检验的method参数可切换校正方式，支持holm、bh等多种校正方法
• cldList的threshold参数对应显著性水平，默认0.05，可根据需求调整

内容的提问来源于stack exchange，提问作者A.Benson