如何结合apply()函数使用rstatix包的pairwise_wilcox_test？

问题背景 #

现有如下数据集data1，需针对其中的微生物丰度列（Firmicutes、Proteobacteria、Verrucomicrobiota）与分组列cls进行两两Wilcoxon检验：

data1 <- structure(list(Firmicutes = c(9916885L, 9923240L, 9897778L, 9887923L, 9832122L, 9717375L, 9820485L, 9805007L, 9676859L, 9636967L, 9581184L, 9908333L, 9624107L, 9899086L, 9827215L),
                        Proteobacteria = c(83115.37, 76759.73, 102222.14, 112077.44, 167423.55, 235007.98, 150719.87, 54276.39, 320811.45, 363032.82, 400989.93, 87339.68, 147003.76, 67276.26, 165133.37),
                        Verrucomicrobiota = c(0, 0, 0, 0, 454.1326, 47616.9546, 28794.7347, 140716.5721, 2329.3241, 0, 17825.6204, 4327.6418, 228889.5762, 33638.1295, 7651.654),
                        cls = c(1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L)),
                   class = "data.frame", row.names = c("1", "10", "13", "16", "19", "11", "14", "17", "2", "20", "12", "15", "18", "21", "24"))

针对单个列执行检验时运行正常：

WIL <- rstatix::pairwise_wilcox_test(Firmicutes ~ cls, data=data1, exact = TRUE, p.adjust.method="bonferroni")

输出结果：

# A tibble: 3 × 9
  .y.        group1 group2    n1    n2 statistic     p p.adj p.adj.signif
* <chr>      <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl> <chr>       
1 Firmicutes 1      2         5     5        15 0.443 1     ns          
2 Firmicutes 1      3         5     5        13 0.977 1     ns          
3 Firmicutes 2      3         5     5         7 0.101 0.303 ns

尝试用apply()批量处理时出现错误：

WIL <- apply(as.matrix(data1),2, function(x){rstatix::pairwise_wilcox_test(x ~ cls, data=data1,exact = TRUE, p.adjust.method="bonferroni")})

错误信息：

ℹ In index: 1.
    ℹ With name: V1.
    Caused by error in `pull()`:
    ! Can't extract columns that don't exist.
    ✖ Column `x` doesn't exist.
    Run `rlang::last_trace()` to see where the error occurred.
    Called from: signal_abort(cnd, .file)

错误原因 #

1. 公式解析逻辑问题：pairwise_wilcox_test的公式参数是在data指定的数据框环境中查找变量的。你传递的x是apply遍历得到的列数值向量，而非数据框中的列名，公式x ~ cls会去data1中查找名为x的列，自然找不到。
2. 不必要的类型转换：将data1转为矩阵后，丢失了数据框的列名与结构信息，且会把分组列cls也纳入遍历范围，这是不需要的。

解决方案 #

方法1：使用lapply遍历列名（base R） #

先提取需要检验的列名（排除分组列cls），再动态构建公式执行检验：

# 获取需要检验的列名
cols_to_test <- setdiff(names(data1), "cls")

# 批量执行检验
WIL <- lapply(cols_to_test, function(col_name) {
  # 动态构建公式：列名 ~ cls
  formula_str <- paste(col_name, "~ cls")
  rstatix::pairwise_wilcox_test(data = data1, formula = as.formula(formula_str),
                                exact = TRUE, p.adjust.method = "bonferroni")
})

# 为结果列表命名，方便查看
names(WIL) <- cols_to_test

方法2：使用purrr::map（tidyverse风格） #

如果你习惯tidyverse工具链，可以用purrr::map实现更简洁的批量处理：

library(purrr)

cols_to_test <- setdiff(names(data1), "cls")

WIL <- map(cols_to_test, ~ rstatix::pairwise_wilcox_test(
  data = data1,
  formula = as.formula(paste(.x, "~ cls")),
  exact = TRUE,
  p.adjust.method = "bonferroni"
))

names(WIL) <- cols_to_test

执行后，WIL是一个列表，每个元素对应一个微生物列的检验结果，可通过WIL$Firmicutes、WIL$Proteobacteria等单独查看。

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