一棵树，怎么长出一个协方差矩阵

我以前对"在树上做统计"最大的抵触，不是我不信树，而是我不信那句很随意的话：

“我们假设性状沿着树做 Brownian motion。”

听起来像一句推脱。像你问一个人为什么迟到，他说：“路上有点堵。”

堵在哪儿？怎么堵的？堵到什么程度？有没有备选路线？——你越问越像在吵架。

但后来我发现，Brownian motion（BM）在这里的意义，和"世界真的在做随机游走"没什么关系。它更像是一种最小的诚实：你承认性状会变、承认变化会积累、承认近缘会更像，然后你就能把"历史"写成一张可以计算的表。

那张表就是协方差矩阵：$\Sigma$。

这篇我想把三件事讲清楚：

1. BM 在系统发育语境下到底是什么（别把它当玄学）
2. 为什么 $\Sigma_{ij}$ 只和"共享祖先到根的那段路"有关
3. Pagel’s $\lambda$ 到底在调什么：不是调"模型拟合"，而是在调"你愿意相信多少历史"

先把 Brownian motion 说成一句人话

BM 在这里可以理解成：

性状的增量（变化）在每一小段时间里都是随机的，方向不固定；但变化会累计，所以时间越长，方差越大。

你不必把它当作宇宙真理。你只要把它当作一个"最低成本"的默认假设：它不要求你知道选择压力、也不要求你知道适应峰值在哪里；它只是说——不解释的东西，就先当成随机漂移。

于是你立刻得到一个很实用的结论：

• 沿着分支走得越久，性状的不确定性越大
• 两个物种如果共享一段历史，它们那段历史里累积的"漂移"是同一份，所以它们会相关

相关性在这里不是"统计技巧"，是"共享账本"。

Brownian Motion 的严格定义

上面的人话版够用了，但如果你想完全理解后面的推导，这里有一份不走捷径的版本。

一维标准 Brownian Motion $W(t)$ 是一个随机过程，满足：

1. 初始条件：$W(0) = 0$

2. 独立增量：对任意 $0 \leq s < t$，增量 $W(t) - W(s)$ 与 $\{W(u): u \leq s\}$ 独立

3. 正态增量：$W(t) - W(s) \sim N(0, \sigma^2(t - s))$

4. 连续路径：$W(t)$ 关于 $t$ 几乎必然连续

第 3 条是核心：增量服从正态分布，方差随时间线性增长。这直接导致"走得越久，方差越大"这个结论。

如果 $\sigma^2 = 1$，就叫标准 Brownian Motion（方差参数为 1）；如果 $\sigma^2 \neq 1$，就是带方差的 Brownian Motion，方差参数会在协方差矩阵里单独出现。

在系统发育的语境下，我们通常关心的是：沿着每条分支，性状演化都是一次独立的 Brownian Motion；分支长度 $t$ 就是这次 Brownian Motion 的时间步长。

协方差的完整推导

有了严格定义，我们可以把 BM 下的协方差算出来。

设物种 $i$ 和 $j$ 的性状值为 $X_i$ 和 $X_j$。把性状值沿树分解：

其中：

• $X_0$ 是根节点的性状值（通常设为常数 0）
• $\Delta_{0 \rightarrow MRCA(i,j)}$ 是从根到最近共同祖先的累积变化（两个物种共享这段）
• $\Delta_{MRCA(i,j) \rightarrow i}$ 是从 MRCA 到物种 $i$ 的独立变化（$j$ 不参与）

由 BM 的独立增量性质，共享段和独立段是独立的。

方差：

由于 $X_0$ 是常数（方差为 0），且每段增量都服从 $N(0, \sigma^2 \cdot t)$：

协方差：

第二项为零（独立增量）；第一项是共享段的方差：

这就是核心结论：协方差等于共享路径长度乘以方差参数。

把它写得更简洁：

对角线 $\Sigma_{ii}$ 对应 $MRCA(i,i) = i$，所以 $t_{0 \rightarrow MRCA} = t_{0 \rightarrow i}$——这和我们刚才算的方差一致。

为什么协方差等于共享路径长度

设想你从根出发，一路走到每个物种的叶子。BM 的性状值可以拆成两部分：

• 从最近共同祖先（MRCA，Most Recent Common Ancestor）的那段：两者共享
• 从 MRCA 分叉之后各自走的那段：各走各的

所以：

• $Var(X_i)$ 随着根到 $i$ 的路径长度增加而增加
• $Cov(X_i, X_j)$ 只来自两者共享的那段路径

如果树的分支长度用"时间"记账，那么"共享路径长度"就真的是共享时间。
如果树的分支长度用"替换数"记账，那么它就是一种"共享进化变化量"的代理。

我很喜欢把 $\Sigma$ 想象成这样一张表：你不是在问"物种 A 和 B 像不像"，你是在问——

你们俩的相似里，有多少是同一段历史造成的。

这张表是可以直接从树算出来的，它就是误差结构的核心。

λ 是一个态度

到这里，很多人的心会突然虚一下：

“好吧，我承认有相关性。但我并不想把一棵树当成完全正确的历史记录。那怎么办？”

这就是 $\lambda$ 出场的方式：它是一种折中。

你可以把它当成一个滑块：

• $\lambda = 0$：你彻底不买账。相关性被关掉，退化成普通 OLS（在误差独立的世界里回归）。
• $\lambda = 1$：你完全照 BM 的相关结构来，相当于"树说什么我就信什么"。
• $0 < \lambda < 1$：你承认有系统发育信号，但它没有强到让整棵树主宰一切。

更具体一点：$\lambda$ 做的事，通常是把 $\Sigma$ 的"非对角元素"按比例缩放，让相关性整体变弱或变强（不同实现细节会略有差异，但直觉是这个直觉）。

所以我更愿意把它理解为：

$\lambda$ 是你对"树解释力"的态度。

（这句话听起来有点哲学，但它在写论文时很实用：你至少知道自己在说什么。）

$\lambda$ 的矩阵构造与估计

$\Sigma_\lambda$ 的构造

Pagel’s $\lambda$ 对原始协方差矩阵 $\Sigma$ 的改造方式，最常见的定义是：

• 对角线保持不变：$(\Sigma_\lambda)_{ii} = \Sigma_{ii}$

• 非对角线按 $\lambda$ 缩放：$(\Sigma_\lambda)_{ij} = \lambda \cdot \Sigma_{ij}$，其中 $i \neq j$

这等价于：

其中 $D$ 是 $\Sigma$ 的对角线矩阵（只保留对角元，其余置零）。

另一种常见表达（需要验证等价性）是：

直觉是：$\lambda = 0$ 时 $\Sigma_\lambda = D$（只有自己的方差，没有相关性）；$\lambda = 1$ 时 $\Sigma_\lambda = \Sigma$（完全按树来）。

最大化似然函数估计 $\lambda$

在实际拟合中，$\lambda$ 是通过 MLE（Maximum Likelihood Estimation，最大化似然函数）来估计的。

模型是：

对给定的 $\lambda$，GLS 估计量是：

对应的对数似然函数（忽略常数项）是：

在实际操作中，通常用数值优化（如 optimize() 或 nlminb()）来找使 $\ell(\lambda)$ 最大化的 $\lambda$。这个过程往往受到 $\lambda \in [0, 1]$ 的约束。

一个常见的问题：当 $\lambda = 0$ 时，$\Sigma_\lambda$ 在边界上（对角矩阵），似然比检验（likelihood ratio test）不能直接用标准的 $\chi^2_1$ 分布。在这种情况下，建议用参数边界上的校正，或者直接报告 $\lambda$ 的置信区间而不是做假设检验。

这一切跟 PGLS 的关系

如果你把性状（或回归残差）看成 BM，那么你自然就得到一个由树导出的协方差矩阵 $\Sigma$。
而 GLS / PGLS 做的事情，就是：在 $\epsilon \sim N(0, \sigma^2 \Sigma)$ 的世界里做回归。

下一篇我们就直接进这个房间：PGLS 到底在拟合什么，它为什么会改变标准误，为什么它经常让你"没那么显著"。

BM 有点像统计学里那种"你先别问太多"的假设。它粗糙、也可能不对。但它把一件最关键的事说清楚了：

历史会累积，累积会制造相关。

至于你要不要完全相信它——这正是 $\lambda$ 这种东西存在的理由：它不是替你决定答案，它只是把"犹豫"也变成了一个可估计的量。

附录

R 示例：用一棵树生成 $\Sigma$（树隐含的协方差结构）

1. 用一棵树生成 $\Sigma$（树隐含的协方差结构）
2. 用相关矩阵的热图看"谁和谁更像"这件事不是均匀随机的

library(ape)

set.seed(1)
tree <- rtree(30)

Sigma <- vcv(tree)      # 树导出的协方差矩阵（与分支长度有关）
R <- cov2cor(Sigma)     # 相关矩阵更直观

heatmap(R, symm = TRUE, margins = c(6, 6))

Python 示例：手动构建 3 物种的 $\Sigma$ 矩阵

用 Python（配合 numpy 和 ete3 或 toytree）可以更清楚地看到矩阵是怎么算出来的。下面是一个手动的版本：

import numpy as np

# 定义一棵3物种树（简化表示）
# 树结构：
#       root
#      /    \
#     A      B
#           /
#          C
# 分支长度：root->A=1, root->B=1, B->C=0.5
# （实际使用时用ape或ete3读取.nwk文件）

# 手动计算协方差矩阵 Sigma
# Sigma[i,j] = 从根到 i 和 j 的 MRCA 的路径长度
# 约定顺序：A, B, C

# 设 sigma^2 = 1
sigma2 = 1.0

# 分支长度
t_root_A = 1.0      # root 到 A
t_root_B = 1.0      # root 到 B
t_B_C = 0.5         # B 到 C
t_root_C = t_root_B + t_B_C  # root 到 C = 1.5

# MRCA关系：
# A-B: MRCA = root, 共享路径 = 0
# A-C: MRCA = root, 共享路径 = 0
# B-C: MRCA = B, 共享路径 = root->B = 1.0

Sigma = np.array([
    [sigma2 * t_root_A, sigma2 * 0,         sigma2 * 0       ],  # A
    [sigma2 * 0,        sigma2 * t_root_C,   sigma2 * t_root_B],  # B
    [sigma2 * 0,        sigma2 * t_root_B,   sigma2 * (t_root_C)]   # C
])

print("Sigma (协方差矩阵):")
print(Sigma)

# 相关矩阵
D_inv = np.diag(1.0 / np.sqrt(np.diag(Sigma)))
R = D_inv @ Sigma @ D_inv
print("\n相关矩阵 (cov2cor):")
print(np.round(R, 3))

$\lambda$ 如何改变 $\Sigma$

下面展示 $\lambda$ 的效果：把非对角线按 $\lambda$ 缩放，相关性整体降低。

import numpy as np

Sigma = np.array([
    [2.0, 1.5, 0.5],
    [1.5, 3.0, 1.2],
    [0.5, 1.2, 2.5]
])

# 对角线矩阵
D = np.diag(np.diag(Sigma))

for lam in [0, 0.5, 1]:
    Sigma_lam = D + lam * (Sigma - D)
    R_lam = Sigma_lam / np.sqrt(np.outer(np.diag(Sigma_lam), np.diag(Sigma_lam)))
    print(f"\nλ = {lam}:")
    print(np.round(R_lam, 3))

你会看到：

• $\lambda = 0$：完全对角（无相关）
• $\lambda = 0.5$：相关性减半
• $\lambda = 1$：原始相关性

这种可视化有助于理解为什么 $\lambda$ 是一个"滑块"：它不是在重塑相关性的结构，只是在调整它的强度。