随机效应模型

第 4 章给出 FE-IV 模型的合并 $\mathrm{RR} = 0.65$，但 $I^2 = 92\%$ 表明 13 项研究并不在估同一个真实效应。本章松开"同一个真实效应"这条假设，把研究间的真实差异用一个新的方差参数 $\tau^2$ 估计出来。BCG 数据是异质性显著的典型案例，随机效应模型给出的合并结果与 FE 差异巨大，是 meta 分析教学的标准范例。

从固定效应到随机效应

固定效应模型假设 $\hat{y}_i = \theta + \varepsilon_i$，所有研究共享同一个真实效应 $\theta$。随机效应模型把这个假设松开：每项研究自己有一个真实效应 $\theta_i$，13 个 $\theta_i$ 不再要求相等，而是从一个总体分布里抽出来的。

把 $\tau^2$ 翻译成生活语言。假设 13 项 BCG 试验各自的真实效应分别是 $\theta_1, \theta_2, \ldots, \theta_{13}$。固定效应模型说这 13 个数都相等（共有一个 $\theta$）；随机效应模型说这 13 个数不相等，但围绕一个总体均值 $\mu$ 分布，分布的宽度就是 $\tau$。BCG 数据里，伦敦 Hart & Sutherland 1977 的真实效应可能是 $\mathrm{log\ RR} = -1.4$，赤道附近 TPT Madras 1980 的真实效应可能是 $0$，二者本来就不一样。$\tau^2$ 把这种本来就有的差异量化出来，让合并 CI 老老实实变宽。

随机效应模型的权重

随机效应模型的权重也按方差倒数计算，但方差包含两部分：研究内的抽样方差 $v_i$ 和研究间的真实方差 $\tau^2$。

把 $\tau^2$ 加进权重的分母会带来两个直接后果。其中一个：所有研究的权重都比 FE 模型小，因为分母变大了，导致合并 CI 比 FE 宽。另一个：原本权重悬殊的研究权重被"摊平"。FE 模型下 TPT Madras 一项就占 41%，到 RE 模型下这个比例会显著缩水，因为 $\tau^2$ 把它的相对优势削弱了。

DerSimonian-Laird 估计量

$\tau^2$ 不是直接观测的，必须从数据里估计。DerSimonian & Laird 1986 提出了最早也是最常用的 $\tau^2$ 矩估计量。

在 BCG 数据上手算一遍。先把第 4 章给出的 13 项研究 IV 权重 $w_i = 1/v_i$ 全部累加：$\sum w_i = 609.70$。再算 $\sum w_i^2 = 94{,}820.58$，所以 $\sum w_i^2 / \sum w_i = 94{,}820.58 / 609.70 = 155.52$。DL 公式分母 $\sum w_i - \sum w_i^2 / \sum w_i = 609.70 - 155.52 = 454.18$。Q 已经在第 4 章算出 $Q = 152.23$，自由度 $k - 1 = 12$，分子 $Q - 12 = 140.23$。代入：

分子 140.23 是观察到的研究间散布超出抽样波动的部分，分母 454.18 是把这个超额量还原到 log RR 平方尺度的换算因子。

library(metafor)
d  <- read.csv("data/bcg.csv")
es <- escalc(measure = "RR", ai = tpos, bi = tneg,
             ci = cpos, di = cneg, data = d)

# 手算 DL τ² 的中间量
w       <- 1 / es$vi
sum_w   <- sum(w)              # 609.70
sum_w2  <- sum(w^2)            # 94820.58
denom   <- sum_w - sum_w2 / sum_w   # 454.18
Q       <- rma(yi, vi, data = es, method = "EE")$QE   # 152.23
tau2_DL <- (Q - 12) / denom
cat("tau^2 (DL, 手算) =", round(tau2_DL, 4), "\n")

# DerSimonian-Laird τ² 估计（与手算应当一致）
re_dl <- rma(yi, vi, data = es, method = "DL")
print(re_dl)

DL 是矩估计——它要求 $Q$ 的期望值在 FE 模型下等于 $k - 1$，超出 $k - 1$ 的部分就归因于研究间方差。这个推导不需要假设 $\theta_i$ 服从正态，是一种"分布无关"的估计。代价是 DL 在小研究数（$k < 20$）和高异质性下偏低：它会低估 $\tau^2$，进而高估权重精度，让 CI 偏窄。

REML 估计

REML（Restricted Maximum Likelihood，限制性最大似然）是基于正态假设的最大似然估计，理论上效率更高，在小研究数下偏倚比 DL 小。Cochrane Handbook 与 metafor 包都把 REML 作为默认 $\tau^2$ 估计量。

REML 的优势在小样本下尤为明显。DL 的矩估计在 $k$ 小的时候 $Q$ 本身的方差就不可忽略，DL 估计波动很大。REML 用整个似然函数的形状信息，对 $\tau^2$ 的估计更稳。代价是 REML 假设 $\theta_i$ 服从正态，DL 不假设。实务上正态假设的违反在 $k$ 较大时影响有限。

re_reml <- rma(yi, vi, data = es, method = "REML")
print(re_reml)

Hartung-Knapp 调整

DL 与 REML 都用大样本正态近似算合并 CI——把 $\hat{\mu}$ 的方差当作正态分布的方差，乘以 1.96 算 95% CI 半宽。这个近似在 $k$ 大时合理，$k$ 少时偏窄。Hartung & Knapp 2001 提出用 $t$ 分布替代正态、自由度 $k - 1$，让 CI 更稳健。

HK 为什么更稳？ 两件事都是把 CI 拉宽。第一件事是用每项研究偏离合并均值的实际散布替换理论方差，让 CI 反映"实际看到的散布"而不是"理论应该有的散布"。第二件事是把 $z_{0.975} = 1.96$ 换成 $t_{k - 1, 0.975}$。当 $k = 13$，$t_{12, 0.975} = 2.179$，比 1.96 大约 11%；$k = 5$ 时 $t_{4, 0.975} = 2.776$，比 1.96 大约 42%。研究数越少，HK 调整把 CI 拉得越宽，避免小样本下的过度自信。

re_hk <- rma(yi, vi, data = es, method = "REML", test = "knha")
print(re_hk)

预测区间 PI

随机效应模型有一个比 CI 更直接对应临床决策的指标——预测区间（prediction interval，简称 PI）。CI 描述合并均值的不确定性，PI 描述未来一项新研究的真实效应可能在哪里。在异质性显著时 PI 比 CI 重要得多。

PI 与 CI 的差别在临床上至关重要。CI 告诉你的是"全球 BCG 的平均效应"；PI 告诉你的是"在某一个具体地区开展 BCG 项目，预期效应的范围"。如果 PI 跨过零线，意味着未来一项研究完全可能给出无效甚至反向的结果，临床推广就要谨慎。

pi_hk <- predict(re_hk)
cat("合并 RR =", round(exp(pi_hk$pred), 3),
    "  95% CI =", round(exp(pi_hk$ci.lb), 3), "to",
                  round(exp(pi_hk$ci.ub), 3), "\n")
cat("95% PI  = [", round(exp(pi_hk$pi.lb), 3), ",",
                   round(exp(pi_hk$pi.ub), 3), "]\n")

Forest plot 带预测区间

随机效应 forest plot 通常在合并菱形下方加一根额外的横线代表 PI，方便临床读者直接读出"未来研究的可能效应范围"。metafor 通过参数 addpred = TRUE 自动添加这一行。

第 5 章到这里收束。三种 RE 估计（DL、REML、REML+HK）给出几乎相同的点估计 $\mathrm{RR} = 0.49$，CI 宽度从 DL 的 $[0.35, 0.70]$ 一路扩到 REML+HK 的 $[0.33, 0.73]$。比 FE-IV 的 $[0.60, 0.70]$ 显著更宽是合理的：RE 模型承认研究本来就不一样，CI 必须把这种差异纳入不确定性。

但合并 $\mathrm{RR} = 0.49$ 与 PI $[0.14, 1.76]$ 跨度极大的反差说明，异质性本身需要被正面量化。第 6 章就用 $Q$、$I^2$、$\tau^2$ 这一组指标系统刻画异质性，PI 不再重复定义，按本章定义直接引用。第 6 章末把所有方法（FE 与 RE）的合并结果汇入一张更新的对比表。