如何在R中对trees数据集子集批量执行Kruskal-Wallis检验并生成结果矩阵？

可以用tidyverse生态工具（dplyr+purrr）配合broom包，实现全自动的子集划分、批量检验和结果整理，无需手动编写每个子集的代码。

步骤1：加载依赖包 #

library(tidyverse)
library(broom)

步骤2：确保分组变量为因子类型（可选但推荐） #

如果plot.type、demog、guild.type还不是因子类型，先做转换：

trees <- trees %>% 
  mutate(across(c(plot.type, demog, guild.type), as.factor))

步骤3：批量执行检验并整理结果 #

# 按plot.type和demog划分子集，批量运行检验并整理成整洁表格
kw_results <- trees %>%
  # 按两个因子的所有组合自动分组
  group_by(plot.type, demog) %>%
  # 嵌套每个分组的数据集，方便批量处理
  nest() %>%
  # 对每个子集的count和prop，以guild.type为分组执行Kruskal-Wallis检验
  mutate(
    kw_count = map(data, ~kruskal.test(count ~ guild.type, data = .x) %>% tidy()),
    kw_prop = map(data, ~kruskal.test(prop ~ guild.type, data = .x) %>% tidy())
  ) %>%
  # 展开嵌套的检验结果，合并为统一表格
  unnest(c(kw_count, kw_prop), names_sep = "_") %>%
  # 筛选并重命名列，让结果更直观
  select(plot.type, demog,
         count_stat = kw_count_statistic, 
         count_p = kw_count_p.value,
         prop_stat = kw_prop_statistic, 
         prop_p = kw_prop_p.value)

步骤4：转换为矩阵格式 #

如果需要将结果转为指定的矩阵格式（行是plot.type+demog组合，列是统计量和p值）：

kw_matrix <- kw_results %>%
  # 将plot.type和demog合并为行名
  unite(row_name, plot.type, demog, sep = "_") %>%
  column_to_rownames(var = "row_name") %>%
  as.matrix()

示例输出（矩阵格式） #

count_stat count_p prop_stat prop_p
one_seed        4.231     0.120    5.672     0.059
one_sap         6.894     0.032    3.110     0.211
one_adult       2.015     0.366    4.553     0.103
ten_seed        7.220     0.027    6.104     0.047
non_seed        1.879     0.391    3.881     0.144
...

异常处理（可选） #

如果某些子集里guild.type的类别不足2个（会导致检验报错），可以用purrr::possibly跳过错误，返回NA值：

# 包装检验函数，遇到错误返回NA
safe_kw <- possibly(function(df, var) {
  kruskal.test(reformulate("guild.type", response = var), data = df) %>% tidy()
}, otherwise = tibble(statistic = NA, p.value = NA))

# 重新运行批量检验
kw_results <- trees %>%
  group_by(plot.type, demog) %>%
  nest() %>%
  mutate(
    kw_count = map(data, safe_kw, var = "count"),
    kw_prop = map(data, safe_kw, var = "prop")
  ) %>%
  unnest(c(kw_count, kw_prop), names_sep = "_") %>%
  select(plot.type, demog, count_stat = kw_count_statistic, count_p = kw_count_p.value,
         prop_stat = kw_prop_statistic, prop_p = kw_prop_p.value)

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