B 树—数据—假设

树—数据—假设

Dec 10, 2024

PGLS 最像一面镜子：它不一定告诉你新东西，但它很擅长把你原来不愿意看的细节照出来。

比如这些：

  • 你的物种名到底有没有对齐
  • 你的分支长度代表什么
  • 你有没有在无意中让几条异常点替你"完成论证"
  • 你拿到的是一棵树，还是一堆不确定的树

我不是在吓人。我只是说：PGLS 的麻烦，大部分不在数学里，在整理里。

这篇我想写成一个"用前自检 + 用后诊断"的清单，但还是按随笔的方式写：因为真正的坑，往往不是你不会做，是你没想到它会坑你。

先把树和数据对齐

这一步没做好，后面都是幻觉。

别信你的眼睛

最常见的错误是：树里叫 Homo_sapiens，表里叫 Homo sapiens。
你以为只是一个空格，代码会很礼貌地替你处理。
但它经常不会，它只会很安静地把样本丢掉。

我建议你在任何拟合之前，强制做三件事：

  • 显式列出：树里有哪些 tip label，数据里有哪些 species
  • 计算交集与差集，并把差集打印出来
  • 最终拟合使用的物种数写进结果里（别让它悄悄变化）

时间，还是替换数

你当然可以说"我用的就是这棵树"。
但 PGLS 会追问：你这棵树的分支长度是什么意思？

  • 时间树（ultrametric）更直观：共享历史 = 共享时间
  • 替换数树有时也能用，但解释要更谨慎：共享历史更像"共享变化量"

（你不写清楚，读者很难判断你 $\Sigma$ 的意义是什么。）

一棵树，还是一堆树

很多时候我们只有一棵"最佳树"，于是就拿它当唯一真相。

但如果你的树来自贝叶斯后验或 bootstrap，你其实已经有一堆树了。
那你至少可以做一件很便宜的事：在多棵树上重复拟合，看看结论是否稳定。

结论不稳定并不可耻。可耻的是你明明可以检查，却假装不知道。

残差里藏了什么

PGLS 的核心假设，是误差结构与树一致（或经 $\lambda$ 调过后足够一致）。
所以你需要看的不是只有系数和 p 值，还有残差。

我通常会做这几件事：

  • 残差 vs 拟合值：有没有明显模式（非线性、异方差）
  • 影响点：有没有某个物种把斜率拽得太狠
  • 系统发育信号：残差里是否仍有明显系统发育结构（意味着模型没吃掉你以为它能吃掉的那部分）

这一步的心态很重要：诊断不是为了证明你对，而是为了找出你错在哪儿。

残差诊断的统计学基础

残差的定义

在 GLS/PGLS 下，残差定义为：

$$e = y - X\hat{\beta}_{GLS}$$

在 OLS 下，残差方差是常数（$\sigma^2$）；但在 GLS 下，由于 $\Sigma$ 不是对角矩阵，原始残差 $e$ 的方差不是常数——离根越近的物种的残差，通常方差更大。

标准化残差

为了使残差可比较，需要标准化。一个常见的形式是：

$$e_{std} = (I - H)^{-1/2}e$$

其中 $H$ 是帽子矩阵（hat matrix）在 GLS 形式下的推广：

$$H = X(X'\Sigma^{-1}X)^{-1}X'\Sigma^{-1}$$

$H$ 的对角元 $h_{ii}$ 叫做杠杆值（leverage），反映第 $i$ 个观测对拟合的影响程度。

如果 $h_{ii}$ 很大（接近 1），说明这个物种的性状值对拟合影响很大——在 PGLS 里，这往往意味着它的相关结构很"孤立"（要么特别近缘、要么特别远缘）。

标准化残差 $e_{std}$ 在正确指定的模型下，应该近似 $N(0, \sigma^2)$。如果它的分布明显偏离这个假设（比如厚尾、偏斜），就说明模型可能有问题。

影响点分析：Cook’s Distance

Cook’s Distance 是衡量第 $i$ 个观测对回归系数影响的经典指标。在 GLS 下，它的形式是：

$$D_i = \frac{e_i^2}{p \cdot \hat{\sigma}^2} \cdot \frac{h_{ii}}{(1 - h_{ii})^2}$$

其中 $p$ 是参数个数（包括截距），$\hat{\sigma}^2$ 是误差方差的估计。

直觉是：$D_i$ 越大，说明删掉第 $i$ 个观测后，$\hat{\beta}$ 的变化越大。

一个常用的阈值是 $D_i > 4/n$ 或 $D_i > 1$，但这些阈值在 PGLS 里不一定适用——你可以用它们作为参考，但更重要的是看哪些点明显偏离其余点。

残差的系统发育信号检验

做完 PGLS 后，你还需要检查：残差本身是否还有系统发育结构。

如果残差的系统发育信号很强（比如用 physig 或 ababntest 检验显著），说明模型没有完全吃掉系统发育相关性——可能有遗漏变量、错误的树、或者 $\lambda$ 设得不对。

踩过的坑

把"控制系统发育"理解成"消除生物学意义"

PGLS 不是在抹掉历史，它是在阻止历史冒充因果。
你依然可以讲生物学故事，只是故事要更诚实：你要区分"共享祖先造成的相似"与"自变量解释的变化"。

把 λ 当成"模型越复杂越好"的装饰

$\lambda$ 不是 KPI。它不是越大越好，也不是越小越现代。
它只是一个态度：你愿意让树在残差结构里说多大声。

你应该做的是报告它、解释它、并检查结果对它是否敏感。

λ 敏感性的检验

似然比检验（LR Test）

如果你的 null hypothesis 是 $\lambda = 1$（完全按 BM 的相关结构来），alternative hypothesis 是 $\lambda$ 自由估计（$\lambda \in [0, 1]$），可以用似然比检验：

$$LRT = 2(\ell_{\hat{\lambda}} - \ell_{\lambda=1}) \sim \chi^2_1$$

其中 $\ell_{\hat{\lambda}}$ 是 $\lambda$ 在最优估计值处的对数似然，$\ell_{\lambda=1}$ 是固定 $\lambda = 1$ 时的对数似然。

这个检验用来判断：允许 $\lambda$ 自由估计，是否显著改善了模型拟合。

边界问题

一个常见的陷阱：当 $\lambda$ 的最优估计接近 0 时（或者正好是 0），上面的 LRT 不适用了。

原因是：$\lambda = 0$ 位于参数空间 $[0, 1]$ 的边界上。在这种"边界约束"的情况下，似然比统计量不再服从标准的 $\chi^2_1$ 分布——它的分布被"截断"了。

一个实用的解决方案：

  • 不依赖 LRT p 值：直接报告 $\lambda$ 的点估计和 95% 置信区间
  • 用轮廓似然法（profile likelihood）画 $\lambda$ 的置信区间
  • 如果置信区间很宽，说明 $\lambda$ 本身估计不稳定，结论对 $\lambda$ 敏感

敏感性分析的实践

做敏感性分析的标准流程：

  1. 在最优 $\hat{\lambda}$ 处报告主结果
  2. 额外报告 $\lambda = 0$ 和 $\lambda = 1$ 时的结果（或者在 $\lambda$ 的 95% 置信区间两端）
  3. 比较：$\beta$ 的方向是否改变、p 值是否跨越显著性阈值

如果上述任意一项对 $\lambda$ 敏感，你的结论应该更保守——或者说，你应该承认"结论对系统发育相关结构的假设是敏感的"。

（这种事写进论文里不好看，但总比发表后被人挑出来强。）

用一棵树给出一个结论，然后忘了树本身也在漂

树来自数据，数据也会变；比对参数会变，建树方法会变；你今天的"最佳树"不一定是明天的"最佳树"。

如果你的结论对树特别敏感，那结论就应该写得更谨慎——哪怕这很不爽。

最少但必须

我更喜欢把这部分写成"最少但必须"的几行：

  • 使用的树来源与构建方法（或引用）
  • 分支长度含义（时间 / 替换数 / 处理方式）
  • PGLS 里使用的相关结构与演化模型（BM / λ 等）
  • 最终纳入分析的物种数（以及缺失处理）
  • 是否做了树不确定性的敏感性分析（做了就一句话说明怎么做，没做也别装作做了）

这些话写进去，不会让你的结果更显著，但会让它更可信。

很多人把 PGLS 当成一种"更高级的回归"，我更愿意把它当成一种提醒：

PGLS 不会替你完成推理，它只会逼你把"树—数据—假设"对齐。

对齐之后，你再去讲故事，才不那么心虚。

附录：“永远别偷懒"的提醒

反复用的骨架：读入树和性状表、对齐物种、拟合、做几个基本诊断。

library(ape)
library(nlme)

# 1) 读入
tree <- read.tree("tree.nwk")
dat <- read.csv("trait.csv", stringsAsFactors = FALSE)

# 假设 dat 里有 species, x, y
rownames(dat) <- dat$species

# 2) 对齐
sp_common <- intersect(tree$tip.label, rownames(dat))
tree2 <- drop.tip(tree, setdiff(tree$tip.label, sp_common))
dat2 <- dat[sp_common, , drop = FALSE]

# 3) 拟合（相关结构来自树）
fit <- gls(y ~ x, data = dat2,
           correlation = corBrownian(phy = tree2),
           method = "REML")
summary(fit)

# 4) 基本诊断
par(mfrow = c(1, 2))
plot(fitted(fit), resid(fit), xlab = "Fitted", ylab = "Residual")
abline(h = 0, lty = 2)
qqnorm(resid(fit)); qqline(resid(fit))

提取 $\Sigma$ 矩阵并检查条件数

协方差矩阵的条件数（condition number）反映它的数值稳定性。如果条件数很大（比如 > 1000），矩阵求逆可能不稳健。

library(ape)

tree <- read.tree("tree.nwk")
Sigma <- vcv(tree)  # 提取协方差矩阵

cat("Sigma 矩阵维度:", dim(Sigma), "\n")
cat("条件数:", kappa(Sigma), "\n")
cat("行列式:", det(Sigma), "\n")

# 热图可视化
heatmap(Sigma, symm = TRUE, margin = c(8, 8))

λ 的似然剖面图

library(ape)
library(nlme)

tree <- rtree(30)
x <- rnorm(30)
y <- 0.5 * x + rTraitCont(tree)[tree$tip.label]

# 在不同 λ 值下计算对数似然
lambda_seq <- seq(0, 1, by = 0.05)
loglik <- sapply(lambda_seq, function(lam) {
  fit <- try(gls(y ~ x,
                  correlation = corPagel(lam, phy = tree),
                  method = "REML"),
              silent = TRUE)
  if (inherits(fit, "gls")) logLik(fit) else NA
})

plot(lambda_seq, loglik, type = "b",
     xlab = expression(lambda),
     ylab = "Log-likelihood",
     main = "Likelihood profile for lambda")
abline(v = lambda_seq[which.max(loglik)], lty = 2)

Cook’s Distance 计算

library(ape)
library(nlme)

tree <- rtree(30)
x <- rnorm(30)
y <- 0.5 * x + rTraitCont(tree)[tree$tip.label]

fit <- gls(y ~ x, correlation = corBrownian(phy = tree), method = "REML")

# 提取必要元素
n <- length(y)
p <- length(coef(fit))
sigma2 <- fit$sigma^2

# 计算帽子矩阵对角元 h_ii
X <- model.matrix(fit)
Sigma_inv <- solve(vcv(tree))
H <- X %*% solve(t(X) %*% Sigma_inv %*% X) %*% t(X) %*% Sigma_inv
h_ii <- diag(H)

# 计算标准化残差
residuals_std <- resid(fit) / sqrt(sigma2 * (1 - h_ii))

# Cook's Distance
D <- (residuals_std^2 / p) * (h_ii / (1 - h_ii)^2)

# 画图
plot(D, type = "h",
     xlab = "Species index",
     ylab = "Cook's Distance",
     main = "Cook's Distance for each observation")
abline(h = 4 / n, col = "red", lty = 2)  # 常用阈值
TouchingFish.top