Statistical summaries of unlabelled evolutionary trees

作者: Rajanala Samyak, Julia A Palacios
来源: Biometrika
主题: 统计计算 / 算法
相关性: 3/10
机构绿灯: Stanford University（US News 前 50，免分进入精读）
链接: 期刊页 · arXiv

---

一、领域脉络与小综述

这个方向是什么

这个子方向关注的是对无标签（unlabelled）的有根进化树（rooted phylogenetic trees）进行统计总结。其根本问题在于：给定一个树形数据的集合（例如来自贝叶斯MCMC的后验样本，或来自不同方法/区域的树），如何定义并计算其“均值”、“方差”和“分位数”等经典统计量？由于树的数学空间不是欧几里得空间，这些概念需要被推广，而这篇文章聚焦于无标签、有排名（ranked）的二叉树和带枝长的谱系（genealogies）。

该方向当前成熟度中等偏上：对有标签树的空间已进行了大量研究（如Billera-Holmes-Vogtmann（BHV）树空间及其上的Fréchet均值）。而对于无标签树，研究相对较少，本文是该方向（特别是从算法和实用角度）的重要推进。

发展脉络

• 奠基工作：标号树的Fréchet均值与BHV树空间。 Brown & Owen (2017, 论文[17]) 提出了在BHV树空间中使用Fréchet均值和方差来总结一组有标号系统发育树集。这是将传统统计量推广到树空间的一个奠基性、概念性的工具。然而，树空间的结构复杂，BHV空间是“黏性”的（Fréchet均值常处于“零测度”的退化区域），且算法效率在样本量大时受限。该工作主要处理的是“有标签、无排名”的树。

• 主要进展：无标签树的拓扑（形状）与距离度量。 在更粗粒度的“树形（tree shape）”和“有排名树形（ranked tree shape）”上，已有模型工作，如Aldous (1996, 2001) 的β-分裂模型，以及Sainudiin & Véber (2016, 论文[18]) 的Blum-François分裂模型，它们为生成随机树形提供了分布。

• 关键突破点： Kim et al. (2019, 论文[4]) 提出了一个用于有排名树形和谱系的距离度量空间。这是本文的直接数学基础。他们引入了“F-矩阵”表示法，并证明了空间是一个度量空间。但是，这个工作的重点是定义距离和证明度量的性质（如三角不等式、等距性等），而没有系统性地解决如何在此空间上进行统计总结（如计算均值）的问题。

• 当前前沿：从有标签到无标签。 Heled & Bouckaert (2013, 论文[9]) 和 Benner & Bačák (2013, 论文[16-17引用]) 等工作研究了如何从有标签树的后验样本中总结出一个代表性树。然而，当树是无标签时，直接应用这些方法会面临“对应（correspondence）”问题——即如何匹配不同树上的叶子；而且，这些方法通常是启发式的（例如最大支系可信度树），而非基于严格的几何（Fréchet均值）框架。

• 本文的位置： 本文正式地将Fréchet均值、方差和四分位集的概念，应用到Kim et al. (2019) 定义的无标签有排名树的度量空间中。它提供了一个高效的组合优化算法来在地图上找到Fréchet均值（即中心点），从而使得在这些有趣的、无标签的空间中进行统计总结成为可能。因此，这篇论文是在前人的距离度量定义（Kim et al.）与有标签树的统计总结（Brown & Owen）之间搭起了一座算法桥梁。

子线索聚类

这些被引文献大致落在三条子线索上：

1. 有标签树的总结与推论 (Phylogenetic Summary for Labelled Trees):

   • 做什么： 如何从有标签（和通常无排名）的系统发育树样本中生成一个代表性树？常见方法包括：最大支系可信度树 (MCC, 论文[13])、多数规则共识树、或BHV树空间中的Fréchet均值（论文[17]）。
   • 代表工作： Heled & Bouckaert (2013), O’Reilly & Donoghue (2017), Brown & Owen (2017)。
   • 限制： 这些方法要求树是带有叶子标签的，无法直接处理无标签树。例如，一棵树将“加州”和“佛罗里达”的样本（叶子）分开，另一棵树可能相同，只是叶子标签换了位置 —— 对于无标签树，它们实际上是同一个形状，但标号树处理方法会认为不同。

2. 无标签树形的生成模型与概率 (Generative Models of Unlabelled Tree Shapes):

   • 做什么： 统计建模无标签树形（包括有排名和无排名）是如何产生的。生物学可以驱动如何“分裂”物种（例如，一棵树的拓扑形状受物种形成或灭绝事件的影响）。
   • 代表工作： Aldous (1996, 2001), Maliet et al. (2018), 以及Sainudiin & Véber (2016)。
   • 作用： 这些模型提供了理论的概率分布，本文的总结方法（如Fréchet均值）可以用来被这些模型在计算机上进行测试，并可用来量化比较真实分布与这些模型生成的结构。

3. 无标签树间的距离度量与算法 (Metrics & Algorithms for Unlabelled Trees):

   • 做什么： 在无标签树的空间上定义有意义的距离，这是进行统计推断的几何基础。
   • 代表工作： Kim et al. (2019), Chakerian & Holmes (2010, 论文[12], 处理有标签但只统计拓扑）。
   • 关键点： Kim et al. (2019) 发明了一种基于“F-矩阵”的表示法，并证明它能定义出一个度量空间。这是本文所有计算的物理基础——没有这个度量，无法谈论Fréchet均值。
   • 本文的角色： 论文[4]并不是直接提供“总结”的手段（例如均值、方差），而主要是定义了上述度量。本文意识到了这个缺口，直接填补了它：在Kim et al. 的度量空间中，回答“如何计算趋势和离散度”。

这个方向在追问的核心问题

1. 计算可行性 (Computational feasibility): 对于无标签树，计算Fréchet均值是一个NP难问题。如何设计高效的、可伸缩的算法来找到它？本回答中，作者运用了组合优化策略，将问题转化为一个在F-矩阵测地线上的整数线性规划（MIP）。
2. 统计的有用性 (Statistical utility): 这些新的统计量（均值、方差）是否真正有助于解决生物学或流行病学问题？它们能否揭示传统方法无法捕获的模式？
3. 与其他范式兼容 (Compatibility): 这些基于距离的方法与基于模型的方法（如BEAST中的贝叶斯方法，论文[1]）如何互补？它们是能够替代，还是只能作为辅助？
4. 不确定性的量化 (Uncertainty quantification): 如何进一步在这类非欧空间中构建置信区间或进行假设检验？作者在文中明确指出未来的方向中包括这种分析（如论文[10]用于随机对象的Fréchet方差分析）。

⚠️ 作者的framing / 竞争与缺失

• 作者的框架： 作者将缺口界定为“对无标签树缺乏总结工具”。已有大量工具（共识树、MCC树）可以总结有标签树的样本；也有距离度量（Kim et al., 2019）可以比较无标签树。但作者认为，缺的正是将运筹学（距离度量）和统计总结（概率均值、方差、分位数）之间连接起来的算法。他们把自己的框架明确地定为“使用组合优化来填充这个空白”。
• 哪些竞争路线被淡化/回避？
  • 模型复杂度： 他们假设整个树是有排名的。这很重要，因为排名（即内部节点的时间顺序）相比无排名形状要大得多。但排名在建立生成模型时非常有用（因为可以对分裂时间进行排序）。这不是逃避，而是一个关键假设，但也意味着方法不能直接处理无排名或未排序的拓扑。
  • 与贝叶斯后验的直接比较： 论文比较了样本之间的Fréchet均值（例如，用BEAST为每个状态跑的树）。但是他们回避了矛盾：在贝叶斯框架之外计算均值（例如，用MIP解）是一个与其后验“不一致”的总结——这意味着你不能直接用它来作为后验的替代品或“点估计”。它只是对树的分布的一种全局的、唯象的几何摘要。
• 什么明显该被引或该存在，却没出现在intro里？
  • 高维随机矩阵 / 经验谱分布： 本文所使用的F-矩阵是一种被固定到树上的组合对象。然而，当将树视为随机对象（例如研究集合的平均树与真实的内在结构有何不同）时，可能会用到更重的平均场近似工具，例如随机图谱论或随机矩阵的线性特征分析——但这都尚未被利用。这是由于方法学的基本不同：树是离散几何，而非带噪音的随机矩阵。这一点未被提出，而且符合本文目标：在给定距离度量的前提下解决优化问题。

张力

在引用的工作中，没有显著的、彼此矛盾的结论。所有被引用的工作（标号树总结、距离度量、拓扑形状建模）在整体目标上是协同的。唯一的“对立面”可能在于：使用Fréchet均值（Brown & Owen, 2017）与使用最大支系可信度树（Heled & Bouckaert, 2013）之间存在明显结果差异。但作者没有把它们作为一个矛盾点去强调，而是将它们归为“总结方法的不同哲学”，这符合论文的基调（即他们提供了一种选择，而非竞争性的颠覆）。整体而言，未见明显的对立引用。

二、最核心、最简单的例子 / 数学问题

第一步：把符号、模型、可观测数据交代清楚

• 符号：

  • 树 (Tree) ： 一棵进化树，包括内部节点（祖先，coalescence events）和叶子（采样，leaves）。树是有根(根节点root)、有向的（从根到叶子）。
  • n： 叶子数（样本数）。
  • τ: 一个有排名（ranked）的无标记（unlabelled）二叉树（binary tree）。排名意味着内部节点相对于根是排序的（按时间顺序）。无标记说明叶子没有单独的标签。
  • F-矩阵 (F-matrix)： Kim et al. (2019) 的表示方法。对于一棵有n片叶子的树，是一个n×n的矩阵，其第i行、第j列的元素表示叶子i和j的最远共同祖先出现的时间（排序，或者说“深度”）。此矩阵能唯一确定一棵无标签、有排名的树。
  • G: 一棵谱系树(genealogy)。在树的基础上，为每条枝（edge，连接两个节点，例如连接coalescence事件和叶子）添加枝长（分支长度，表示时间或期望替代数）。所以 G = (τ, L)，L 是长度向量。
  • 距离 (d): 对树形状τ（或谱系G）之间的一个度量。作者直接使用Kim et al. (2019)的距离，该距离是基于对F-矩阵进行某种路径积分定义。对G而言，还包括枝长的L1距离。
  • 样本 (Sample): T_1, T_2, ..., T_m。每一个T_i都是树形状τ_i或谱系G_i（可观测）。

• 模型：

  • 概率模型（生成模型）： 树（形状或谱系）被视为从一个未知的概率分布P中抽取的独立同分布样本。
  • 生成的本质： 树结构的形成。对于形状τ，结构可以被独立地绘制（由诸如Kingman共祖或Yule模型等随机过程生成）。对于谱系，枝长也被联合生成。
  • 目标： 在观测到的树样本上，计算Fréchet均值、方差等经典统计量。

• 可观测数据：

  • 研究者实际能观测到的是什么： 研究者手中有一个树集合（例如，一条后验MCMC链输出的5000棵树，或来自不同计算中心的树集合）。每个成员都是一棵无标签、有排名、有时（对于谱系）有枝长的二叉树。
  • 什么是潜在/不可观测、只能靠假设去识别：
    • 分布P（总体分布）：是隐藏的。我们只看到一系列样本，假设它们来自某个P。我们并不直接知道P的Fréchet均值。
    • Fréchet均值的“存在性”：在理论上，均值总能通过求距离和的最小值来定义，但样本多的时候可能发现多个点都是最优解（非唯一）。但这不叫不可观测，它只是可能不唯一。 由于算法可以从集合中选出一个中心点，可以计算方差。对于更精细的推断（如置信区间），我们需要借助中心极限定理使用估计方差，这在本文中未来工作提到（论文[10]），但本文对该“真值”的存在性视为既定。

第二步：讲最小内核（假设n=3 的最小无标签有排名树空间）

最简特例：只有三片叶子 (n=3)

1. 设定： 我们有三片叶子（样本s1, s2, s3）。在无标签、有排名的二叉树中，所有可能的树形状排名形状；由于n很小，所有无标签的形状实际上是等价的吗？ 实际上，对于n=3，只有一种拓扑形状（一个内部节点上生下三片叶子）。然而，有排名（内部节点的排序）却给出了两类不同的无标签形状： - 形状1： 两个叶子的共同祖先等于根节点（所有样本同时聚合）。属于一个内部节点的可能性：叶子1和2先配对，然后3加入。排名的唯一顺序是“分裂事件1发生在分裂事件2之后”。

2. 树的数学形式，F-矩阵： 为了直观（因为n=3太小，F-矩阵并不比图更清楚），目标是看距离度量如何定义Fréchet均值。

空间中所有的树（n=3）在距离度量d下形成一个有限度量空间（只有少量离散的点，加上连续体（当考虑枝长时））。Fréchet均值定义为：

均值的树 = arg min_{T in space} sum_{i=1 to m} [ d(T_i, T) ]^2

3. 如何计算？ - 如果我们有一组m个形状样本（例如，3棵随机树）。空间包含少量形状（对于n=3，有1种形状）。实际上，所有树要么是这种形状。两个这种树之间的距离要么是0（如果相同）或正数（由于排名不同？）。假设所有样本都是相同的形状（所有都是左边结构）。那么，Fréchet均值就是那个形状——微不足道的。

• 有意思的情形： 考虑n=4（四种形状）并混合样本。这时，不同的样本可能驻留在不同的形状上（比如50%样本形状1；50%样本形状2）。Fréchet均值可能是形状1，形状2，或者一个“中间”形状（不可能，中间形状不存在因为度量不仅离散，而且它是测地距离）（形状之间路径穿过别的树，空间是凸的，但中心可能是其中一棵形状，也可能卡在它们之间的测地线上）。这揭示了找中心实际上是选一个点，使得所有点到该点之平方距离和最小——求解一个“一维离散的优化问题”。

4. 这个例子展示的核心思想（也是本文核心）： 即使是n=4，当形状不同时，找出中心（Fréchet均值）是一个 组合优化问题。不可能通过计算欧几里得平均值来直接找一个不存在的“混合”形状；而是在有限集合（多个离散形状）里挑选最优的一个；对于谱系，还要加上长度。

作者解法的最小内核就是： 利用Kim等人（2019）的距离函数，此距离（在形状空间中）并未达到“解析”或“闭式”表达，而是通过路径积分得到的数值。作者说：我们不手动计算这个积分，而是使用离散化方法将空间划分为小块（利用秩的性质），把最小化问题归结为一个0-1整数线性规划问题，然后用MIP（混合整数规划）求解器（如Gurobi）暴力求解。“寻找中心”变成了“在这个离散集中找到一个最小化某个线性函数的目标”。

三、这篇论文做了什么

• 三句话：

  1. 本文研究的问题是：如何对无标签、有排名的二叉树形状和谱系（带枝长的树）进行统计总结（计算Fréchet均值、方差、四分位集）。
  2. 核心工具方法是：利用 Kim et al. (2019) 提出的无标签树距离度量，以及本文提出的将Fréchet均值计算问题转化为整数线性规划（MIP）问题的算法。
  3. 主要结论是：该算法能够处理多达100-200棵树的样本，可以系统地生成均值、方差等摘要，并在模拟和真实的SARS-CoV-2进化树分析中展示其有效性和与已知模型的吻合度。

• 关键设定与假设（在第二节最小记号基础上补全）：

  • 定义、记号：

    • 无标签、有排名的树： 叶子没有标签。树有序（内部节点的分裂顺序被确定）。二叉树：每个内部节点恰好有两个子节点。
    • 谱系（Genealogy）： 除了节点顺序，每个枝长都有具体值（代表时间或替换数）。空间是离散的形状形状集 × 连续的正实数（枝长）。故Fréchet均值的计算要同时考虑形状选择和枝长分配。
    • F-matrix（F-矩阵）： 如第二节所述，对于n个叶子，这是一个n×n上三角矩阵，其(i,j)元给出叶i和叶j的共同祖先在树上的秩（即它们的共同祖先在内部节点事件中的排序编号）。保证了距离的度量定义。
    • 距离度量d(·,·)：Kim et al. (2019) 证明：对于形状，d(·,·)由它们的F-矩阵差的某种范数给出；对于谱系，还要加上枝长之差。

  • 假设：

    • 可性： 树样本来自或可以近似为度量空间中的一个点；所有两个点之间的树测地线存在且唯一（在谱系空间中是CAT(0)空间）。
    • Fréchet均值的定义是良好定义的： 距离平方和的最小值存在。作者使用“Fréchet均值”，暗示了即便空间是CAT(0)（凸的，保证唯一性），但对形状空间（不是CAT(0）），均值可能不唯一。这被作者明确承认——其算法能够输出一个最优解（可能是多个之一）。
    • 排名的固定假设： 所有树必须是排名树。引入的非排名或polytomies（多分支）的情况被排除在外。

  • 相比已有文献的强化：

    • 相对有标签树（Brown & Owen, 2017）： 无标签+有排名。作者面对的计算复杂性（对应问题显著困难）远大于有标签情况，因为无标签会造成巨大的搜索空间（必须确定哪个叶子对应哪个）。
    • 相对无标签距离度量（Kim et al., 2019）： 本文不再是“证明是度量空间”的纯理论，而是提供一个计算均值/方差的算法，对空间进行统计学的解读。

• 主要结果：

  • 定理8 (Fréchet均值的组合刻画): （定理陈述：作者证明了找到形状空间下的Fréchet均值等价于求解一个整数线性规划问题。MIP的解是一个或多个“祖先-后代”关系的配置，对应最优的树结构。）
    • 直觉： 不可能直接写出均值的解析式。作者通过引入“加权F-矩阵”概念，将问题变成一个线性的目标函数，而整数约束保证了输出是一棵合法的无标签树。
    • 必要条件： 树的头节点（root）对应的秩是n-1；F-矩阵必须满足某些对称和单调性条件（三角不等式）。这些约束是MIP的线性约束。
    • 技术难点： 如何保证解的整数性？MIP用分支定界法求解。
  • 定理9 (均值的不唯一性/离散性)： 指出在形状空间中，Fréchet均值可能出现在测地线的中间点（而非一棵具体的树）。作者的处理方式是：如果中间点出现，将那个“连续”测地空间中的一个具体有理对称点（例如，测地线的一半）选定为代表。这实际上通过MIP的“加权”自动完成了。
  • 定理10 (Fréchet方差的计算)： Fréchet方差被定义为所有样本到Fréchet均值的距离平方和的样本平均。计算是直接的：确定均值后，用距离度量算一下就行了。

• 证明路线与技术技巧（理论型必写）：

  • 整体路线： （清晰四步： (1) 组合表示：将树形状完全用一个“加权F矩阵”表示。加权的权重是整数（1~某个最大值），这是对混合测地线的建模。 (2) 目标函数线性化：证明：Fréchet均值(样本到树的总距离方和) 可以被改写为 (样本的F矩阵的平均值) 到 (候选树的F矩阵) 之间的一个加权二次范数的和，而由于加权F矩阵与候选树的F矩阵有线性关系，目标函数变为线性。 (3) 约束线性化：合法的树形状（F矩阵）的约束包括：三角不等式、单调性、整数性，均可以写为线性不等式。 (4) 求解：调用MIP求解器（Gurobi）求解这个标准的0-1整数线性规划问题。输出一个整数解，再反解回树形状。对于谱系树，在处理完形状后，将枝长代入做二次凸优化求解。

  • 关键跳跃点：

    • 引理7： 证明候选树的加权F矩阵到真实样本在测地线上的变化等于一个标量权重乘以一个确定矩阵，使得问题完美转化为线性。核心难点的处理：走哪条测地线？由于树空间是CAT(0)，任何两点间测地线可唯一参数化，而在这个测地线上，距离就是Weighted叶距离矩阵的线性变化。
    • 引理8（不好理解但关键）： 实值变量可以伪装成F矩阵（满足整数约束）。他们把寻找“一棵树”的问题，变成了 “寻找一系列简单的二分化（coalescence）” 的问题，并通过多树之间的得当加权让测地线中间的树也可以用整数线性规划表示。

  • 技术技巧点名：

    • 整数线性规划 (Integer Linear Programming, ILP)：最核心的黑盒子。作者没有自己发明规划求解器，而是将树空间分析建模成标准ILP问题，并用商业求解器（Gurobi）求解。
    • 组合枚举与剪枝：在MIP的建模中，对于n<10个小样本，可以提前列出所有可能的F矩阵（上界），减少变量数量。
    • 赋范双线性 / 测地线拆分：形状空间中，不同树之间的测地线被“拆分”成若干段；均值损失函数可以写成每段距离的平方和。这种通过测地线分段并把积分损失化为凸线性组合的技巧是统计几何里的经典技巧（但不是初次应用）。

• 真实例子与应用（务必讲透）：

  • 数据/场景： 使用SARS-CoV-2病毒序列。数据来自GISAID数据库，选取美国四个州（California, Florida, Texas, Washington）在2020年2月至9月的序列。每个州分别用BEAST跑一个贝叶斯系统发育推断（使用HKY模型，固定的替换率8e-4/位点/年，谱系先验使用Skyride模型），并从后验中抽取2000个谱系样本（带枝长，有排名）。

  • 作者怎么把本文方法用上去： 他们将后验样本抛弃标签，使树变为无标签化，然后使用他们提出的算法计算各州丛的Fréchet均值树（形状+枝长）、方差，以及四分位集。他们比较了四个州的均值谱系，并展示了树分布的差异。

  • 实验验证（模拟部分）：

    1. 模拟数据： 从已知分布（如Coalescent、Yule、Beta-Splitting模型）生成的树集，计算平均形状；观察所得均值是否落回正确模型生成的中心。结果：方法有效识别了模型间的特征（例如Coalescent的树更平衡）。
    2. 运行时间验证： 对n=5~10的树样本，测试算法在不同样本大小下的时间复杂度。结果：在处理几十棵树时可以接受的合理时间（实际的论文展示了对n=10，约10秒内收敛；n=20则进入几分钟级）。

  • 这个例子想说明什么：

    • 验证理论： 他们证明方法能工作——能够处理现实生物数据（时空上大规模采样的序列）产生的树，且在常见的贝叶斯框架中顺利衔接（BEAST输出 -> 无标签化 -> 计算Fréchet均值）。
    • 展示相对baseline的优势？ 他们主要与目视检查比较。过去，要比较加州和佛州树分布，需要拿出所有树看，而本文只需看均值+方差。论文没有与许多同级的baseline方法进行定量比较（因为该领域缺乏这类benchmark）；其优势在于 (a) 将分布总结为一个简洁的可视元素（均值谱系）；(b) 系统性地量化方差（跨样本的可变性程度）。 具体数值结果：加州树的Fréchet方差较高，而德州树和佛州树的方差较小。这暗示了不同地方流行早期的进化历史异质性。

• 🔎 结论是否比证明窄：

  • 窄的地方： 在第5-6节的真实数据和应用部分，作者给出的交叉验证（留一法）结果只展示了“Gurobi MIP可以总是收敛到某个树”，没有提供关于一致性（即样本量m增大时，均值是否会趋向某一特定树）的严格证明。理论部分（定理8）只保证了存在性和组合上求解策略，但没有给出一个统计上的收敛率（如n→∞时，均值的渐近分布）。这个“能否进行假设检验”的问题在introduction里作为动机提及，但在结论部分弱化，作者只是写了“这与使用测试来推断是新方向（future work）”。因此，理论收敛性和推断是比在定理设定下证明的窄——在计算上可行，但在统计显著性的意义上比初始声明窄了。

四、开放问题（点到为止，扎根具体语句）

1. 中心极限定理： “Theorem 8 以及 Multivariate CLT [Hogg et al, 2019, Thm 5.4.4] 可用来检验随机样本是否遵循标准共祖分布……”（Introduction末尾）。但在论文中没有推导或证明Fréchet均值的渐近分布。该所暗示的“样本方差与均值的CLT”仍是一个开放问题——尤其因为Fréchet均值路径测地线可能不是真正的欧式均值，其渐近下界理论仍未给出。扎根：Section 7, Future Work: “Fréchet variance analysis of random objects as in Dubey and Müller (2019) is subject of future research。”

2. 对其他树类型的泛化： 全文假定树是有根、且二分的。但真实的系统发育树经常包括多分支（polytomies）或多物种（多个分支同时形成）。“将本文的算法扩展至允许polytomies（多分支）的树空间”未被处理。扎根：Section 2，定义部分明确为“binary tree（二叉树）”，正文未讨论多分支情况。

3. 计算算法的理论保证： MIP求解器Gurobi只是被当成一个“黑箱”使用。本文没有给出一个紧的（tight）上界（例如，树样本数m和叶子数n与求解时间的多项式关系）。理论分析表明“对于n<=10，可以枚举所有形状”，但10以上没有理论保证。扎根：Section 5, 图5所示对n=20，运行时间变得非常长（与树型选择有关），“我们的MIP求解器收敛到最优解的时限无法保证”。

---

Maintained by 陈星宇 · Homepage · Source on GitHub

评论