R语言实现二维分箱热图数据多元正态性与中心位置检验

以下实现完全适配tidy数据格式，可直接配合dplyr::group_by完成批量热图检验，针对你给出的8*8模拟数据格式（对应热图如下）可直接运行：

依赖包安装与加载 #

首次运行请先安装所需依赖，加载后即可直接调用函数：

# install.packages(c("tidyverse", "MVN"))
library(tidyverse)
library(MVN)

核心检验函数 #

函数输入为单张热图的整洁数据（需包含x/y/intensity三列），输出为单行结果表，包含两个检验的经-log10转换的p值：

• center_test_log10p：分布中心是否等于指定点的检验结果，值越大越拒绝「围绕指定点中心化」的原假设
• mvn_test_log10p：分布是否服从多元正态的检验结果，值越大越拒绝「服从多元正态分布」的原假设

calc_heatmap_tests <- function(df, center = c(3.5, 3.5)) {
  # 过滤强度为0的无效像素
  valid_df <- df %>% filter(intensity > 0)
  n_total <- sum(valid_df$intensity)
  
  # 有效样本量不足时返回空值避免流程中断
  if(n_total < 3) {
    return(tibble(
      center_test_log10p = NA_real_,
      mvn_test_log10p = NA_real_
    ))
  }
  
  # 计算加权均值、加权协方差矩阵（权重为intensity计数）
  w <- valid_df$intensity
  mat <- as.matrix(valid_df[, c("x", "y")])
  mu_hat <- colSums(mat * w) / n_total
  mat_centered <- sweep(mat, 2, mu_hat)
  sigma_hat <- crossprod(mat_centered * sqrt(w)) / (n_total - 1)
  
  # 1. 中心化Hotelling T²检验
  diff_mu <- mu_hat - center
  t2 <- n_total * t(diff_mu) %*% solve(sigma_hat) %*% diff_mu
  f_stat <- t2 * (n_total - 2) / (2 * (n_total - 1))
  p_center <- pf(f_stat, df1 = 2, df2 = n_total - 2, lower.tail = FALSE)
  log10p_center <- -log10(p_center)
  
  # 2. 多元正态性Mardia加权检验
  mvn_res <- mvn(
    data = mat,
    weights = w,
    mvnTest = "mardia",
    desc = FALSE,
    showOutliers = FALSE,
    showNewData = FALSE
  )
  # 取偏度、峰度检验p值的最小值做保守判断
  p_mvn <- min(mvn_res$multivariateNormality$p)
  log10p_mvn <- -log10(p_mvn)
  
  # 返回标准化结果
  tibble(
    center_test_log10p = as.numeric(log10p_center),
    mvn_test_log10p = as.numeric(log10p_mvn)
  )
}

批量调用方式 #

如果所有热图数据存储在同一个数据框中，用heatmap_id类的列区分不同热图，直接搭配group_by+group_modify即可批量计算，完全适配整洁数据流：

# 构造2张测试热图
set.seed(42)
multi_heatmap_df <- bind_rows(
  tibble(
    heatmap_id = "random_map",
    x = rep(0:7, each = 8),
    y = rep(0:7, 8),
    intensity = sample(0:10, 64, replace = TRUE)
  ),
  tibble(
    heatmap_id = "skew_map",
    x = rep(0:7, each = 8),
    y = rep(0:7, 8),
    # 模拟单侧偏斜的非正态热图
    intensity = dbinom(x, 7, 0.2) * dbinom(y, 7, 0.8) * 1000
  )
)

# 批量执行检验
result <- multi_heatmap_df %>%
  group_by(heatmap_id) %>%
  group_modify(~calc_heatmap_tests(.x, center = c(3.5, 3.5))) %>%
  ungroup()

print(result)

使用说明 #

• 若intensity为连续权重而非计数，代码无需修改可直接运行
• 若需要更稳健的多元正态检验，可将mvn函数的mvnTest参数改为"energy"或"royston"，其中energy检验精度更高但速度稍慢，适合热图总量不大的场景
• 检验返回的-log10(p)可直接用于后续筛选，通常阈值设为2（对应p<0.01）即可判定存在显著偏斜/显著偏离正态分布

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