如何使用 R 语言绘制相关系数热力图？交互式的更好？

之前有个小伙伴问过这样一个问题：

老师，问下怎么画不同变量间的相关性矩阵图，然后在图上每一个方格里面加上相关系数和对应的P值啊?

我们今天就一起来计算下相关系数矩阵并用图展示下，在展示手段上，我将分别使用 ggplot2 和 highcharter。

首先我们需要计算相关系数矩阵和 p 值矩阵，Hmisc 包的 rcorr 函数可以实现，下面的代码中我将使用 mtcars 数据进行演示：

library(tidyverse)
mtcars %>%
  as_tibble() -> df

Hmisc::rcorr(as.matrix(df), type = "pearson") -> corrlist

得到的 corrlist 是个 list，里面有三个数据框，我们需要的是 r 和 P：

names(corrlist)
#> [1] "r" "n" "P"

首先提取相关系数矩阵：

corrlist$r %>%
  as_tibble() %>%
  mutate(v = colnames(.)) %>%
  select(v, everything()) %>%
  pivot_longer(2:12) -> corrdf

#> # A tibble: 121 x 3
#>    v     name   value
#>    <chr> <chr>  <dbl>
#>  1 mpg   mpg    1
#>  2 mpg   cyl   -0.852
#>  3 mpg   disp  -0.848
#>  4 mpg   hp    -0.776
#>  5 mpg   drat   0.681
#>  6 mpg   wt    -0.868
#>  7 mpg   qsec   0.419
#>  8 mpg   vs     0.664
#>  9 mpg   am     0.600
#> 10 mpg   gear   0.480
#> # … with 111 more rows

然后是提取 p 值矩阵：

corrlist$P %>%
  as_tibble() %>%
  mutate(v = colnames(.)) %>%
  select(v, everything()) %>%
  pivot_longer(2:12) %>%
  mutate(label = case_when(
    is.na(value) ~ "",
    value <= 0.001 ~ "\n***",
    between(value, 0.001, 0.01) ~ "\n**",
    between(value, 0.01, 0.1) ~ "\n*",
    T ~ ""
  )) -> pdf
pdf

#> # A tibble: 121 x 4
#>    v     name      value label
#>    <chr> <chr>     <dbl> <chr>
#>  1 mpg   mpg   NA        ""
#>  2 mpg   cyl    6.11e-10 "\n***"
#>  3 mpg   disp   9.38e-10 "\n***"
#>  4 mpg   hp     1.79e- 7 "\n***"
#>  5 mpg   drat   1.78e- 5 "\n***"
#>  6 mpg   wt     1.29e-10 "\n***"
#>  7 mpg   qsec   1.71e- 2 "\n*"
#>  8 mpg   vs     3.42e- 5 "\n***"
#>  9 mpg   am     2.85e- 4 "\n***"
#> 10 mpg   gear   5.40e- 3 "\n**"
#> # … with 111 more rows

这里我生成了一列新的变量 label 表示显著性，这里的 \n 是用于接下来绘图的时候换行的。

然后我们将这两个数据框连接起来：

corrdf %>%
  left_join(pdf, by = c("v", "name")) %>%
  rename(corr = value.x, p = value.y) %>%
  mutate(corr = round(corr, 2)) -> corrdf

corrdf

#> # A tibble: 121 x 5
#>    v     name   corr         p label
#>    <chr> <chr> <dbl>     <dbl> <chr>
#>  1 mpg   mpg    1    NA        ""
#>  2 mpg   cyl   -0.85  6.11e-10 "\n***"
#>  3 mpg   disp  -0.85  9.38e-10 "\n***"
#>  4 mpg   hp    -0.78  1.79e- 7 "\n***"
#>  5 mpg   drat   0.68  1.78e- 5 "\n***"
#>  6 mpg   wt    -0.87  1.29e-10 "\n***"
#>  7 mpg   qsec   0.42  1.71e- 2 "\n*"
#>  8 mpg   vs     0.66  3.42e- 5 "\n***"
#>  9 mpg   am     0.6   2.85e- 4 "\n***"
#> 10 mpg   gear   0.48  5.40e- 3 "\n**"
#> # … with 111 more rows

最后就可以绘图了，我们需要使用的是 geom_tile 图层：

corrdf %>%
  mutate(v = forcats::fct_reorder(v, corr),
         name = forcats::fct_reorder(name, corr)) %>%
  ggplot(aes(x = v, y = name)) +
  geom_tile(aes(fill = corr)) +
  geom_text(aes(label = paste(corr, label, sep = "")),
            family = cnfont, size = 2.8) +
  ggthemes::scale_fill_gradient2_tableau("Red-Blue Diverging") +
  labs(x = "\n注：* <0.1 **<0.01, ***<0.01",
       y = "",
       title = "mtcars 数据集相关系数矩阵",
       subtitle = "数据来源：datasets 包",
       caption = "绘图：微信公众号 RStata")

ggsave("pic1.png", width = 8, height = 6, dpi = 300)

图1

另外还可以使用 highcharter 绘制一幅炫酷的交互式热力图：

library(highcharter)
corrdf %>%
  mutate(v = forcats::fct_reorder(v, corr),
         name = forcats::fct_reorder(name, corr)) %>%
  hchart("heatmap",
         hcaes(x = v, y = name, value = corr),
         dataLabels = list(
           enabled = TRUE,
           format = "{point.corr}<br>{point.label}",
           align = "center",
           style = list(fontFamily = "STSong")
         )) %>%
  hc_colorAxis(dataClasses = JS('
      [{from: -1, to: -0.5, color: "#EFF3FF"},
      {from: -0.5, to: 0, color: "#BDD7E7"},
      {from: 0, to: 0.5, color: "#6BAED6"},
      {from: 0.5, to: 1, color: "#2171B5"}]'),
               labels = list(style = list(fontFamily = "STSong"))) %>%
  hc_xAxis(title = list(text = "注：* <0.1 **<0.01, ***<0.01")) %>%
  hc_yAxis(title = list(text = JS('null'))) %>%
  hc_legend(title = list(text = "相关系数",
                         style = list(fontFamily = "STSong")),
            itemStyle = list(fontFamily = "STSong")) %>%
  hc_title(text = "mtcars 数据集相关系数矩阵",
           style = list(fontFamily = "STSong")) %>%
  hc_subtitle(text = "数据来源：datasets 包",
              style = list(fontFamily = "STSong")) %>%
  hc_add_theme(hc_theme_google()) %>%
  hc_credits(enabled = TRUE,
             text = "微信公众号 RStata") %>%
  hc_tooltip(headerFormat = "",
             pointFormat = "<b>{point.v} 和 {point.name}变量的相关性为：</b>{point.corr}",
             borderRadius = 5,
             style = list(fontFamily = "STSong")) %>%
  hc_chart(margin = c(70, 70, 100, 70)) %>%
  hc_exporting(enabled = T)

图2

点击这里跳转到 RStata 短书平台获取附件：如何使用 R 语言绘制相关系数热力图？交互式的更好？