如何将Kruskal-Wallis检验及两两比较结果存入DataFrame

需求背景 #

已实现对多分组（B1/B2/B3）的循环统计分析，包括计算各年份中位数、执行Kruskal-Wallis检验、BH校正的pairwise.wilcox.test两两比较，现需将所有统计结果整合为结构化DataFrame。

测试数据 #

scaledDaily <- structure(list(loc = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L), levels = c("B1", "B2", "B3"), class = "factor"), 
    mAODscale = c(0.561948848, 0.851883432, 0.293151829, 0.683807808, 
    0.864472706, 0.380303934, 0.324501709, 0.765022304, 0.900772901, 
    0.204731715, 0.877104175, 0.56367162, 0.206528162, 0.353350116, 
    0.219628257, 0.840723901, 0.716389918, 0.569798858, 0.707583441, 
    0.120064246, 0.275325307, 0.438391155, 0.308969668, 0.350156436, 
    0.886955315, 0.693416677, 0.710065022), year = structure(c(1L, 
    1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 
    3L, 3L, 1L, 1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L), levels = c("1", 
    "2", "3"), class = "factor")), row.names = c(NA, -27L), class = "data.frame")

修改后的代码（整合结果为DataFrame） #

library(dplyr)
library(tidyr)

# 初始化空列表存储每个分组的结果
result_list <- list()

names <- c('B1','B2','B3')

for (name in names){
  mereLoc <- scaledDaily %>%
    filter(loc == name) %>%
    mutate(year = as.numeric(year))
  
  # 1. 整理各年份中位数为DataFrame
  median_df <- tapply(mereLoc$mAODscale, mereLoc$year, median, na.rm=T) %>%
    as.data.frame() %>%
    rename(median = ".") %>%
    mutate(year = rownames(.), loc = name) %>%
    relocate(loc, year, median)
  
  # 2. 提取Kruskal-Wallis检验关键结果
  krusk <- kruskal.test(mAODscale~year, data=mereLoc)
  krusk_df <- data.frame(
    loc = name,
    kruskal_stat = krusk$statistic,
    kruskal_p = krusk$p.value,
    df = krusk$parameter
  )
  
  # 3. 整理两两Wilcox检验结果（BH校正）
  pairs <- pairwise.wilcox.test(mereLoc$mAODscale, mereLoc$year, p.adjust.method='BH')
  pairs_df <- pairs$p.value %>%
    as.data.frame() %>%
    mutate(year1 = rownames(.)) %>%
    pivot_longer(cols = -year1, names_to = "year2", values_to = "wilcox_p_bh") %>%
    filter(!is.na(wilcox_p_bh)) %>%
    mutate(loc = name) %>%
    relocate(loc, year1, year2, wilcox_p_bh)
  
  # 合并当前分组的所有结果
  loc_result <- median_df %>%
    left_join(krusk_df, by = "loc") %>%
    left_join(pairs_df, by = c("loc", "year" = "year1"))
  
  # 添加到结果列表
  result_list[[name]] <- loc_result
}

# 合并所有分组结果为最终DataFrame
final_result <- bind_rows(result_list)

# 查看最终结果
print(final_result)

代码说明 #

• 中位数整理：将tapply输出的向量转为结构化DataFrame，补充loc和year标识，保证数据格式统一。
• Kruskal-Wallis结果提取：直接从检验对象中提取统计量、p值和自由度，封装为单行DataFrame便于合并。
• 两两比较结果整理：将pairwise检验的矩阵结果通过pivot_longer转为长格式，过滤重复对比的NA值，补充分组标识。
• 结果合并：通过left_join将中位数、Kruskal检验结果、两两比较结果按分组和年份关联，最后合并所有分组数据。

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