ComBat 的 Location-Scale 模型：从原理到调参直觉

做批次校正的人，十有八九绕不开 ComBat。但说实话，很多人只把它当黑盒调参，参数调来调去，最后一问“你到底在调整什么？”，答不上来。我个人习惯是，用工具前先搞懂它背后的数学模型，不然你连报错都读不懂。

ComBat 的核心假设很简单：每个批次的生产数据，可以看作是从一个“共性分布”里加上批次特有的位置偏移（Location）和尺度伸缩（Scale）得来的。用统计语言讲，这就是一个 Location-Scale 模型。你在做的 Batch Effect 校正，本质上就是在估计并消除这两个扰动量。

Location-Scale 模型的数学骨架

我们先把公式摆出来，你会发现它其实不复杂。假设你有 \(i\) 个特征（比如基因表达矩阵里的基因）， \(j\) 个样本，\(k\) 个批次。原始数据 \(Y_{ijk}\) 被建模为：

Y_{ijk} = α_i + β_i * X_j + γ_{ik} + δ_{ik} * ε_{ijk}

其中：
  α_i   —— 特征的全局基线均值
  β_i   —— 协变量（如年龄、性别）对特征的影响系数
  X_j   —— 样本 j 的协变量向量
  γ_{ik} —— 批次 k 对特征 i 的**位置偏移**（均值偏移）
  δ_{ik} —— 批次 k 对特征 i 的**尺度因子**（方差缩放）
  ε_{ijk} —— 残差，服从 N(0, σ_i²)

注意看，这里的 γ_ik 和 δ_ik 才是我们真正要校正的目标。 在校正过的数据 Y_ijk^adj 里，它们会被减去或除掉：

Y_{ijk}^{adj} = (Y_{ijk} - γ_{ik}) / δ_{ik} + γ_{ik}^{global}

γ_ik^global 是全局均值，目的是让校正后数据的整体中心不跑偏。嗯，这里要注意—— 很多使用者会忘记加回这个全局均值，导致校正后的数据整体平移了一个常数。

贝叶斯先验：不让人多嘴的经验法则

直接对每个批次每个特征单独估计 γ 和 δ 是不行的。样本少时，估计值方差大，校正完反而引入噪声。ComBat 干的聪明事是：把所有特征的 γ 和 δ 放到一个层次模型里，用经验贝叶斯（Empirical Bayes）方法共享信息。

简单说，它们假设：

• γ_ik 来自一个正态分布，均值为 0，方差为 τ_k²（先验方差）。
• δ_ik 来自一个逆伽马分布，形状参数 λ_k，尺度参数 θ_k。

这些先验参数 τ_k²、λ_k、θ_k 是从所有特征的数据中估计出来的。你想想看，这就相当于把所有特征的“批次效应强弱”聚在一起看，然后给每个单独的估计“拉一把”。

用代码表达这个先验估计过程（简化版）：

# 伪代码：对每个批次 k 计算先验参数
for k in batches:
    # 先算每个特征的原始偏移和方差
    gamma_hat[k, i] = mean(Y[i, j in batch k]) - global_mean[i]
    delta_hat[k, i] = var(Y[i, j in batch k]) / pool_var[i]

    # 然后从所有特征中求先验
    tau_k_sq = var(gamma_hat[k, :])   # 位置偏移的方差
    lambda_k, theta_k = fit_inv_gamma(delta_hat[k, :])  # 尺度因子的先验

后验估计与校正：贝叶斯把力往中间“拉”

有了先验，下一步就是计算后验均值 E[γ_ik | data] 和 E[δ_ik | data]。 对于正态-正态共轭模型，后验公式是：

γ_{ik}^{post} = ( tau_k² / ( tau_k² + sigma_i² )) * γ_{ik}^{hat}

也就是：后验均值 = 先验方差 / (先验方差 + 残差方差) * 原始估计

这个比值叫做“收缩因子”。当某个特征的原始估计方差（sigma_i²）很大时， 收缩因子接近 0，后验值被拉向先验均值 0，相当于告诉系统“这个特征的单批样本太少， 别信它的偏移，用全局信息补”。反之，如果特征本身很稳定（sigma_i² 小）， 就基本保留原始估计。这个机制，说白了就是贝叶斯在替你做权衡。

同理，δ_ik 的后验期望也由逆伽马共轭公式给出，你会看到同样的收缩效果—— 特征数少时，尺度因子被拉向全局先验，避免方差被错误放大或缩小。

一张图看懂整个流程

下面这张 SVG 把 Location-Scale 模型的核心数据流画了出来。从左到右， 你能清楚看到：原始数据 → 估计位置/尺度 → 贝叶斯先验 → 后验收缩 → 校正输出。