不爱解释的人 / Phylogenetic GLS

同一份数据，同一个自变量，同一个因变量。

OLS 说：显著。

PGLS 说：不显著。

很想把锅甩给树：是不是树不对？是不是分支长度乱了？

还是哪里写错了代码？

PGLS 只是站在那儿，像一个不爱解释的人。

这篇我想用一种不那么"公式"的方式，讲清楚 PGLS 到底在做什么。

你只需要抓住一个中心句：

PGLS 不是换了一个更玄的回归，它只是承认误差有相关结构。

误差像一盘散沙

普通线性回归（OLS，Ordinary Least Squares，普通最小二乘法）通常默认：

这里 $I$ 的意思很朴素：误差之间互不相关，方差还都一样。

你可以把它理解成：每个物种的"解释不掉的那部分"，互相之间不应该传染。

但系统发育数据里，这个默认经常站不住脚：近缘物种共享历史，解释不掉的那部分会一起漂。

于是你进入另一个世界。

OLS 估计量的完整推导

$y$ 是因变量的向量（每个物种一个值），$X$ 是设计矩阵（每行一个物种，每列一个自变量，第一列通常是全 1，代表截距），$\beta$ 是我们要估计的系数向量。符号 $'$ 表示转置。

OLS 想做的事情很简单：找一组 $\beta$，让预测值 $X\beta$ 和观测值 $y$ 之间的残差平方和最小。目标函数写成矩阵形式就是：

把它展开，是为了后面能对 $\beta$ 求导：

中间两项看起来不一样，其实 $y'X\beta$ 是一个标量（1×1 的数），标量的转置等于自己，所以 $\beta'X'y = (y'X\beta)' = y'X\beta$。合并之后：

接下来对 $\beta$ 求导。这里用到两条向量求导规则：$\partial(\beta' a)/\partial\beta = a$（线性项的导数就是系数本身），$\partial(\beta' A\beta)/\partial\beta = (A + A')\beta$（二次项的导数，当 $A$ 对称时简化为 $2A\beta$）。代入：

令导数为零，得到正规方程（normal equation）：

两边左乘 $(X'X)^{-1}$（假设 $X$ 列满秩，所以 $X'X$ 可逆）：

这就是 OLS 的解。它的方差是：

这个方差的来源：$\beta_{OLS} = (X'X)^{-1}X'y$ 是 $y$ 的线性函数，而 $y$ 的协方差是 $\sigma^2 I$，线性变换的协方差会传递，于是

（这套推导的关键在于：OLS 假设误差协方差是 $\sigma^2 I$，所以它只关心残差的总平方，不关心谁和谁相关。）

误差有了骨架

GLS（Generalized Least Squares，广义最小二乘法）把假设改成：

这张 $\Sigma$ 不是凭空捏出来的。做 PGLS 时，它通常来自树（以及你选择的性状演化模型，比如 BM，再加上可能的 $\lambda$）。

所以 PGLS 的意思其实是：

同样是回归，我只是把你样本之间的"亲疏关系"算进了误差结构。

换一个坐标系再看

理解 PGLS 最好用的一个比喻，是"换坐标系"。

如果误差有相关性，你可以想象在原空间里它们拉着彼此的手，不容易看清每个人的真实波动。
GLS 做的事情，等价于把数据做一个线性变换，把相关性"抹平"，让新的误差更接近独立同分布——也就是把手松开。

这就是很多教材里说的 whitening（白化）直觉：
在变换后的世界里，OLS 才变得合理。

所以你会看到一种常见现象：

• 系数估计（$\beta$）可能变化不大
• 但标准误（SE）会变，进而 p 值会变

也就是说，PGLS 经常不是在跟你争"方向"，而是在跟你争"你到底有多少有效信息"。

GLS 估计量的完整推导

GLS 的问题设置是：

其中 $\Sigma \neq I$（样本之间有相关结构）。

前面说"换坐标系"，现在来兑现。核心思想是：找到一个变换矩阵 $P$，使得 $P\epsilon$ 的协方差变成 $\sigma^2 I$——把拉着手的误差松开，然后在变换后的世界里安心用 OLS。

步骤 1：构造变换矩阵 $P$

$\Sigma$ 是对称正定矩阵（协方差矩阵天然如此），所以它可以做 Cholesky 分解——你可以把它理解成"把一个方块拆成两个三角形的乘积"：

其中 $L$ 是下三角矩阵。令 $P = L^{-1}$，也就是 $L$ 的逆矩阵。

验证白化效果——把 $P$ 作用到 $\Sigma$ 上：

中间两步分别是 $L^{-1}L = I$ 和 $L'L'^{-1} = I$。结果就是：变换后的误差协方差变成了单位矩阵，手松开了。

步骤 2：变换模型

对原模型 $y = X\beta + \epsilon$ 两边左乘 $P$：

记 $y^* = Py$，$X^* = PX$，$\epsilon^* = P\epsilon$，模型变成：

其中 $\epsilon^* \sim N(0, \sigma^2 I)$。注意 $\beta$ 没变——我们换的是坐标系，不是要估计的东西。

步骤 3：在变换后的世界用 OLS

既然 $\epsilon^*$ 已经是独立同分布了，直接套 OLS 公式：

把 $X^* = PX$ 和 $y^* = Py$ 代回去（注意 $(PX)' = X'P'$）：

这里需要一个关键性质：$P'P = \Sigma^{-1}$。

推导：从白化条件 $P\Sigma P' = I$ 出发，右乘 $P'^{-1}$ 得 $P\Sigma = P'^{-1}$，再左乘 $P'$ 得 $P'P\Sigma = I$，所以 $P'P = \Sigma^{-1}$。

代入：

这就是 GLS 的闭式解。和 OLS 的 $\beta_{OLS} = (X'X)^{-1}X'y$ 对比，唯一的区别就是 $\Sigma^{-1}$ 出现在了所有该出现的位置——它就是那个"去相关的权重"。

GLS 估计量的方差：

对比 OLS 的方差 $Var(\beta_{OLS}) = \sigma^2(X'X)^{-1}$，差异来自 $\Sigma^{-1}$ 对 $X'X$ 的"加权"作用。

你以为有 N 个样本

当样本之间高度相关时，“有效样本量"会下降。

你表面上有 $n$ 个物种，但其中很多物种共享很长一段历史，它们提供的信息不是独立新增的。
OLS 假装这些信息都独立，于是它往往会低估不确定性，让你更容易显著。

PGLS 把相关结构纳入后，常见的结果就是：不确定性变大，显著性变弱。

这不是"保守”，这是"把历史算进了账"。

为什么 SE 会变

OLS 和 GLS 的方差对比：

当样本独立时（$\Sigma = I$），两者完全一样。
当样本相关时（$\Sigma$ 有较大的非对角元），$\Sigma^{-1}$ 的效果是"惩罚"那些来自相关样本的信息。

直觉是：如果物种 A 和 B 高度相关，它们俩合在一起提供的新信息量，小于两个独立物种的信息量之和。
$\Sigma^{-1}$ 就像一个"去相关"的权重，把重复计算的信息压下去。

结果就是 $(X'\Sigma^{-1}X)^{-1}$ 的对角元（也就是 $\beta$ 的方差）变大。
方差变大 → SE 变大 → t 值变小 → p 值变大。

这就是为什么 PGLS 经常让你"没那么显著"：不是因为你的效应方向变了，而是因为你的有效信息量被你高估了。

（如果你曾经在 OLS 里看到漂亮的小 p 值，在 PGLS 里突然变大了——现在你至少知道它去哪儿了。）

历史不替 x 背锅

很多文章会写一句话：we controlled for phylogeny。

我一直觉得这句话很危险，因为它容易被理解成：我把生物学消掉了。
但 PGLS 真正控制的不是生物学，而是残差结构：它不让"共同祖先造成的相似"悄悄冒充成"自变量的解释力"。

你依然在问同一个问题：$x$ 和 $y$ 的关系是什么。
只不过你现在更谨慎：你不允许历史替 $x$ 偷走功劳。

写在最后

PGLS 最讨人嫌的一点在于：它经常让你没那么显著。
但我后来慢慢能接受这件事，是因为我意识到——显著性本来就不该那么廉价。

当你研究的是有家谱的对象，你就得承认：它们之间的相似，不完全是你模型里那根 $x$ 的功劳。

PGLS 只是把这种不舒服，写进了标准误里。

附录

R 示例：PGLS 与 OLS 的对比

这段代码不追求完整的 PGLS 工程化流程，只想让你看见：

1. 同样的回归式子
2. 在独立误差 vs 相关误差的假设下，会出现不同的 SE / p 值

library(ape)
library(nlme)

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

# 造一点数据：x 独立生成；y 带一点系统发育相关的误差
x <- rnorm(40)
names(x) <- tree$tip.label

e <- rTraitCont(tree, model = "BM")
e <- e[names(x)]

y <- 0.5 * x + scale(e)[, 1]
dat <- data.frame(species = names(x), x = x, y = y)
rownames(dat) <- dat$species

# OLS
summary(lm(y ~ x, data = dat))

# 相关误差（近似 PGLS 的核心：误差相关结构由树给出）
fit_gls <- gls(y ~ x, data = dat,
               correlation = corBrownian(phy = tree),
               method = "REML")
summary(fit_gls)

Python 示例：$(X'X)^{-1}$ vs $(X'\Sigma^{-1}X)^{-1}$ 的数值对比

下面用一个小例子直观展示 OLS 和 GLS 方差矩阵的差异：

import numpy as np

np.random.seed(1)

# 构造一个强相关的协方差矩阵（模拟近缘物种）
Sigma = np.array([
    [1.0, 0.9, 0.8, 0.0],
    [0.9, 1.0, 0.85, 0.0],
    [0.8, 0.85, 1.0, 0.0],
    [0.0, 0.0, 0.0, 1.0]
])

# 设计矩阵 X（简单：截距 + 一个协变量）
X = np.column_stack([
    np.ones(4),          # 截距
    [1, 2, 3, 4]         # 协变量 x
])

# OLS 的方差矩阵
XtX_inv = np.linalg.inv(X.T @ X)
print("OLS Var(beta):")
print(np.round(XtX_inv, 4))

# GLS 的方差矩阵
Sigma_inv = np.linalg.inv(Sigma)
XtSigmaX_inv = np.linalg.inv(X.T @ Sigma_inv @ X)
print("\nGLS Var(beta):")
print(np.round(XtSigmaX_inv, 4))

# 观察：对角元变大了
print("\n放大倍数 (GLS/OLS):")
print(np.round(np.diag(XtSigmaX_inv) / np.diag(XtX_inv), 2))

你会看到 GLS 的方差对角元比 OLS 大——特别是当样本相关性更强的时候。

Cholesky 分解的手动计算

Cholesky 分解 $\Sigma = LL'$ 是 GLS “白化” 变换的核心。下面是一个 3×3 矩阵的手动分解：

import numpy as np

Sigma = np.array([
    [2.0, 1.0, 0.5],
    [1.0, 2.0, 0.8],
    [0.5, 0.8, 1.5]
])

# 检查对称正定性
print("特征值:", np.linalg.eigvalsh(Sigma))  # 应该全部为正

# Cholesky 分解
L = np.linalg.cholesky(Sigma)
print("\nL (下三角):")
print(np.round(L, 3))

# 验证：L @ L.T = Sigma
print("\nL @ L.T:")
print(np.round(L @ L.T, 3))

# 白化变换矩阵 P = L^{-1}
P = np.linalg.inv(L)
print("\nP = L^{-1}:")
print(np.round(P, 3))

# 验证白化效果
print("\nP @ Sigma @ P.T (应该接近单位矩阵):")
print(np.round(P @ Sigma @ P.T, 3))

Cholesky 分解的好处是：$L$ 是下三角的，求逆比一般的对称矩阵快得多。在实际的大规模 PGLS 计算中，这个分解是效率的关键。