如何在rstatix::wilcox_test中使用变量作为公式参数？

尝试通过循环绘制箱线图，并添加rstatix::wilcox_test()计算的p值，将公式所需列名存入变量后触发报错：

Error in `pull()`:
! Can't extract columns that don't exist.
✖ Column `rs_id` doesn't exist.

原代码 #

rs_ids <- c('rs2160669', 'rs964184')
phenos <- c('HBA1C', 'TG')
for (rs_id in rs_ids){
  for(phen in phenos){
    bxp <- ggboxplot(d.k, x = rs_id, y = phen)
    stat.test <- d.k %>%  wilcox_test(.data[[`phen`]] ~ .data[[`rs_id`]]) %>%  add_significance()
    stat.test <- stat.test %>% add_xy_position(x = .data[[rs_id]]) %>% mutate(myformatted.p = ifelse(p < 0.001, 'p<0.001', paste0('p=',signif(p, digits = 2))))
    
    bxp <-  bxp +   geom_jitter(position=position_jitter(0.2), alpha = 0.5)   +  theme_bw() + theme(legend.position="none", axis.text=element_text(size=12), axis.title=element_text(size=14,face="bold")) + xlab(rs_id) + ylab(phen)
    plot <- bxp + stat_pvalue_manual(stat.test, label = "{myformatted.p}", tip.length = 0.01, step.increase = 0.1)
  }
}

示例数据 #

structure(list(HBA1C = c(4.5, 5.9, 8.02, 5.5, 5.4, 6, 7.1, 9.9, 
5.58, 5, 7.6, 8, 5.5, 6.4, 6.8, 9, 6.2, 4.9, 5.7, 6.2, 6.3, 5.7, 
5.9, 6.9, 7.3, 5.7, 7.2, 10.4, 5.4, 5.3, 4.8, 5.5, 7.1, 6.2, 
7, 7.4, 11.4, 5.4, 5.5, 5.3, 5.5, 5.9, 5.5, 5.79, 4.8, 6.16, 
5.63, 7.41, 7.68, 5.82), TG = c(0.83, 0.95, 2.48, 0.78, 0.48, 
1.16, 1.45, 1.32, 1.03, 3.77, 3.87, 2.34, 1.98, 2.59, 2.92, 2.71, 
6.88, 1.25, 3.68, 1.15, 2.33, 3.37, 0.61, 1.02, 1.63, 1.32, 1.21, 
1.35, 0.85, 3.97, 4.04, 3.63, 2.62, 1.46, 2.33, 2.46, 1.09, 1.46, 
2.77, 3.1, 3.13, 2.55, 1.91, 0.97, 0.87, 1.46, 1.45, 1.15, 2.61, 
2.15), rs2160669 = c(3, 2, 3, 3, 3, 3, 3, 3, 3, 2, 3, 2, 3, 2, 
1, 2, 3, 3, 3, 3, 2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 
2, 3, 3, 3, 2, 3, NA, 3, 3, 3, 3, 3, 3, 2, 2), rs964184 = c(1, 
2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 3, 2, 1, 1, 1, 1, 2, 2, 
1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, NA, 
2, 1, 3, 1, 1, 2, 2, 2)), row.names = c(NA, -50L), class = c("tbl_df", 
"tbl", "data.frame"))

报错核心原因是：在wilcox_test和add_xy_position中错误地将变量名用反引号包裹，导致R把rs_id、phen当成字面列名而非变量。以下是修正方案：

修正要点 #

1. wilcox_test公式写法：直接使用.data[[phen]] ~ .data[[rs_id]]（无需反引号），或用as.formula()动态生成公式，后者更符合R公式语法习惯。
2. add_xy_position参数：直接传入字符串变量rs_id即可，该函数支持接收列名字符串，不需要.data语法。

修正后完整代码 #

library(ggplot2)
library(rstatix)
library(ggpubr)

# 加载示例数据
d.k <- structure(list(HBA1C = c(4.5, 5.9, 8.02, 5.5, 5.4, 6, 7.1, 9.9, 
5.58, 5, 7.6, 8, 5.5, 6.4, 6.8, 9, 6.2, 4.9, 5.7, 6.2, 6.3, 5.7, 
5.9, 6.9, 7.3, 5.7, 7.2, 10.4, 5.4, 5.3, 4.8, 5.5, 7.1, 6.2, 
7, 7.4, 11.4, 5.4, 5.5, 5.3, 5.5, 5.9, 5.5, 5.79, 4.8, 6.16, 
5.63, 7.41, 7.68, 5.82), TG = c(0.83, 0.95, 2.48, 0.78, 0.48, 
1.16, 1.45, 1.32, 1.03, 3.77, 3.87, 2.34, 1.98, 2.59, 2.92, 2.71, 
6.88, 1.25, 3.68, 1.15, 2.33, 3.37, 0.61, 1.02, 1.63, 1.32, 1.21, 
1.35, 0.85, 3.97, 4.04, 3.63, 2.62, 1.46, 2.33, 2.46, 1.09, 1.46, 
2.77, 3.1, 3.13, 2.55, 1.91, 0.97, 0.87, 1.46, 1.45, 1.15, 2.61, 
2.15), rs2160669 = c(3, 2, 3, 3, 3, 3, 3, 3, 3, 2, 3, 2, 3, 2, 
1, 2, 3, 3, 3, 3, 2, 2, 3, 3, 3, 3, 3, 2, 3, 3, 3, 3, 3, 3, 3, 
2, 3, 3, 3, 2, 3, NA, 3, 3, 3, 3, 3, 3, 2, 2), rs964184 = c(1, 
2, 1, 2, 1, 2, 1, 1, 1, 2, 2, 2, 1, 2, 3, 2, 1, 1, 1, 1, 2, 2, 
1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, NA, 
2, 1, 3, 1, 1, 2, 2, 2)), row.names = c(NA, -50L), class = c("tbl_df", 
"tbl", "data.frame"))

rs_ids <- c('rs2160669', 'rs964184')
phenos <- c('HBA1C', 'TG')

# 存储所有生成的图，避免循环中被覆盖
plots_list <- list()

for (rs_id in rs_ids){
  for(phen in phenos){
    # 绘制基础箱线图
    bxp <- ggboxplot(d.k, x = rs_id, y = phen) +
      geom_jitter(position=position_jitter(0.2), alpha = 0.5) +
      theme_bw() +
      theme(legend.position="none", 
            axis.text=element_text(size=12), 
            axis.title=element_text(size=14,face="bold")) +
      xlab(rs_id) + ylab(phen)
    
    # 方法1：用.data[[变量名]]构建公式
    stat.test <- d.k %>%  
      wilcox_test(.data[[phen]] ~ .data[[rs_id]]) %>%  
      add_significance()
    
    # 方法2：用as.formula动态生成公式（可选，效果一致）
    # stat.test <- d.k %>%  
    #   wilcox_test(as.formula(paste(phen, "~", rs_id))) %>%  
    #   add_significance()
    
    # 添加位置信息并格式化p值
    stat.test <- stat.test %>% 
      add_xy_position(x = rs_id) %>% 
      mutate(myformatted.p = ifelse(p < 0.001, 'p<0.001', paste0('p=',signif(p, digits = 2))))
    
    # 将p值添加到图中
    plot <- bxp + stat_pvalue_manual(stat.test, label = "{myformatted.p}", tip.length = 0.01, step.increase = 0.1)
    
    # 存入列表方便后续查看/导出
    plots_list[[paste(rs_id, phen, sep = "_")]] <- plot
  }
}

# 示例：查看其中一个图
# print(plots_list[["rs2160669_HBA1C"]])

补充说明 #

• 新增plots_list存储所有生成的图，避免循环中后续图覆盖之前的结果。
• 两种公式写法可任选，as.formula更贴近传统R公式语法，适合习惯此类写法的用户。

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