051

Scaffold Split:为什么随机划分会高估模型能力

随机划分会严重高估模型能力。这篇讲清骨架划分的原理、实现与更严格的替代方案。

骨架划分是分子机器学习中最重要的评测实践之一。它解决的问题很具体:药物数据是成系列产生的,随机划分会让测试集分子的「近亲」留在训练集中——模型只需查表插值就能得高分。

问题的根源

# 【药物数据的产生方式】:
#   化学家围绕一个骨架合成几十到几百个类似物
#   → 这些分子彼此高度相似
#   → 【数据点不是独立同分布的】
#
# 【随机划分的后果】:
#   同一系列的分子被随机分到训练集与测试集
#   → 测试集中的每个分子,训练集中都有它的近亲
#   → 模型只需「找到最相似的训练分子,输出它的活性」
#   → 【这不是泛化,是记忆】
#
# 【典型的性能差距】:
#   同一模型、同一数据:
#     随机划分   R² = 0.85
#     骨架划分   R² = 0.55
#     时间划分   R² = 0.40
#
# 【而实际部署时面对的是什么】:
#   新合成的分子(可能是新系列)
#   → 更接近骨架划分或时间划分的情形
#   → 【随机划分的数字对实际决策毫无参考价值】
#
# 【这个问题有多普遍】:
#   相当比例的已发表工作仍以随机划分为主要结果
#   → 【读论文时必须查这一点】(见 169)

实现

from rdkit import Chem
from rdkit.Chem.Scaffolds import MurckoScaffold
from collections import defaultdict
import numpy as np

def scaffold_split(smiles_list, frac_train=0.8, frac_valid=0.1,
                   generic=False, include_chirality=False, seed=0):
    """按 Bemis-Murcko 骨架划分数据集"""
    groups = defaultdict(list)
    for i, smi in enumerate(smiles_list):
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            continue
        scaf = MurckoScaffold.GetScaffoldForMol(mol)
        if generic:
            try:
                scaf = MurckoScaffold.MakeScaffoldGeneric(scaf)
            except Exception:
                pass
        key = Chem.MolToSmiles(scaf)
        # 【无环分子单独处理,避免全部归为一类】
        if not key:
            key = f"__acyclic_{i}__"
        groups[key].append(i)

    # 【大骨架优先进训练集】
    #   → 让测试集包含更多罕见骨架,测试更严格
    sets = sorted(groups.values(), key=len, reverse=True)

    n = len(smiles_list)
    n_train, n_valid = frac_train * n, frac_valid * n
    train, valid, test = [], [], []
    for s in sets:
        if len(train) + len(s) <= n_train:
            train += s
        elif len(valid) + len(s) <= n_valid:
            valid += s
        else:
            test += s
    return train, valid, test

# 【变体:随机骨架划分】
#   把骨架组随机分配,而非按大小排序
#   → 可以跑多个种子,得到方差估计
def random_scaffold_split(smiles_list, frac_train=0.8, seed=0):
    groups = defaultdict(list)
    for i, smi in enumerate(smiles_list):
        mol = Chem.MolFromSmiles(smi)
        if mol is None:
            continue
        key = MurckoScaffold.MurckoScaffoldSmiles(mol=mol) or f"__acyclic_{i}__"
        groups[key].append(i)
    rng = np.random.default_rng(seed)
    sets = list(groups.values())
    rng.shuffle(sets)
    n_train = frac_train * len(smiles_list)
    train, test = [], []
    for s in sets:
        (train if len(train) + len(s) <= n_train else test).extend(s)
    return train, test
# 【推荐用这个】:能跑多个种子,报告均值 ± 标准差

验证划分是否足够严格

# 【骨架划分不一定就够严格】
#   不同的骨架可能仍然很相似
#   (如苯并咪唑 vs 苯并噻唑)
#
# 【必做的检查】:看测试集与训练集的相似度分布

from rdkit import DataStructs
from rdkit.Chem import rdFingerprintGenerator

gen = rdFingerprintGenerator.GetMorganGenerator(radius=2, fpSize=2048)

def split_severity(train_smiles, test_smiles):
    tr_fps = [gen.GetFingerprint(Chem.MolFromSmiles(s))
              for s in train_smiles if Chem.MolFromSmiles(s)]
    max_sims = []
    for s in test_smiles:
        m = Chem.MolFromSmiles(s)
        if m is None:
            continue
        sims = DataStructs.BulkTanimotoSimilarity(gen.GetFingerprint(m), tr_fps)
        max_sims.append(max(sims))
    max_sims = np.array(max_sims)

    print("测试集分子与训练集的最大相似度:")
    print(f"  中位数     {np.median(max_sims):.3f}")
    print(f"  > 0.7 比例 {(max_sims > 0.7).mean():.1%}  ← 【应该很低】")
    print(f"  > 0.5 比例 {(max_sims > 0.5).mean():.1%}")
    print(f"  < 0.3 比例 {(max_sims < 0.3).mean():.1%}  ← 真正的外推")
    return max_sims

# 【判读】:
#   > 0.7 的比例超过 20% → 【划分不够严格】
#   → 考虑用 generic 骨架或聚类划分
#
# 【更有信息量的报告方式】:
#   按相似度区间分层报告性能

def stratified_performance(y_true, y_pred, max_sims):
    bins = [(0.0, 0.3), (0.3, 0.5), (0.5, 0.7), (0.7, 1.0)]
    print("\n按与训练集的相似度分层:")
    for lo, hi in bins:
        mask = (max_sims >= lo) & (max_sims < hi)
        if mask.sum() < 10:
            continue
        rmse = np.sqrt(((y_pred[mask] - y_true[mask]) ** 2).mean())
        print(f"  相似度 [{lo}, {hi}): n={mask.sum():4d}, RMSE={rmse:.3f}")

# 【这张表直接回答】:
#   「模型在距离训练集多远的地方还可信?」
#   → 【这是最有实用价值的评测输出】

各种划分方式的严格程度

方式 严格度 对应场景 需要
随机 最松 几乎不对应任何真实场景
Murcko 骨架 遇到新化学系列
Generic 骨架 较严 杂环变体也算同类
指纹聚类 整个相似子空间未见过 050《分子相似性聚类》
时间划分 最接近真实 用历史预测未来 时间戳
跨项目/跨靶点 最严 模型迁移 多项目数据

时间划分:最诚实的评测

# 【如果数据有时间戳,就用时间划分】

def temporal_split(df, date_col="date", cutoff="2024-01-01"):
    train = df[df[date_col] < cutoff]
    test = df[df[date_col] >= cutoff]
    print(f"训练 {len(train)} (至 {cutoff}),测试 {len(test)}")
    return train, test

# 【为什么最诚实】:
#   它完全模拟了生产环境:
#     用截至今天的数据训练,
#     预测明天合成的分子
#   → 【自然地包含了化学空间的漂移】
#
# 【滚动评测:更完整的图景】
def rolling_evaluation(df, date_col="date", n_folds=5):
    df = df.sort_values(date_col)
    fold_size = len(df) // (n_folds + 1)
    results = []
    for k in range(1, n_folds + 1):
        train = df.iloc[:k * fold_size]
        test = df.iloc[k * fold_size:(k + 1) * fold_size]
        # ... 训练与评测 ...
        results.append({"fold": k, "n_train": len(train),
                        "n_test": len(test)})
    return results
# 【能看到「数据越多模型越好吗」这个关键问题的答案】

# 【如果内部数据没有时间戳】:
#   → 【现在就开始记录】
#     合成日期、测定日期、录入日期
#   → 这是低成本高回报的数据治理
#   → 【一年后你会感谢自己】

报告评测结果的正确方式

# 【一份完整的评测报告应该包含】:
#
# 1) 【多种划分的结果】
#    随机划分(作为上界参照)
#    骨架划分(主要结果)
#    时间划分(如果有时间戳)
#
# 2) 【多个随机种子】
#    均值 ± 标准差
#    → 【提升必须超过 2 倍标准差才算显著】
#
# 3) 【测试集与训练集的相似度分布】
#
# 4) 【分相似度区间的性能】
#    ← 【最有信息量的部分】
#
# 5) 【与强基线的比较】
#    ECFP + LightGBM,同等调优预算(见 131)
#
# 6) 【适用域的界定】
#    模型在什么范围内可信
#
# 【一个反面例子】:
#   「我们的模型在 XX 数据集上达到 R² = 0.92」
#   → 缺少:什么划分?几个种子?基线是多少?
#   → 【这个数字无法解读】
#
# 【一个正面例子】:
#   「骨架划分(5 个随机种子),R² = 0.55 ± 0.04;
#    ECFP+LightGBM 基线 0.52 ± 0.05;
#    在与训练集相似度 < 0.4 的测试分子上,
#    R² 降至 0.31。
#    因此模型适用于与已有系列相似度 > 0.4 的分子。」
#   → 【这个结论可以直接支撑决策】

关键要点

  • 药物数据是成系列产生的,不是独立同分布的——这是随机划分失效的根源;
  • 无环分子的骨架为空,必须单独处理否则会全部落在同一边;
  • 骨架划分后要检查测试集与训练集的最大相似度分布,确认是否足够严格;
  • 分相似度区间报告性能是最有实用价值的输出——它界定了模型的适用域。

延伸资源