求可展示双向ANOVA原始均值、SE及大小写显著性字母的R包

针对你的需求（原始均值、原始标准误、大小写字母分别标注列/行方向的显著性），可以通过dplyr分组统计原始值 + agricolae包做分组显著性检验 + tidyr整理表格结构的组合方案实现，具体步骤如下：

1. 调整模型（移除交互项） #

由于交互效应不显著，重新拟合主效应模型：

fit <- aov(log(response) ~ treatments + species, data = sample_data)

2. 计算原始均值与标准误 #

直接从原始数据分组计算，避免使用模型估算值：

library(dplyr)
# 分组统计原始均值、标准误
summary_stats <- sample_data %>%
  group_by(treatments, species) %>%
  summarise(
    mean_raw = mean(response),          # 原始响应变量的均值
    se_raw = sd(response)/sqrt(n()),    # 原始标准误
    .groups = "drop"
  )

3. 生成双向显著性字母标注 #

使用agricolae包的HSD.test分别对列（species）、行（treatments）方向做Tukey检验，生成大小写字母：

library(agricolae)

# 列向比较（每个treatments下的species差异）：生成大写字母
species_labels <- sample_data %>%
  group_by(treatments) %>%
  group_modify(~ {
    # 对当前treatments组内的species做ANOVA+Tukey检验
    species_fit <- aov(log(response) ~ species, data = .x)
    hsd_result <- HSD.test(species_fit, "species", group = TRUE)
    # 提取分组字母并转为大写
    tibble(species = rownames(hsd_result$groups), col_letter = toupper(hsd_result$groups$groups))
  })

# 横向比较（每个species下的treatments差异）：生成小写字母
treatment_labels <- sample_data %>%
  group_by(species) %>%
  group_modify(~ {
    # 对当前species组内的treatments做ANOVA+Tukey检验
    treat_fit <- aov(log(response) ~ treatments, data = .x)
    hsd_result <- HSD.test(treat_fit, "treatments", group = TRUE)
    # 提取分组字母并转为小写
    tibble(treatments = rownames(hsd_result$groups), row_letter = tolower(hsd_result$groups$groups))
  })

4. 合并数据并整理为目标表格 #

将统计量与字母标注合并，转换为行是treatments、列是species的宽表：

library(tidyr)

final_table <- summary_stats %>%
  # 合并列向、横向字母标注
  left_join(species_labels, by = c("treatments", "species")) %>%
  left_join(treatment_labels, by = c("treatments", "species")) %>%
  # 格式化显示内容：均值±SE + 大小写字母
  mutate(
    mean_se = paste0(round(mean_raw, 2), " ± ", round(se_raw, 2)),
    cell_content = paste0(mean_se, " ", col_letter, row_letter)
  ) %>%
  # 转换为宽表结构
  pivot_wider(
    id_cols = treatments,
    names_from = species,
    values_from = cell_content
  )

5. 美化表格（可选） #

用kableExtra生成更美观的可导出表格：

library(kableExtra)
final_table %>%
  kbl(caption = "双向ANOVA结果：原始均值±SE与显著性标注（大写=列向species差异，小写=横向treatments差异）") %>%
  kable_classic(full_width = FALSE)

方案说明 #

• 原始统计量：直接从sample_data分组计算，完全符合你"非估算值"的要求
• 显著性标注：通过分组做Tukey检验，分别生成大小写字母，精准对应列/行方向的比较需求
• 灵活性：可以根据需求调整小数位数、字母格式或表格样式

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