在成对比较后用 p 值制作堆图

Question

我正在运行对分组因素进行多重成对比较，我想用热图表示生成的 p 值。

这是我的数据集的一个小例子，但实际上要复杂得多。我的因素是地点（2 个级别）和生物体（3 个级别）。

site = c("SITE1","SITE1","SITE1","SITE1","SITE1","SITE1",
          "SITE1","SITE1","SITE1","SITE2","SITE2","SITE2",
          "SITE2","SITE2","SITE2","SITE2","SITE2","SITE2")
organism = c("Insects","Insects","Insects","Mammals","Mammals",
          "Mammals","Reptiles","Reptiles","Reptiles","Insects",
          "Insects","Insects","Mammals","Mammals","Mammals",
          "Reptiles","Reptiles","Reptiles")
variable = c(5,6,7,12,13,14,1,2,3,7,8,9,22,24,25,11,12,14)


data = data.frame(site, organism, variable)
head(data)

  site organism variable
1 SITE1  Insects        5
2 SITE1  Insects        6
3 SITE1  Insects        7
4 SITE1  Mammals       12
5 SITE1  Mammals       13
6 SITE1  Mammals       14

这就是我在每个单独站点计算生物对之间的成对测试的方式：

data %>%
  group_by(site) %>%
  t_test(variable ~ organism)

   site      .y.  group1   group2 n1 n2  statistic       df        p    p.adj p.adj.signif
1 SITE1 variable Insects  Mammals  3  3  -8.573214 4.000000 0.001000 0.002000           **
2 SITE1 variable Insects Reptiles  3  3   4.898979 4.000000 0.008000 0.008000           **
3 SITE1 variable Mammals Reptiles  3  3  13.472194 4.000000 0.000176 0.000528          ***
4 SITE2 variable Insects  Mammals  3  3 -14.862705 3.448276 0.000300 0.000900          ***
5 SITE2 variable Insects Reptiles  3  3  -4.110961 3.448276 0.020000 0.020000            *
6 SITE2 variable Mammals Reptiles  3  3   9.086882 4.000000 0.000813 0.002000           **

我想以一种可以制作类似于此的热图的方式组织我的数据：

我认为包 emmeans 中的函数 pwpm 做了类似的事情，但它只适用于 emmeans。我找不到其他任何东西。

此外，如果也可以将p.values细分为p<0.05，p<0.01，p<0.001，用于另一个热图，这样就更容易看出意义了。

有人可以帮我解决这个问题吗？我找了好几个小时都没找到办法。

Answer 1

实现您想要的结果的一个选项是通过 ggplot2 和分面。由于您的数据已经采用整洁的数据格式，因此可以通过 geom_tile + facet_wrap 获得基本的热图。剩下的就是样式，比如通过 geom_text 添加 p 值，通过 scale_fill_gradient 设置颜色或使用例如一些额外的数据整理tidyr::complete 添加“缺失”组类别。

library(dplyr)
library(rstatix)
library(ggplot2)

d <- data %>%
  group_by(site) %>%
  t_test(variable ~ organism) %>% 
  tidyr::complete(group1 = unique(data$organism), group2 = unique(data$organism), site = unique(data$site))

ggplot(d, aes(group2, rev(group1), fill = p)) +
  geom_tile() +
  geom_text(aes(label = scales::number(p, accuracy = 1e-6))) +
  scale_fill_gradient(low = "red", high = "green", na.value = NA) +
  facet_wrap(~site, ncol = 1)
#> Warning: Removed 12 rows containing missing values (geom_text).

在成对比较后用 p 值制作堆图

Making a heapmap with p values after pairwise comparisons

r

ggplot2

matrix

heatmap

pairwise