首先帮你理清两个核心问题：

1. 你不小心把stat_compare_means()写成了stat_compare_nea()，这是函数名笔误，直接导致统计检验代码无效；
2. 当你给ggboxplot传入多个y变量时，它本质是自动把数据转成长格式并生成分面图，默认的stat_compare_means()不会自动在每个分面里独立计算组间差异，需要额外配置参数适配这个场景。

修正后的完整解决方案 #

下面是针对你的需求调整后的代码，包含了正确的函数调用和分面适配：

library(palmerpenguins)
library(tidyverse)
library(ggplot2)
library(ggpubr)

# 加载并清洗数据
data(package = 'palmerpenguins')
df_clean <- na.omit(penguins)

# 定义你需要的组间比较列表
comp_list <- list(c("Chinstrap", "Adelie"), 
                  c("Chinstrap", "Gentoo"), 
                  c("Adelie", "Gentoo"))

# 生成带p值的多指标箱线图
df_boxplot <- ggboxplot(df_clean, 
                        x = "species", 
                        y = c("bill_length_mm", "bill_depth_mm", "flipper_length_mm", "body_mass_g"), 
                        ylab = "Measurement Value",  # 调整y轴标签，适配多个不同指标
                        xlab = "Species", 
                        color = "species", 
                        fill = "species", 
                        notch = TRUE, 
                        alpha = 0.5, 
                        ggtheme = theme_pubr()) +
  # 关键：为每个分面添加指定组的两两比较p值
  stat_compare_means(comparisons = comp_list, 
                     show.legend = FALSE,
                     label.y = c(60, 22, 240, 6500))  # 给每个分面单独设置p值的垂直位置，适配各指标的数值范围
df_boxplot

额外优化：手动控制分面（可选） #

如果你想更灵活地控制分面样式（比如自由缩放y轴），可以先把数据转成长格式，再手动用facet_wrap实现分面，逻辑会更清晰：

# 将宽格式数据转为长格式
df_long <- df_clean %>% 
  pivot_longer(cols = c(bill_length_mm, bill_depth_mm, flipper_length_mm, body_mass_g),
               names_to = "measurement",
               values_to = "value")

# 生成箱线图
ggboxplot(df_long,
          x = "species",
          y = "value",
          color = "species",
          fill = "species",
          notch = TRUE,
          alpha = 0.5,
          ggtheme = theme_pubr()) +
  facet_wrap(~measurement, scales = "free_y") +  # 每个分面独立缩放y轴
  stat_compare_means(comparisons = comp_list, show.legend = FALSE)

为什么你的原始代码不工作？ #

1. 函数名笔误：stat_compare_nea()是不存在的函数，必须改成正确的stat_compare_means()；
2. 分面适配问题：当ggboxplot处理多个y变量时，会自动生成分面，但默认的stat_compare_means()如果不指定comparisons，只会计算整体ANOVA的p值，而且不会自动在每个分面重复渲染。通过传入comparisons参数，并为每个分面设置合适的label.y，就能让p值正确显示在每个子图中。

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