使用 ggplot2 将显着性级别添加到矩阵相关热图

Posted 2023-02-18

【中文标题】使用 ggplot2 将显着性级别添加到矩阵相关热图

【英文标题】：Significance level added to matrix correlation heatmap using ggplot2

【发布时间】：2012-08-25 04:11:03

【问题描述】：

我想知道如何在矩阵相关热图中添加另一层重要且需要的复杂性，例如除了 R2 值（-1 到 1）之外的显着性水平星的方式之后的 p 值？ 在这个问题中，不打算将显着性水平星号或 p 值作为文本放在矩阵的每个正方形上，而是在矩阵的每个正方形上以显着性水平的开箱即用图形表示形式显示这一点。我认为只有那些享受创新思维祝福的人才能赢得掌声来解开这种解决方案，以便有最好的方式来表示我们的“半真半真矩阵相关热图”中增加的复杂性组件。我用谷歌搜索了很多，但从未见过合适的，或者我会说一种“眼睛友好”的方式来表示显着性水平加上反映 R 系数的标准色调。 可重现的数据集在这里找到：http://learnr.wordpress.com/2010/01/26/ggplot2-quick-heatmap-plotting/ R代码请在下面找到：

library(ggplot2)
library(plyr) # might be not needed here anyway it is a must-have package I think in R 
library(reshape2) # to "melt" your dataset
library (scales) # it has a "rescale" function which is needed in heatmaps 
library(RColorBrewer) # for convenience of heatmap colors, it reflects your mood sometimes
nba <- read.csv("http://datasets.flowingdata.com/ppg2008.csv")
nba <- as.data.frame(cor(nba[2:ncol(nba)])) # convert the matrix correlations to a dataframe 
nba.m <- data.frame(row=rownames(nba),nba) # create a column called "row"
rownames(nba) <- NULL #get rid of row names
nba <- melt(nba)
nba.m$value<-cut(nba.m$value,breaks=c(-1,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1),include.lowest=TRUE,label=c("(-0.75,-1)","(-0.5,-0.75)","(-0.25,-0.5)","(0,-0.25)","(0,0.25)","(0.25,0.5)","(0.5,0.75)","(0.75,1)")) # this can be customized to put the correlations in categories using the "cut" function with appropriate labels to show them in the legend, this column now would be discrete and not continuous
nba.m$row <- factor(nba.m$row, levels=rev(unique(as.character(nba.m$variable)))) # reorder the "row" column which would be used as the x axis in the plot after converting it to a factor and ordered now
#now plotting
ggplot(nba.m, aes(row, variable)) +
geom_tile(aes(fill=value),colour="black") +
scale_fill_brewer(palette = "RdYlGn",name="Correlation")  # here comes the RColorBrewer package, now if you ask me why did you choose this palette colour I would say look at your battery charge indicator of your mobile for example your shaver, won't be red when gets low? and back to green when charged? This was the inspiration to choose this colour set.

矩阵相关热图应如下所示：

增强解决方案的提示和想法： - 此代码可能有助于了解从该网站获取的显着性水平星：http://ohiodata.blogspot.de/2012/06/correlation-tables-in-r-flagged-with.html R代码：

mystars <- ifelse(p < .001, "***", ifelse(p < .01, "** ", ifelse(p < .05, "* ", " "))) # so 4 categories  

- 可以像 alpha 美学一样将显着性级别作为颜色强度添加到每个正方形，但我认为这不容易解释和捕捉 - 另一个想法是有 4 个不同大小的正方形对应于星星，当然给最小的给不重要的，如果最高的星星增加到全尺寸的正方形 - 在这些重要的正方形内包含一个圆圈的另一个想法，圆圈的线的粗细对应于一种颜色的重要性级别（剩下的 3 个类别） - 与上述相同，但固定线条粗细，同时为剩余的 3 个重要级别提供 3 种颜色 - 可能你想出了更好的主意，谁知道呢？

【问题讨论】：

你的代码启发我用ggplot2重写arm::corrplot函数：rpubs.com/briatte/ggcorr

效果很好！您能否扩展此功能以使那些不显着的相关性（例如

感谢您的反馈。我对这里（和其他地方）使用 $p$-values 持怀疑态度，但我会尝试找出一些方法来标记无关紧要的系数。

上面引用的函数现在是 GGally 包的一部分，包的维护者进行了更正和添加。

(-1, -0,75) 颜色在哪里？？使用 c(-1,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1)，我们应该有 8 个区间和 8 种颜色，而不是 7...

【参考方案1】：

这只是对最终解决方案的增强尝试，我在这里绘制了星星作为解决方案的指标，但正如我所说的，目的是找到一个比星星更能说明问题的图形解决方案。我只是使用 geom_point 和 alpha 来指示显着性水平，但是 NA（也包括非显着值）会像三颗星的显着性水平一样出现问题，如何解决这个问题？我认为在使用多种颜色时使用一种颜色可能对眼睛更友好，并且避免用眼睛解决的许多细节给情节带来负担。提前致谢。 这是我第一次尝试的情节：

或者可能会更好？！

我认为到目前为止最好的是下面的那个，直到你想出更好的东西！

根据要求，以下代码用于最后一张热图：

# Function to get the probability into a whole matrix not half, here is Spearman you can change it to Kendall or Pearson
cor.prob.all <- function (X, dfr = nrow(X) - 2) 
R <- cor(X, use="pairwise.complete.obs",method="spearman")
r2 <- R^2
Fstat <- r2 * dfr/(1 - r2)
R<- 1 - pf(Fstat, 1, dfr)
R[row(R) == col(R)] <- NA
R

# Change matrices to dataframes
nbar<- as.data.frame(cor(nba[2:ncol(nba)]),method="spearman") # to a dataframe for r^2
nbap<- as.data.frame(cor.prob.all(nba[2:ncol(nba)])) # to a dataframe for p values
# Reset rownames
nbar <- data.frame(row=rownames(nbar),nbar) # create a column called "row" 
rownames(nbar) <- NULL
nbap <- data.frame(row=rownames(nbap),nbap) # create a column called "row" 
rownames(nbap) <- NULL
# Melt
nbar.m <- melt(nbar)
nbap.m <- melt(nbap)
# Classify (you can classify differently for nbar and for nbap also)         
nbar.m$value2<-cut(nbar.m$value,breaks=c(-1,-0.75,-0.5,-0.25,0,0.25,0.5,0.75,1),include.lowest=TRUE, label=c("(-0.75,-1)","(-0.5,-0.75)","(-0.25,-0.5)","(0,-0.25)","(0,0.25)","(0.25,0.5)","(0.5,0.75)","(0.75,1)")) # the label for the legend
nbap.m$value2<-cut(nbap.m$value,breaks=c(-Inf, 0.001, 0.01, 0.05),label=c("***", "** ", "*  ")) 
nbar.m<-cbind.data.frame(nbar.m,nbap.m$value,nbap.m$value2) # adding the p value and its cut to the first dataset of R coefficients
names(nbar.m)[5]<-paste("valuep") # change the column names of the dataframe 
names(nbar.m)[6]<-paste("signif.")
nbar.m$row <- factor(nbar.m$row, levels=rev(unique(as.character(nbar.m$variable)))) # reorder the variable factor
# Plotting the matrix correlation heatmap
# Set options for a blank panel
po.nopanel <-list(opts(panel.background=theme_blank(),panel.grid.minor=theme_blank(),panel.grid.major=theme_blank()))
pa<-ggplot(nbar.m, aes(row, variable)) +
geom_tile(aes(fill=value2),colour="white") +
scale_fill_brewer(palette = "RdYlGn",name="Correlation")+ # RColorBrewer package
opts(axis.text.x=theme_text(angle=-90))+
po.nopanel
pa # check the first plot
# Adding the significance level stars using geom_text 
pp<- pa +
geom_text(aes(label=signif.),size=2,na.rm=TRUE) # you can play with the size
# Workaround for the alpha aesthetics if it is good to represent significance level, the same workaround can be applied for size aesthetics in ggplot2 as well. Applying the alpha aesthetics to show significance is a little bit problematic, because we want the alpha to be low while the p value is high, and vice verse which can't be done without a workaround
nbar.m$signif.<-rescale(as.numeric(nbar.m$signif.),to=c(0.1,0.9)) # I tried to use to=c(0.1,0.9) argument as you might expect, but to avoid problems with the next step of reciprocal values when dividing over one, this is needed for the alpha aesthetics as a workaround
nbar.m$signif.<-as.factor(0.09/nbar.m$signif.) # the alpha now behaves as wanted  except for the NAs values stil show as if with three stars level, how to fix that?
# Adding the alpha aesthetics in geom_point in a shape of squares (you can improve here)
pp<- pa +
geom_point(data=nbar.m,aes(alpha=signif.),shape=22,size=5,colour="darkgreen",na.rm=TRUE,legend=FALSE) # you can remove this step, the result of this step is seen in one of the layers in the above green heatmap, the shape used is 22 which is again a square but the size you can play with it accordingly  

我希望这可以成为向前迈出的一步！请注意： - 一些人建议对 R^2 进行不同的分类或切割，好吧，我们当然可以这样做，但我们仍然希望以图形方式向观众展示显着性水平，而不是用星级水平来打扰眼睛。我们能否在原则上实现这一目标？ - 一些人建议以不同的方式削减 p 值，好的，这可以是在显示 3 个显着性水平失败后的选择，而不会打扰眼睛。那么最好在没有水平的情况下显示显着/非显着 - 对于 ggplot2 中的上述解决方法，您可能会想出一个更好的主意，以获得 alpha 和尺寸美学，希望很快收到您的来信！ - 这个问题还没有回答，等待一个创新的解决方案！ - 有趣的是，“corrplot”包做到了！我通过这个包得出了下面的这个图表，PS：交叉的方块不是显着的，signif = 0.05的水平。但是我们怎么能把它翻译成ggplot2，可以吗？！

-或者你可以做圆圈并隐藏那些不重要的？如何在 ggplot2 中做到这一点？！

【讨论】：

很棒的图，但是人们应该意识到，在这种情况下，p 值可能并不意味着他们对相关系数热图的期望（或任何东西）。如果你有这么多可能的相关性（$n^2-n/2$），甚至 >.99 的 p 值实际上开始变得有点可能。在这种情况下过度依赖 p 值可能被认为是 p-hacking 的另一个名称。这个xkcd 解释得很漂亮。

opts() 很久以前就被弃用了。使用 theme() 选项。

【参考方案2】：

library("corrplot")
nba <- as.matrix(read.csv("https://raw.githubusercontent.com/Shicheng-Guo/Shicheng-Guo.Github.io/master/data/ppg2008.csv")[-1])
res1 <- cor.mtest(nba, conf.level = .95)
par(mfrow=c(2,2))

# correlation and P-value
corrplot(cor(nba), p.mat = res1$p, insig = "label_sig",sig.level = c(.001, .01, .05), pch.cex = 0.8, pch.col = "white",tl.cex=0.8)

# correlation and hclust
corrplot(cor(nba), method = "shade", outline = T, addgrid.col = "darkgray", order="hclust", 
         mar = c(4,0,4,0), addrect = 4, rect.col = "black", rect.lwd = 5,cl.pos = "b", tl.col = "indianred4", 
         tl.cex = 0.8, cl.cex = 0.8)

【参考方案3】：

要沿估计的相关系数表示显着性，您可以改变着色量 - 使用 alpha 或仅填充每个图块的子集：

# install.packages("fdrtool")
# install.packages("data.table")
library(ggplot2)
library(data.table)

#download dataset
nba <- as.matrix(read.csv("http://datasets.flowingdata.com/ppg2008.csv")[-1])
m <- ncol(nba)
# compute corellation and p.values for all combinations of columns
dt <- CJ(i=seq_len(m), j=seq_len(m))[i<j]
dt[, c("p.value"):=(cor.test(nba[,i],nba[,j])$p.value), by=.(i,j)]
dt[, c("corr"):=(cor(nba[,i],nba[,j])), by=.(i,j)]

# estimate local false discovery rate
dt[,lfdr:=fdrtool::fdrtool(p.value, statistic="pvalue")$lfdr]

dt <- rbind(dt, setnames(copy(dt),c("i","j"),c("j","i")), data.table(i=seq_len(m),j=seq_len(m), corr=1, p.value=0, lfdr=0))


#use alpha
ggplot(dt, aes(x=i,y=j, fill=corr, alpha=1-lfdr)) + 
  geom_tile()+
  scale_fill_distiller(palette = "RdYlGn", direction=1, limits=c(-1,1),name="Correlation") +
  scale_x_continuous("variable", breaks = seq_len(m), labels = colnames(nba)) +
  scale_y_continuous("variable", breaks = seq_len(m), labels = colnames(nba), trans="reverse") +
  coord_fixed() +
  theme(axis.text.x=element_text(angle=90, vjust=0.5),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
  )

#use area
ggplot(dt, aes(x=i,y=j, fill=corr,  height=sqrt(1-lfdr),  width=sqrt(1-lfdr))) + 
  geom_tile()+
  scale_fill_distiller(palette = "RdYlGn", direction=1, limits=c(-1,1),name="Correlation") +
  scale_color_distiller(palette = "RdYlGn", direction=1, limits=c(-1,1),name="Correlation") +
  scale_x_continuous("variable", breaks = seq_len(m), labels = colnames(nba)) +
  scale_y_continuous("variable", breaks = seq_len(m), labels = colnames(nba), trans="reverse") +
  coord_fixed() +
  theme(axis.text.x=element_text(angle=90, vjust=0.5),
        panel.background=element_blank(),
        panel.grid.minor=element_blank(),
        panel.grid.major=element_blank(),
  )

这里的关键是 p.values 的缩放：为了获得仅在相关区域中显示出较大变化的易于解释的值，我使用 @ 提供的局部错误发现 (lfdr) 的上限估计987654328@ 代替。 即，图块的 alpha 值可能小于或等于该相关性不同于 0 的概率。