如何在ggpairs配对散点图中拆分相关系数框展示两种处理的相关值

我使用如下R语言代码基于GGally包的ggpairs函数绘制散点矩阵图，现有代码仅针对全量数据集计算相关系数矩阵，未区分处理组别。我的数据集包含名为Si的分组列，对应2种处理水平：0mM、4mM，下三角区域我已经通过自定义lowerFn函数实现了按Si分组展示散点与线性拟合线。目前我希望将上三角的相关系数展示框拆分为两部分，分别对应两个Si处理组的相关系数，方便组间对比。

library(GGally) 
leg <- grab_legend(ggplot(data=data1, aes(x=NA, y=NA, colour=Si)) +
                     geom_line(size=2))

my_fn <- function(data1, mapping, method="p", use="pairwise", ...){
  
  # grab data
  x <- eval_data_col(data1, mapping$x)
  y <- eval_data_col(data1, mapping$y)
  
  # calculate correlation
  corr <- cor(x, y, method=method, use=use)
  
  # calculate colour based on correlation value
  # Here I have set a correlation of minus one to blue, 
  # zero to white, and one to red 
  # Change this to suit: possibly extend to add as an argument of `my_fn`
  colFn <- colorRampPalette(c("blue", "white", "red"), interpolate ='spline')
  fill <- colFn(100)[findInterval(corr, seq(-1, 1, length=100))]
  
  ggally_cor(data=data1, size=5, digits=2, stars=TRUE, mapping=mapping, ...) + 
    theme_void() +
    theme(panel.background=element_rect(fill=fill))
}

lowerFn <- function(data1, mapping, emap=NULL, method = "lm", ...) {
  # mapping <- c(mapping, emap)
  # class(mapping) = "uneval" # need this to combine the two aes
  # Can use this instead
  mapping <- ggplot2:::new_aes( c(mapping, emap))
  p <- ggplot(data = data1, mapping = mapping) +
    geom_point(data = data1, alpha = 0.8, size = 3, shape = 16) +
    geom_smooth(method = method, ...) +
    theme_gray() # to get the white background and prominent axis
  p
}

ggpairs(
  data1, columns=4:6, legend=leg,
  upper = list(continuous=my_fn),
  lower = list(continuous = 
                 wrap(lowerFn, 
                      method = "lm", 
                      emap=aes(color=Si),
                      fullrange=TRUE, 
                      se=FALSE))) +
  theme(legend.position='top')

你只需要重写上三角的自定义相关系数计算函数，按Si分组分别计算两组的相关系数，再将单元格拆分为上下两个区域分别展示即可，修改后的完整代码如下：

library(GGally)
library(ggplot2)

# 提取图例代码不变
leg <- grab_legend(ggplot(data=data1, aes(x=NA, y=NA, colour=Si)) +
                     geom_line(size=2))

# 重写分组相关系数展示的上三角函数
grouped_cor_fn <- function(data, mapping, method="p", use="pairwise", ...){
  # 定义颜色映射函数不变
  colFn <- colorRampPalette(c("blue", "white", "red"), interpolate ='spline')
  
  # 按Si分组计算两组的相关系数
  cor_res <- by(data, data$Si, function(subdat){
    x <- eval_data_col(subdat, mapping$x)
    y <- eval_data_col(subdat, mapping$y)
    corr <- cor(x, y, method=method, use=use)
    # 计算显著性p值
    p_val <- cor.test(x, y, method=method)$p.value
    # 生成显著性星号
    stars <- ifelse(p_val < 0.001, "***",
                    ifelse(p_val < 0.01, "**",
                           ifelse(p_val < 0.05, "*", "")))
    return(list(corr=round(corr,2), stars=stars, fill=colFn(100)[findInterval(corr, seq(-1, 1, length=100))]))
  })
  
  # 提取两组结果
  cor_0mM <- cor_res[["0mM"]]
  cor_4mM <- cor_res[["4mM"]]
  
  # 绘制拆分的相关系数面板
  p <- ggplot() +
    # 上半部分放4mM结果
    annotate("rect", xmin=0, xmax=1, ymin=0.5, ymax=1, fill=cor_4mM$fill) +
    annotate("text", x=0.5, y=0.75, label=paste0("4mM: ", cor_4mM$corr, cor_4mM$stars), size=4) +
    # 下半部分放0mM结果
    annotate("rect", xmin=0, xmax=1, ymin=0, ymax=0.5, fill=cor_0mM$fill) +
    annotate("text", x=0.5, y=0.25, label=paste0("0mM: ", cor_0mM$corr, cor_0mM$stars), size=4) +
    # 设置坐标系和主题
    coord_cartesian(expand=FALSE) +
    theme_void()
  return(p)
}

# 下三角函数不变
lowerFn <- function(data1, mapping, emap=NULL, method = "lm", ...) {
  mapping <- ggplot2:::new_aes( c(mapping, emap))
  p <- ggplot(data = data1, mapping = mapping) +
    geom_point(alpha = 0.8, size = 3, shape = 16) +
    geom_smooth(method = method, ...) +
    theme_gray()
  p
}

# 调用ggpairs时替换upper的自定义函数即可
ggpairs(
  data1, columns=4:6, legend=leg,
  upper = list(continuous=grouped_cor_fn),
  lower = list(continuous = 
                 wrap(lowerFn, 
                      method = "lm", 
                      emap=aes(color=Si),
                      fullrange=TRUE, 
                      se=FALSE))) +
  theme(legend.position='top')

• 你可以根据需求调整文本大小、颜色映射规则、显著性判定阈值等参数
• 如果Si的水平顺序和示例不同，对应调整cor_res提取的分组名即可

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