XGBoost-SHAP框架用于筛位点

用XGBoost-SHAP框架可以得到每个特征的贡献值，因此可以用来筛选位点。

本文内容主要来自于AI。

XGBoost-SHAP 框架理论

在现代数据科学与传统科研的交叉领域，存在一个核心痛点：“高精度模型往往不可解释，而可解释的模型往往精度不足”。XGBoost-SHAP 框架正是打破这一僵局的“黄金搭档”。如果说 XGBoost 是榨干数据最后一点预测潜力的“性能猛兽”，那么 SHAP 就是照亮其内部逻辑的“高亮探照灯”。

下面的理论部分都是 AI 写的，我也不知道对不对。我的理解就是通过新提出的 TreeSHAP 算法可以高校对这些决策树算法计算得到特征的贡献值，不然正常是没法计算特征的贡献值的。

前提：二者作用与分工

可以把整个分析流程看作一场考古发掘：

结论：XGBoost 是因，SHAP 是果。没有 XGBoost 的强力预测，SHAP 无据可依；没有 SHAP 的归因解释，XGBoost 只是个黑箱。

1. XGBoost 内置重要性 feature_importances_：负责“粗筛”

• 角色：特征过滤器 (Filter)。
• 原理：基于树模型训练过程中的统计（分裂次数、增益）。
• 优点：极快。对于几万个 SNP，几秒钟就能排出名次。
• 缺点：
  • 不稳定（尤其是小样本）。
  • 无法处理特征共线性（LD）导致的贡献分散。
  • 只有分数，没有方向（不知道是增效还是减效）。
• 你的用法：用来把 2 千万 $\rightarrow$ 5 万 $\rightarrow$ 5 千。这是它最擅长的工作。

2. SHAP 值：负责“精修”与“解释”

• 角色：归因解释器 (Explainer)。
• 原理：基于博弈论，计算特征在各种组合下的边际贡献。
• 优点：
  • 更公平：能正确处理 LD 区域，不会因为两个 SNP 高度相关就胡乱分配权重。
  • 有方向：告诉你 0$\rightarrow$2 是增加还是减少表型。
  • 抗噪：结合交叉验证后，能筛掉“偶尔重要”的噪声位点。
• 缺点：极慢。计算 5 千个 SNP 的 SHAP 值，可能需要几分钟到几小时。

3. 核心区别对照表

但是我看别人论文和网页，并没有使用 XGBoost feature_importances_ 进行初筛位点的操作。

第一部分：XGBoost 的数学本质

XGBoost (eXtreme Gradient Boosting) 的核心思想是加法模型（Additive Model）与前向分步算法。它不是一次性训练一个复杂的模型，而是一步步地添加新的弱分类器（决策树）来修正之前的残差。

1. 目标函数（Objective Function）

XGBoost 的目标函数由两部分组成：损失函数（衡量预测准不准）和 正则化项（控制模型复杂度，防止过拟合）。

给定 $n$ 个样本，训练 $K$ 棵树，预测结果为 $\hat{y}_i$：

$$\text{Obj} = \sum_{i=1}^n l(y_i, \hat{y}i) + \sum{k=1}^K \Omega(f_k)$$

• $l(y_i, \hat{y}_i)$：损失函数（如均方误差、Logistic Loss）。
• $\Omega(f_k)$：第 $k$ 棵树的正则化项。

2. 加法模型与泰勒展开

在第 $t$ 轮迭代中，模型预测为前 $t-1$ 轮的模型加上一棵新树 $f_t$：
$$\hat{y}_i^{(t)} = \hat{y}_i^{(t-1)} + f_t(x_i)$$

代入目标函数：
$$\text{Obj}^{(t)} = \sum_{i=1}^n l(y_i, \hat{y}_i^{(t-1)} + f_t(x_i)) + \Omega(f_t)$$

由于 $l(y_i, \hat{y}i^{(t-1)})$ 是常数，为了优化 $f_t$，XGBoost 使用了二阶泰勒展开来近似损失函数：
$$\text{Obj}^{(t)} \approx \sum{i=1}^n \left[ l(y_i, \hat{y}_i^{(t-1)}) + g_i f_t(x_i) + \frac{1}{2} h_i f_t^2(x_i) \right] + \Omega(f_t)$$

其中：

• $g_i = \partial_{\hat{y}^{(t-1)}} l(y_i, \hat{y}^{(t-1)})$：损失函数的一阶导数（梯度）。
• $h_i = \partial^2_{\hat{y}^{(t-1)}} l(y_i, \hat{y}^{(t-1)})$：损失函数的二阶导数（海森矩阵）。

3. 树的参数化与最优权重

定义树 $f_t(x) = w_{q(x)}$，其中 $q(x)$ 将样本映射到叶子节点，$w_j$ 是第 $j$ 个叶子节点的分数（权重）。
正则化项定义为：
$$\Omega(f_t) = \gamma T + \frac{1}{2} \lambda \sum_{j=1}^T w_j^2$$

• $T$：叶子节点数。
• $\gamma$：惩罚分裂的复杂度（Gamma）。
• $\lambda$：L2 正则化系数。

将样本按叶子节点分组（$I_j = {i | q(x_i)=j}$），目标函数简化为：
$$\text{Obj}^{(t)} = \sum_{j=1}^T \left[ (\sum_{i \in I_j} g_i) w_j + \frac{1}{2} (\sum_{i \in I_j} h_i + \lambda) w_j^2 \right] + \gamma T$$

对 $w_j$ 求导并令其为 0，得到最优叶子权重：
$$w_j^* = -\frac{\sum_{i \in I_j} g_i}{\sum_{i \in I_j} h_i + \lambda}$$

代入后得到衡量树结构好坏的评分函数（Gain）：
$$\text{Gain} = \frac{1}{2} \left[ \frac{(\sum_{i \in I_L} g_i)^2}{\sum_{i \in I_L} h_i + \lambda} + \frac{(\sum_{i \in I_R} g_i)^2}{\sum_{i \in I_R} h_i + \lambda} - \frac{(\sum_{i \in I} g_i)^2}{\sum_{i \in I} h_i + \lambda} \right] - \gamma$$

这就是 XGBoost 选择分裂节点的数学依据：最大化这个 Gain 值。

第二部分：SHAP 的公理化推导

SHAP (SHapley Additive exPlanations) 将博弈论中的 Shapley 值 引入机器学习。其核心在于将模型的预测值 $f(x)$ 解释为基值 $E[f(X)]$ 加上各特征贡献 $\phi_i$ 的和。

1. SHAP 值的定义

对于任意样本 $x$，其 SHAP 值 $\phi_i$ 满足：
$$f(x) = \phi_0 + \sum_{i=1}^M \phi_i$$

• $\phi_0$：模型在训练集上的平均预测值（Base Value）。
• $\phi_i$：特征 $i$ 对当前样本预测值的贡献。

2. 基于博弈论的 Shapley 公式

对于特征 $i$，其 Shapley 值是通过计算它在所有可能特征组合（联盟 $S$）中的边际贡献加权平均得到的：

$$\phi_i = \sum_{S \subseteq N \setminus {i}} \frac{|S|! (M - |S| - 1)!}{M!} \left[ f_{S \cup {i}}(x_{S \cup {i}}) - f_S(x_S) \right]$$

• $N$：所有特征的集合。
• $S$：不包含特征 $i$ 的特征子集。
• $|S|$：子集 $S$ 的大小。
• $f_{S \cup {i}}$：仅使用子集 $S$ 加上特征 $i$ 时的模型预测值。

问题：这个公式的计算复杂度是 $O(2^M)$，对于特征多的模型是不可行的。

3. TreeSHAP：针对树模型的高效算法

Lundberg 等人提出的 TreeSHAP 算法，利用树的结构特性，将复杂度降至多项式时间 $O(TLD^2)$ （其中 T 是树的数量，L 是叶节点数，D 是树的最大深度）。

核心思想：
对于树模型，预测值是样本落入叶子节点的权重之和。TreeSHAP 通过追踪样本路径，计算“包含特征 $i$ 的路径”在整棵树中的概率分布。

对于单棵树，特征 $i$ 的贡献 $\phi_{i,j}$ 计算为：
$$\phi_{i,j} = \sum_{k=1}^{T_j} \frac{\text{PathCoverage}(i, k)}{T_j} \cdot \Delta \text{Value}$$

• PathCoverage：特征 $i$ 在节点 $k$ 处对样本路径的覆盖率（权重）。
• $\Delta \text{Value}$：由于分裂带来的预测值变化。

最终 SHAP 值是所有树的贡献之和：
$$\phi_i = \sum_{j=1}^K \phi_{i,j}$$

(这也是为什么 SHAP 对树模型（XGBoost, LightGBM, CatBoost）的支持如此无脑且快速的原因)

第三部分：XGBoost + SHAP 的结合

交互效应的数学表达：
SHAP 还可以计算特征 $i$ 和 $j$ 的交互值 $\phi_{i,j}$：
$$\phi_{i,j} = \sum_{S \subseteq N \setminus {i,j}} \frac{|S|! (M - |S| - 2)!}{M!} \left[ f_{S \cup {i,j}} - f_{S \cup {i}} - f_{S \cup {j}} + f_S \right]$$

这解释了为什么在你的 GWAS 分析中，使用 SHAP 依赖图（Dependence Plot）能看到两个基因位点如何共同作用影响表型（例如：当基因 A 突变时，基因 B 的影响才会显现）。

总结

XGBoost 通过二阶优化保证了预测的高精度，而 SHAP 通过公理化归因保证了逻辑的可解释性。二者结合，实现了从“知其然”（预测准确）到“知其所以然”（机制解析）的跨越。

XGBoost-SHAP 框架技术路线

🟢 阶段一：数据工程与预处理 (Data Engineering)

这是决定模型上限的基础环节，主要完成从原始数据到机器学习可用特征的转化。

• 多源数据整合：收集自变量（如遥感影像的波段反射率、建筑形态指标、SNP位点基因型）与因变量（如LST、表型性状）。
• 数据清洗与标准化：处理缺失值与异常值；对不同量纲的特征（如温度与面积）进行标准化或归一化处理。
• 数据集划分：按比例（如 8:2 或 7:3）切分为训练集 (Training Set) 与 测试集 (Test Set)，确保模型评估的无偏性。

🟡 阶段二：XGBoost 模型构建与优化 (Black-box Prediction)

此阶段致力于打磨一个高精度、强泛化能力的“黑箱”预测模型。

• 模型初始化与训练：调用 XGBRegressor或 XGBClassifier（要求因变量是0 1 2 形式的整数，又可分为二分类和多分类），通过加法模型与前向分步算法，迭代训练出由多棵决策树组成的集成模型。
• 超参数调优 (Hyperparameter Tuning)：
  • 利用网格搜索 (Grid Search) 或贝叶斯优化，确定最佳参数组合（如 max_depth, learning_rate, n_estimators, reg_alpha等）。
  • 早停法 (Early Stopping)：监控验证集误差，防止模型过拟合。
• 性能评估：在测试集上使用 R2、RMSE（回归）或 Accuracy、AUC（分类）等指标评估模型预测精度。(此时模型已具备高预测能力，但缺乏可解释性)。

🔵 阶段三：SHAP 机制解析与归因计算 (Interpretability)

这是连接“黑箱”与“科学机理”的桥梁，利用博弈论拆解模型的内在逻辑。

• 构建 SHAP 解释器：调用 shap.TreeExplainer(model)，基于 TreeSHAP 算法 O(TLD2)的高效复杂度，计算测试集中每个样本、每个特征的 Shapley 值。
• 全局特征重要性量化：聚合所有样本的 SHAP 值绝对值，得出各变量对模型预测的总体贡献度排名，筛选出关键驱动因子。
• 局部预测逻辑拆解：针对极端样本或特定个案，提取其特征的 SHAP 值，还原模型对该样本的具体决策路径。

具体解释

阶段三是整个框架的**“灵魂”。它将第一阶段（数据）和第二阶段（模型）的成果，转化为人类可理解的因果逻辑**。这一阶段不再关注“预测得准不准”，而是关注“为什么准”以及“是怎么决策的”。

1. 核心操作流程

在代码中，这个阶段通常对应以下几步：

1
2
3
4
5

# 1. 初始化解释器（绑定训练好的黑盒模型）
explainer = shap.TreeExplainer(model)  # model 是第二阶段训练好的 XGBoost

# 2. 计算 SHAP 值（这是最关键的一步）
shap_values = explainer(X_test)  # 传入数据

这里的 shap_values是一个数组，形状为 (样本数, 特征数)。例如，如果你有 1000 个测试样本，20 个特征，那么 shap_values就有 1000 行、20 列，每一个数字代表**“某个特征对某个样本的预测贡献了多少”** 。

采用 “平均绝对 SHAP 值” 对特征进行排序（这里实际算的是该特征在所有测试样本中的贡献值绝对值的均值，而不是该特征每增加1个单位的增加的效应值），这是 SHAP 的标准做法。

这种做法其实衡量的是“这个特征在群体中平均有多大的话语权”，而不是“这个特征每变动一个单位会带来多少线性变化”。

1
2

# 计算每个 SNP 的平均贡献度（跨所有样本）
importance = np.abs(shap_values).mean(axis=0)

我问了 AI ，这种对于我这种非0/1类的特征（如 SNP 和 连续变量），这种做法也适用吗？AI 说是，而且有文章里面就有连续变量的特征。AI 的解释我看了都不太通，只有一个点比较通，那就是效应的非线性，强行算单位变化会丢失信息。

2. 灵魂拷问：为什么要算“测试集”？训练集不能算吗？

答案是：技术上完全可以算训练集，但在科学分析和论文发表中，强烈建议使用测试集。

以下是为什么我们要坚持用测试集（Test Set） 来计算 SHAP 值的深层原因：

原因一：模拟“真实世界”的决策逻辑（泛化性）

• 训练集：模型见过无数次，甚至已经把噪声都记下来了（过拟合）。在训练集上解释模型，可能解释的是“模型如何记住了数据”，而不是“模型学到了规律”。
• 测试集：模型从未见过这些数据。在这里计算 SHAP 值，反映的是模型面对未知新样本时的真实决策逻辑。这才是科研关心的：“如果我用这个模型预测新的病例/新的城市/新的基因数据，它会依据什么下判断？”

原因二：避免“数据泄露”导致的逻辑谬误

如果你用训练集计算 SHAP 值，可能会出现一种假象：某个特征看起来很重要，但实际上只是因为它在训练集中恰好和标签高度相关（甚至是巧合）。

• 举例：在训练集中，可能因为采样偏差，“穿红色衣服”的人患病的多。模型在训练集上会给“红色衣服”很高的 SHAP 值。但在测试集（真实世界）中，这两者无关。用测试集计算就能揭露这种虚假关联。

原因三：符合“训练-验证-测试”的科研范式

在《Building and Environment》那篇论文中，逻辑链条必须是闭环的：

1. 训练集：用来拟合参数（让模型聪明）。
2. 测试集：用来评估性能（证明模型聪明），同时也用来解释模型（证明模型聪明的理由是合理的）。

如果解释用的是训练集，审稿人会质疑：“你只是解释了它怎么背答案，没解释它怎么解题。”

3. 什么时候可以用训练集算 SHAP？

虽然推荐用测试集，但在以下特定场景下，训练集也是有价值的：

4. 结合GWAS分析

在 GWAS（全基因组关联分析）中，这一点尤为重要：

• 如果你用 训练集 SNP 计算 SHAP，可能会发现一些 SNP 很重要，但它们可能只是恰好在你采集的样本群体中出现频率高（假阳性）。
• 如果你用 测试集 SNP 计算 SHAP，发现这些 SNP 依然有高贡献，这就能强有力地证明：这些基因位点确实通过生物学机制影响了表型，而不是统计噪音。

总结：

计算测试集的 SHAP 值，是为了确保我们解释的是模型的**“智慧”（泛化规律），而不是模型的“记忆”**（过拟合噪声）。这是从“机器学习跑通了”迈向“科学研究完成了”的关键一步。

🔴 阶段四：科学机理可视化与洞察 (Visualization & Insights)

将高维的 SHAP 数值转化为直观的图表，从而推导出符合领域知识的科学结论。

• 特征重要性摘要图 (Summary Plot)：宏观展示哪些因素主导了目标变量的变化。
• 依赖图与交互效应 (Dependence Plot)：揭示单一特征的边际效应（如：随着建筑容积率增加，LST如何非线性变化），并挖掘特征间的协同或拮抗作用（如：风速与湿度对降温的联合影响）。
• 个体瀑布图 (Waterfall Plot) / 力图 (Force Plot)：微观解释单个预测值是如何由基线值（E[f(X)]）加上各特征的贡献值（ϕi）累加而成的。

升级技术路线-双层交叉验证

当样本量较小时，单次划分训练集/测试集确实存在较大的随机性风险。在 GWAS 这种高维小样本场景中，5×交叉验证（5-fold CV）不仅是“锦上添花”，更是确保科学结论稳健性的必要条件。

以下是针对小样本优化的 XGBoost-SHAP 技术路线升级版，将交叉验证逻辑深度融入各个阶段：

🟢 阶段一：数据工程与预处理（不变）

• 标准化特征，确保数据质量。

🟡 阶段二：XGBoost 模型构建与优化（核心改动点）

不再使用单次划分，而是引入 嵌套交叉验证（Nested CV）：

1. 外层循环（模型评估）：
   • 将数据分为 5 折（5 folds）。
   • 每次取 4 折训练，1 折测试。重复 5 次，确保每一个样本都被预测过一次。
2. 内层循环（超参数调优）：
   • 在每一折的训练集内部，再次使用 5 折交叉验证（或贝叶斯优化）来搜索最佳的 max_depth、learning_rate等参数。
   • 目的：防止在小样本上，因为超参数刚好适配了某一次特定的训练/测试划分而产生“伪最优”。

🔵 阶段三：SHAP 机制解析与归因计算（关键升级）

这是解决你痛点的核心。不再只计算一个测试集的 SHAP 值，而是计算 5 个。

1. 循环计算：
   • 在 5 折交叉验证的每一轮中，使用训练好的模型计算当前测试折的 SHAP 值。
   • 你会得到 5 份 shap_values数据（对应 5 个不同的测试子集）。
2. 结果聚合（Aggregation）：
   • 均值法（最常用）：在每一次交叉验证折（Fold）中，计算该折测试集上所有样本的 SHAP 值绝对值的平均值（即该折的特征重要性排名）。计算这 5 个数字的 均值（Mean） 和 标准差（Std）。
   • 中位数法：如果数据中有极端异常值，使用中位数更稳健。
   • 一致性检验：统计某个 SNP 在 5 次实验中是否都排在 Top 10。只有那些**“次次都稳居前列”**的 SNP 才被认为是真正的强关联位点。

🔴 阶段四：科学机理可视化与洞察（基于聚合结果）

• 使用聚合后的 SHAP 值（而不是单次的结果）来绘制 Summary Plot 或 Dependence Plot。
• 在论文中，你可以自信地陈述：“经过 5 次交叉验证，特征 X1 的平均 SHAP 值为 0.85 ± 0.12，表明其对表型具有稳定且显著的正向影响。”

挖掘表型关联SNP位点技术路线

针对挖掘表型关联SNP位点，由于一般是样本少，位点多的情况。AI 的建议是使用机器学习进行验证位点，而不是发现位点。

1. 先用 GEMMA、GCTA 或 Plink 做传统的 GWAS 关联分析，算出 P 值。
2. 进行 clumping ，基于 P 值结果过滤 LD ，剔除连锁位点
3. 过滤后，选取 Top 1000 个 P 值最小的 SNP 。
4. 将剩下的位点丢给上面的 XGBoost-SHAP 框架，基于 SHAP 值对SNP位点进行排序，进行最终的位点筛选

最终位点的筛选最简单的思路就是人为设定一个数目，选取top n 的位点，但是会引发一个问题，比如选前 50，别人会问你为什么是 50 ？

别人的网页没有这方面内容，他们都没有筛选特征，只是用 shap 画图。结合与 AI 的问答，我自己想了一个比较自动的标准。首先必须要做 5× 交叉验证，得到 SHAP_Mean 和 SHAP_Std ，并计算变异系数 （std/mean, CV）。但是由于不同 fold 之间可能 SNP 的基因型分布差异比较大（随机划分的），因此 CV 更多反映 fold 间样本结构差异 ，而不是 SNP 本身是否稳定（小样本量会加剧这一问题），因此看 CV 没有用（当样本量很大时，fold间样本结构就趋于一致，此时的 CV 就反映的snp本身的稳定性）。

筛选位点的唯一金标准就是 SHAP_Mean 大于给定阈值（可以看看每个 fold 中的 shap_mean ，看看是不是只在1折里不是0，其他折都是0，这种可能就是噪音，进行剔除）。大样本量时可以关注变异系数 CV 。

为什么选Top 1000

针对你的 109 个样本 和 重测序数据，结合传统 GWAS 先验知识，以下是具体的参数建议和原理解析。

结论：选 1000 比选 5000 更安全、更科学。

建议：优先选择 Top 1000。如果担心漏掉信号，可以在做完 SHAP 后，观察第 1000 名的 P 值和 SHAP 值是否断崖式下跌，如果没有，再放宽到 2000。

为什么要过滤 LD 吗？

结论：必须过滤，且这是最关键的步骤之一。

即便你已经用 GEMMA/GCTA 选出了 Top 1000 个 P 值最小的 SNP，LD 过滤依然不可或缺。原因如下：

1. 消除冗余信号（去重）

   • 传统 GWAS 的 P 值往往成簇出现。比如 1 号染色体上 100 个 SNP 都显著（因为 LD 高），但它们代表的是同一个因果变异。
   • 如果不去 LD，你的 Top 1000 列表里可能有 300 个都来自同一个基因区域。这浪费了你本就稀缺的“1000 个名额”。

2. 保障 XGBoost 稳定性

   • XGBoost 在高度相关的特征面前会“选择困难”，导致每次训练选中的特征不同，使得 feature_importances_ 和 SHAP 值波动极大。

3. 提升 SHAP 归因准确性

   • SHAP 虽然能处理共线性，但在小样本下，去除 LD 能让贡献度更清晰地分配给“代表位点”，而不是分散在连锁位点中。

3. 针对你的数据的“黄金标准”流程

请严格按照以下顺序操作，以确保结果能被审稿人接受：

4. 为什么这样做最稳妥？

• 生物学上：你利用了传统遗传学的先验知识（P 值），保证了候选位点有遗传基础。
• 统计学上：你剔除了由于连锁造成的冗余，保证了特征独立性。
• 机器学习上：109 个样本对应几百个特征，模型自由度充足，不会过拟合。

总结：
不要直接把 Top 1000 丢给 XGBoost。一定要对这 1000 个再做一次 LD 修剪。 这一步能防止你的“关键基因列表”里塞满同一个基因的重复拷贝，是产出干净、可信结果的决定性环节。

回归问题实操

输入文件：

1. 基因型为 012格式的 qc_prune.raw
2. 表型文件 blue_pick.txt （个体号和 blue 值，无标题）

参考文献

1. https://mp.weixin.qq.com/s/Gv5GTBX5lKXEd_4u_U7NPA
2. https://www.bilibili.com/read/cv39836525/?opus_fallback=1