如何用R语言基于TukeyHSD绘制事后检验p值三角热图

我正在分析30种职业（每种职业约50名受试者）之间的血压差异。通过单因素ANOVA检验确认不同职业间血压存在显著差异，已完成TukeyHSD事后分析，代码如下：

AOV <- aov(Blood_pressure ~ Job, data = myData)
Tukey <- TukeyHSD(AOV)

TukeyHSD()返回的结果包含difference（组间血压差异值）和post-hoc p value（事后检验校正p值）的数据框。我需要在R中基于该结果绘制三角热图，要求按p值区间（<0.05、<0.01、<0.001、NS（无显著性））显示不同颜色。TukeyHSD结果的结构如下（保留5位小数）：

structure(c(0.74955, 0.35787, ...))

以下是实现三角热图的完整步骤：

1. 整理TukeyHSD结果为整洁数据框 #

首先将TukeyHSD返回的列表转换为易于处理的数据框，并拆分职业对比组：

# 加载所需包
library(dplyr)
library(tidyr)

# 提取TukeyHSD的职业对比结果
tukey_raw <- as.data.frame(Tukey$Job)
# 添加对比组名称列，并拆分为两个职业
tukey_raw$comparison <- rownames(tukey_raw)
tukey_df <- separate(tukey_raw, comparison, into = c("Job1", "Job2"), sep = "-")
# 保留核心列：两个职业、校正p值、差异值
tukey_df <- tukey_df %>%
  select(Job1, Job2, `p adj`, diff) %>%
  rename(p_value = `p adj`, difference = diff)

2. 标记p值的显著性区间 #

根据p值给每个对比组打上对应的显著性标签：

tukey_df <- tukey_df %>%
  mutate(
    sig_level = case_when(
      p_value < 0.001 ~ "*** (<0.001)",
      p_value < 0.01 ~ "** (<0.01)",
      p_value < 0.05 ~ "* (<0.05)",
      TRUE ~ "NS"
    )
  )
# 统一职业的因子水平，保证热图中顺序一致
all_jobs <- unique(sort(c(tukey_df$Job1, tukey_df$Job2)))
tukey_df$Job1 <- factor(tukey_df$Job1, levels = all_jobs)
tukey_df$Job2 <- factor(tukey_df$Job2, levels = all_jobs)

3. 绘制三角热图 #

使用ggplot2绘制仅显示上三角的热图，避免重复的组间对比：

library(ggplot2)

ggplot(tukey_df, aes(x = Job1, y = Job2)) +
  # 仅绘制上三角区域（Job1的排序在前，避免重复对比）
  geom_tile(aes(fill = sig_level), data = filter(tukey_df, as.numeric(Job1) < as.numeric(Job2))) +
  # 自定义显著性区间的颜色
  scale_fill_manual(
    values = c(
      "*** (<0.001)" = "#d73027", # 深红色
      "** (<0.01)" = "#fc8d59",   # 橙色
      "* (<0.05)" = "#fee08b",    # 浅黄色
      "NS" = "#e0f3f8"            # 浅蓝灰色
    ),
    name = "显著性水平"
  ) +
  # 调整主题，适配职业名称显示
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 8),
    axis.text.y = element_text(size = 8),
    axis.title = element_blank(),
    legend.position = "bottom"
  )

可选优化 #

• 显示血压差异值：如果需要在热图块中显示组间血压差异，可添加geom_text层：

  ggplot(...) +
    # 原geom_tile代码...
    geom_text(
      aes(label = round(difference, 2)),
      data = filter(tukey_df, as.numeric(Job1) < as.numeric(Job2)),
      size = 3
    )
  

• 调整职业名称大小：根据职业名称长度，修改axis.text.x/y中的size参数
• 安装依赖包：若未安装所需包，先运行：

  install.packages(c("dplyr", "tidyr", "ggplot2"))
  

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