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

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,