correlation_heatmap

Introduction

两个数值变量看散点图就够了。但 10 个变量有 45 对关系，30 个变量有 435 对——逐一画散点做不到。

相关性热图把所有两两相关系数压进一个矩阵，用颜色表达方向和强度。它是探索入口，不是结论：先用它发现哪里可能有结构，再回到具体变量对上画散点、查离群点、做模型。

读的时候按这个顺序：

1. 先看方向。 一端是负相关、另一端是正相关，中间色接近 0。
2. 再看深浅。 越深，相关系数绝对值越大。
3. 看色块。 聚类排序后相似的变量被排到一起，会出现连续的块状结构。
4. 看星号。 它说的是"在这个样本量下相关性是否显著偏离 0"，不是效应大小。强相关看颜色，显著性看星号，两件事。
5. 回到散点图。 对最要紧的那几对变量一定要再画散点确认——热图只压缩，不解释。

Example Data

相关性热图要的输入是：3 个以上数值变量，同一批样本，一行一个样本。变量必须是数值——分类变量编码成 1/2/3 也能算出系数，但那个数没有意义。

mtcars 是 32 辆车的 11 项指标。这里取其中 10 个连续或近似连续的变量：mpg 和 wt、disp、hp、cyl 明显负相关——这不能读成"重量导致油耗降低"，只能说这组变量和燃油效率有系统性的反向变化。

# 公开地址，和数据集页上「下载 CSV」给的是同一个 —— 不写仓库相对路径：
# 那个目录不进仓库，读者 clone 下来也没有这个文件，这段代码就跑不了
d <- read.csv("https://assets.evanzhou.org/tessera/csv/mtcars.csv")
vars <- c("mpg", "cyl", "disp", "hp", "drat", "wt", "qsec", "vs", "am", "gear")
df <- d[, vars]

cor_mat <- cor(
  df,
  use = "pairwise.complete.obs",  # 有缺失时逐对使用完整观测，不整行删掉
  method = "pearson"              # 线性相关；变量分布很偏时改 "spearman"
)

星号来自 p 值，得单独算一张同样大小的矩阵。这个函数只依赖 base R：

cor_pmat <- function(x, method = "pearson") {
  mat <- as.matrix(x)
  n <- ncol(mat)
  p <- matrix(NA_real_, n, n)
  # 行列名必须跟着带上：下面两种画法都靠名字对齐，靠位置迟早错位
  colnames(p) <- rownames(p) <- colnames(mat)
  for (i in seq_len(n)) {
    for (j in seq_len(n)) {
      p[i, j] <- cor.test(mat[, i], mat[, j], method = method)$p.value
    }
  }
  p
}

p_mat <- cor_pmat(df)

Palettes

相关系数有明确方向，所以必须用发散配色：一端负相关、另一端正相关，而中性的那一档必须正好落在 0 上。锚点偏了，读者会把"没关系"读成"弱正相关"——这是配色能造成的最实质的一种误读。

用的是 walter_white（冷蓝 → 冷白 → 深橄榄）。换成红蓝也行，关键不是挑哪两个色，而是读者能不能一眼分清负、零、正三个位置。

library(biopalette)

# 五色发散色板里取两端和正中间三个锚点：低 / 中 / 高
cor_cols <- get_palette("walter_white", type = "diverging")[c(1, 3, 5)]

Recipe

1 · corrplot

一个函数出图，聚类排序、图例、星号都在里面。改颜色改图例不用翻文档。

#| fig: corrplot
#| fig-width: 8
#| fig-height: 6
library(corrplot)

corrplot(
  corr        = cor_mat,
  method      = "color",      # 用色块填格；也可以是 circle / square / number
  type        = "full",       # 画满；upper / lower 只画一半，省掉重复信息
  order       = "hclust",     # 按层次聚类重排变量 —— 块状结构就是这么来的
  hclust.method = "complete", # 聚类的连接方式，影响分块边界怎么切
  col         = colorRampPalette(cor_cols)(200),
  col.lim     = c(-1, 1),     # 固定色标范围，否则色深随数据变、两张图没法比
  addgrid.col = "white",      # 格线颜色，白线让色块之间断开
  cl.pos      = "r",          # 色标（colour legend）放右边
  cl.ratio    = 0.15,         # 色标占整图宽度的比例
  tl.pos      = "l",          # 变量名（text label）只放左边，不重复放上面
  tl.col      = "black",
  tl.srt      = 0,            # 标签旋转角度，0 = 水平
  p.mat       = p_mat,        # 传入 p 值矩阵，星号才画得出来
  sig.level   = c(0.001, 0.01, 0.05),
  insig       = "label_sig",  # 不显著的格子留白；改成 "blank" 会整格涂掉
  pch.cex     = 0.9,          # 星号大小
  pch.col     = "black",
  title       = "Correlation Matrix of mtcars",
  mar         = c(0, 0, 1.8, 0)  # 上边留白，否则标题会被裁掉（base graphics 的老问题）
)

2 · ggcorrplot

返回的是 ggplot 对象，所以能继续叠图层，字体和主题也能和同一套报告里的其它图对齐。

星号这里自己叠一层 geom_text，而不是用包内置的选项——分级、颜色、大小都能自己控制。

#| fig: ggcorrplot
#| fig-width: 8
#| fig-height: 6
library(ggcorrplot)
library(ggplot2)

gg <- ggcorrplot(
  corr          = cor_mat,
  method        = "square",   # 方块；"circle" 会用面积额外编码强度
  type          = "full",
  colors        = cor_cols,   # 低 / 中 / 高，三个锚点
  ggtheme       = ggplot2::theme_void,  # 从空主题起步，自己加要的部分
  hc.order      = TRUE,       # 同 corrplot 的 order = "hclust"
  hc.method     = "complete",
  outline.color = "grey85"    # 格子描边；设 NA 可以完全去掉
) +
  labs(title = "Correlation Matrix of mtcars") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    # 显式给白底：theme_void 是透明的，导出 PNG 时会变成透明背景
    panel.background = element_rect(fill = "white", colour = NA),
    plot.background  = element_rect(fill = "white", colour = NA)
  )

# 星号单独一层：先把 p 值矩阵摊成长表，再按聚类后的顺序对齐
star_mat <- ifelse(p_mat <= 0.001, "***",
             ifelse(p_mat <= 0.01, "**",
              ifelse(p_mat <= 0.05, "*", "")))
diag(star_mat) <- ""   # 对角线是变量和自己，永远显著，标了没意义

star_df <- as.data.frame(as.table(star_mat), stringsAsFactors = FALSE)
names(star_df) <- c("Var1", "Var2", "stars")
star_df <- subset(star_df, nzchar(stars))

# 关键一步：因子水平必须取 ggcorrplot 聚类后的顺序，否则星号会错位到别的格子 ——
# 而错位之后图看起来完全正常
star_df$Var1 <- factor(star_df$Var1, levels = levels(gg$data$Var1))
star_df$Var2 <- factor(star_df$Var2, levels = levels(gg$data$Var2))

p <- gg + geom_text(
  data        = star_df,
  aes(x = Var2, y = Var1, label = stars),
  inherit.aes = FALSE,   # 不继承底图的 fill 映射，否则星号会被上色
  size        = 3.8
)

p

导出。 页面上这两张是 144 dpi，够看。要投稿或进幻灯片得自己导一遍，而两种画法导出方式根本不同——这也是 base graphics 和 ggplot2 最实际的一处区别。

base graphics 的图不是对象：corrplot() 一执行，图就画在当前设备上了，没有东西可以先存起来再交给别的设备。所以三种格式只能各开一次设备、把同一段绘图代码重跑一遍。

#| eval: false
library(ragg)

# 包成函数只是为了不把 corrplot(...) 那一长串抄三遍
draw <- function() {
  corrplot(corr = cor_mat, method = "color", order = "hclust", ...)  # 同上，原样
}

agg_png("cor_mtcars.png", width = 8, height = 6, units = "in",
        res = 300, background = "white")
draw()
dev.off()   # 忘了这句，文件是空的或者被下一张图接着往里画

# 期刊投稿通常要 TIFF，300 dpi 起，LZW 是无损压缩
agg_tiff("cor_mtcars.tiff", width = 8, height = 6, units = "in",
         res = 300, compression = "lzw", background = "white")
draw()
dev.off()

ggplot 是对象，导出就是把对象交给设备，不用重画：

#| eval: false
ggsave("ggcor_mtcars.png", p, width = 8, height = 6, dpi = 300,
       device = ragg::agg_png)   # 指定 ragg：默认设备在 Windows 上字体渲染很糙

ggsave("ggcor_mtcars.tiff", p, width = 8, height = 6, dpi = 300,
       device = function(...) ragg::agg_tiff(..., compression = "lzw"))

ggsave() 不需要 background——前面在 theme() 里显式给的白底就是它的背景。这两件事是同一个：theme_void 不补白底的话，这里导出来就是透明的。

别用 JPEG 存热图。 热图是大块纯色加硬边，JPEG 的有损压缩正好在这种边缘上出振铃，色块交界会糊出一圈杂色——省下的那点体积不值得。JPEG 适合照片，不适合色块图。

Constraints

• 相关不是因果。 一起变化不说明谁导致谁。
• Pearson 只看线性。 一条明显的曲线关系可以给出接近 0 的相关系数——热图上那一格看起来就像"没关系"。
• 离群点和小样本都能让系数跳。 一个极端点就可能显著改变系数和 p 值；样本量小时两者都不稳。
• 检验多了星号就多。 30 个变量是 435 次检验，星号需要谨慎解读。
• 聚类只是描述。 它帮你读结构，不等于发现了机制。

两个 recipe 怎么比

1. 只是想看一眼有没有结构 → recipe 1。参数都在一处。
2. 这张图要进一套统一风格的报告，或者要在图上叠自定义标注 → recipe 2。它是 ggplot 对象；recipe 1 的图层系统是封闭的，加东西只能靠 base graphics 的 text() 硬算坐标。

两张图表达的是同一组系数，选哪个不影响结论，只影响后面好不好改。