在成对比较后用 p 值制作堆图
Making a heapmap with p values after pairwise comparisons
我正在 运行 对分组因素进行多重成对比较,我想用热图表示生成的 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,用于另一个热图,这样就更容易看出意义了。
有人可以帮我解决这个问题吗?我找了好几个小时都没找到办法。
实现您想要的结果的一个选项是通过 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 值。
这是我的数据集的一个小例子,但实际上要复杂得多。 我的因素是地点(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,用于另一个热图,这样就更容易看出意义了。
有人可以帮我解决这个问题吗?我找了好几个小时都没找到办法。
实现您想要的结果的一个选项是通过 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).