在带facet_wrap与colour映射的3因子ggplot箱线图添加CLD

一、给分面箱线图添加紧凑字母显示（CLD） #

geom_tukey默认处理全局检验，无法直接适配分面的分组比较，需手动生成每个分面内的CLD结果再与ggplot结合，步骤如下：

1. 生成CLD标注数据 #

假设你的数据框为banana_data，包含变量：DAY（成熟天数）、TREATMENT（灌溉液/树类型组合或单因子）、taste_score（美味度）。用dplyr按DAY分组，在每个组内完成事后检验并生成字母标注：

library(dplyr)
library(emmeans)
library(multcompView)

cld_data <- banana_data %>%
  group_by(DAY) %>%
  do({
    # 拟合组内模型（正态分布用aov，非参数可替换为kruskal适配）
    model <- aov(taste_score ~ TREATMENT, data = .)
    # 计算均值并生成CLD
    emm <- emmeans(model, ~ TREATMENT)
    cld <- multcompLetters(emm$emmean, emm$SE, adjust = "tukey")$Letters
    # 整理成适配ggplot的格式
    data.frame(TREATMENT = names(cld), 
               mean_score = emm$emmean,
               cld = as.character(cld),
               DAY = unique(.$DAY))
  }) %>%
  ungroup()

2. 将CLD添加到箱线图 #

在原有箱线图代码后加入geom_text，把字母定位到对应组的上方：

library(ggplot2)

# 原有箱线图基础
p <- ggplot(banana_data, aes(x = TREATMENT, y = taste_score, colour = TREATMENT)) +
  geom_boxplot() +
  facet_wrap(~ DAY) +
  theme_bw()

# 叠加CLD标注
p + geom_text(data = cld_data, 
              aes(x = TREATMENT, y = mean_score + 0.5, label = cld),
              colour = "black", size = 4, vjust = 0)

注：y = mean_score + 0.5的偏移量可根据美味度的数值范围调整，避免字母遮挡箱线；若数据不符合正态，可将aov替换为kruskal.test后，用emmeans的type = "response"参数适配。

二、适合汇报给同事的3因子非参数检验方法 #

针对美味度这类主观评分（常不符合正态分布），推荐两种简单易解释的方法：

1. 分层双因子非参数检验（按成熟天数拆分） #

将3因子拆分为每个成熟天数下，灌溉液×香蕉树类型的双因子检验，用coin包的independence_test可同时检验主效应与交互效应，结果直观：

library(coin)

# 按成熟天数分层检验
banana_data %>%
  group_by(DAY) %>%
  do({
    test <- independence_test(taste_score ~ irrigation * tree_type, 
                              data = ., 
                              teststat = "quadratic")
    # 整理成易读的p值表格
    p_vals <- data.frame(
      Effect = c("灌溉液", "香蕉树类型", "灌溉液×树类型交互"),
      p_value = c(pvalue(test)[1], pvalue(test)[2], pvalue(test)[3])
    )
    cbind(DAY = unique(.$DAY), p_vals)
  }) %>%
  ungroup()

汇报时可直接展示分层的p值表格，比如：“我们按成熟天数分层检验，结果显示第7天灌溉液对美味度有显著影响（p<0.05），第14天灌溉液与树类型的交互作用显著（p<0.03）”，同事易快速理解。

2. 非参数混合效应模型（全局检验） #

若需要整体结论，用lme4结合Bootstrap做全局3因子检验，适合展示跨天数的整体效应：

library(lme4)
library(boot)

# 拟合混合效应模型，纳入所有因子及交互项
model <- lmer(taste_score ~ irrigation * tree_type * DAY + (1|DAY), data = banana_data)
# 用Bootstrap计算p值
boot_p <- bootMer(model, FUN = function(x) fixef(x), nsim = 1000)
# 提取各效应的p值
p_vals_boot <- sapply(1:length(fixef(model)), function(i) {
  mean(abs(boot_p$t[,i]) > abs(fixef(model)[i]))
})
names(p_vals_boot) <- names(fixef(model))

汇报时聚焦显著效应，比如：“全局检验显示，成熟天数与灌溉液的交互作用对香蕉美味度有显著影响（Bootstrap p=0.03）”，适合需要整体结论的场景。

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