set.seed(2027)
# 模拟 8 项独立研究的数据
meta_data <- data.frame(
study = c("Wang 2019", "Li 2020", "Zhang 2018", "Chen 2021",
"Liu 2017", "Huang 2022", "Yang 2020", "Zhou 2019"),
ecosystem = c("grassland", "forest", "cropland", "grassland",
"wetland", "forest", "cropland", "grassland"),
n_t = c(15, 12, 20, 10, 18, 14, 16, 12),
n_c = c(15, 12, 20, 10, 18, 14, 16, 12),
mean_t = c(3.2, 2.9, 3.5, 2.8, 3.0, 3.3, 3.1, 2.7),
mean_c = c(3.8, 3.1, 3.9, 3.6, 3.4, 3.2, 3.7, 3.3),
sd_t = c(0.5, 0.4, 0.6, 0.3, 0.5, 0.4, 0.5, 0.3),
sd_c = c(0.6, 0.5, 0.7, 0.4, 0.6, 0.5, 0.6, 0.4)
)
meta_data13 荟萃分析入门
14 荟萃分析入门
在生态学研究中，同一个科学问题往往有许多独立研究给出了不同的答案。例如，“施肥是否会改变土壤微生物多样性？”这个问题，可能有几十篇论文分别在不同地点、不同土壤类型、不同施肥方案下进行了实验，得到了各不相同的结论——有的发现多样性显著下降，有的发现没有变化，还有的发现反而增加了。面对这些看似矛盾的结果，我们该如何得出一个可靠的总体结论？
这就是荟萃分析（Meta-analysis）要解决的问题。
14.1 什么是荟萃分析
荟萃分析是一种定量的文献综合方法，它通过统计手段将多个独立研究的结果合并在一起，计算出一个总体效应量（overall effect size），从而得到比任何单一研究都更可靠的结论。
与传统的文献综述（narrative review）不同，荟萃分析有三个核心特征：
- 定量化：不是简单地说”大多数研究发现了正效应”，而是精确计算效应的大小和置信区间
- 系统性：遵循严格的文献检索和筛选流程，尽量减少选择偏差
- 加权合并：每个研究的贡献不是等权的——样本量大、精度高的研究获得更大的权重
14.1.1 荟萃分析的典型应用场景
在生态学中，荟萃分析被广泛用于以下场景：
| 应用场景 | 示例问题 |
|---|---|
| 处理效应评估 | 氮沉降对土壤碳储量的影响有多大？ |
| 方法比较 | 免耕与传统耕作，哪种更有利于土壤微生物？ |
| 全球格局总结 | 全球范围内，增温对植物物候的影响是否一致？ |
| 争议问题澄清 | 生物多样性与生态系统功能之间是正相关还是无关？ |
14.1.2 荟萃分析的基本流程
一个完整的荟萃分析通常包括以下步骤：
提出研究问题 → 系统检索文献 → 筛选纳入研究 → 提取数据 → 计算效应量 → 合并分析 → 结果可视化本章将聚焦于后半段——从数据提取到结果可视化，帮助你掌握荟萃分析的核心技术。
14.2 效应量：荟萃分析的”通用货币”
不同研究可能使用不同的测量指标、不同的单位，甚至不同的实验设计。为了让这些研究的结果能够放在一起比较和合并，我们需要将它们转换为一个标准化的度量——这就是效应量（effect size）。
14.2.1 常用效应量类型
生态学荟萃分析中最常用的效应量包括：
| 效应量 | 公式 | 适用场景 |
|---|---|---|
| 对数响应比（lnRR） | \(\ln(\bar{X}_T / \bar{X}_C)\) | 处理组 vs 对照组的比值，生态学最常用 |
| 标准化均值差（Hedges’ g） | \(J \cdot ({\bar{X}_T - \bar{X}_C})/{S_p}\) | 两组均值差异，适合不同量纲的比较 |
| 相关系数（r） | Pearson 相关系数 | 两个连续变量的关联强度 |
其中，\(\bar{X}_T\) 和 \(\bar{X}_C\) 分别是处理组和对照组的均值，\(S_p\) 是合并标准差，\(J = 1 - \frac{3}{4(n_T + n_C - 2) - 1}\) 是 Hedges 小样本校正因子（当样本量较大时 \(J \approx 1\)）。
Tip为什么生态学偏爱对数响应比？
对数响应比（log response ratio, lnRR）在生态学中特别流行，原因有三：（1）它衡量的是相对变化（百分比变化），比绝对差异更有生态学意义；（2）对数变换使其近似正态分布，满足统计假设；（3）解释直观——lnRR = 0.2 意味着处理组比对照组高约 22%（\(e^{0.2} - 1 \approx 0.22\)）。
14.2.2 效应量的方差
每个效应量都有一个对应的方差（variance），它反映了该估计的不确定性。方差越小，说明这个研究的估计越精确，在合并分析中获得的权重也越大。
对于对数响应比，其方差的近似公式为：
\[ V_{lnRR} = \frac{SD_T^2}{n_T \cdot \bar{X}_T^2} + \frac{SD_C^2}{n_C \cdot \bar{X}_C^2} \]
其中 \(SD\) 是标准差，\(n\) 是样本量。
14.3 文献数据提取
在实际的荟萃分析中，数据提取是最耗时的环节之一。你需要从每篇纳入的论文中提取以下关键信息：
14.3.1 必须提取的数据
- 处理组和对照组的均值（mean）
- 标准差或标准误（SD 或 SE）
- 样本量（n）
- 研究的基本信息：作者、年份、研究地点、生态系统类型等
14.3.2 数据提取的注意事项
Warning数据提取中的常见陷阱
- 标准差 vs 标准误：论文中报告的可能是 SD，也可能是 SE。两者的关系是 \(SE = SD / \sqrt{n}\)，提取时务必区分清楚
- 从图中提取数据：当论文只提供图形而没有数值时，可以使用 WebPlotDigitizer 等工具从图中读取数据
- 多个处理水平：如果一篇论文包含多个施肥水平（低、中、高），每个水平与对照的比较可以作为独立的效应量，但需要注意非独立性问题
14.3.3 构建数据表
提取完成后，通常将数据整理为如下格式的表格：
| study | year | ecosystem | treatment | mean_t | sd_t | n_t | mean_c | sd_c | n_c |
|---|---|---|---|---|---|---|---|---|---|
| Wang | 2019 | grassland | N_addition | 3.2 | 0.5 | 15 | 3.8 | 0.6 | 15 |
| Li | 2020 | forest | N_addition | 2.9 | 0.4 | 12 | 3.1 | 0.5 | 12 |
14.4 R 实战：使用 metafor 包进行荟萃分析
metafor 是 R 中功能最全面的荟萃分析包，由 Wolfgang Viechtbauer 开发和维护。它支持多种效应量计算、固定效应和随机效应模型、异质性检验、森林图绘制等功能。
14.4.1 准备示例数据
我们构造一个模拟数据集：8 项研究分别考察了氮添加对土壤微生物 Shannon 多样性指数的影响。
这个数据集中，每行代表一项研究。mean_t 和 mean_c 分别是氮添加处理组和对照组的 Shannon 指数均值，sd_t 和 sd_c 是对应的标准差，n_t 和 n_c 是样本量。
14.4.2 计算效应量
使用 metafor::escalc() 函数计算对数响应比（lnRR）：
library(metafor)
# 计算对数响应比 (log response ratio)
meta_es <- escalc(
measure = "ROM", # ROM = Ratio of Means (即 lnRR)
m1i = mean_t, # 处理组均值
m2i = mean_c, # 对照组均值
sd1i = sd_t, # 处理组标准差
sd2i = sd_c, # 对照组标准差
n1i = n_t, # 处理组样本量
n2i = n_c, # 对照组样本量
data = meta_data,
slab = study # 用研究名称作为标签
)
# 查看计算结果
meta_es[, c("study", "yi", "vi")]escalc() 在原数据框中新增了两列：
yi：每个研究的效应量（lnRR）vi：每个效应量的方差
负值的 yi 表示氮添加降低了微生物多样性，正值表示增加了多样性。
14.4.3 拟合随机效应模型
在生态学荟萃分析中，我们几乎总是使用随机效应模型（random-effects model），因为不同研究之间的真实效应量很可能不完全相同——不同地点、不同生态系统、不同施肥量都会导致效应的异质性。
# 拟合随机效应模型（默认使用 REML 估计）
res <- rma(yi, vi, data = meta_es)
# 查看结果
summary(res)结果解读的关键信息：
- estimate（总体效应量）：所有研究合并后的 lnRR 估计值
- pval（p 值）：总体效应是否显著不为零
- ci.lb / ci.ub：95% 置信区间
- tau²：研究间方差，反映异质性大小
- I²：异质性占总变异的百分比。I² > 75% 通常认为异质性较高
- Q 检验：检验各研究的效应量是否同质（H₀: 所有研究的真实效应相同）
Note固定效应 vs 随机效应
固定效应模型假设所有研究估计的是同一个真实效应，研究间的差异仅来自抽样误差。随机效应模型则假设每个研究有自己的真实效应，这些效应围绕一个总体均值波动。在生态学中，由于研究条件的多样性，随机效应模型几乎总是更合理的选择。
14.4.4 将效应量转换为百分比变化
对数响应比的数值不太直观，我们通常将其转换为百分比变化来报告：
# 提取总体效应量和置信区间
overall <- coef(res)
# 转换为百分比变化
pct_change <- (exp(overall) - 1) * 100
pct_ci_lb <- (exp(res$ci.lb) - 1) * 100
pct_ci_ub <- (exp(res$ci.ub) - 1) * 100
cat(sprintf(
"氮添加使土壤微生物多样性平均变化了 %.1f%% (95%% CI: %.1f%% 到 %.1f%%)\n",
pct_change, pct_ci_lb, pct_ci_ub
))14.5 森林图：荟萃分析的标志性可视化
森林图（forest plot）是荟萃分析中最重要的图形，它将每个研究的效应量及其置信区间以图形方式展示出来，并在底部显示合并后的总体效应。
forest(
res,
header = TRUE, # 显示列标题
xlab = "对数响应比 (lnRR)",
mlab = "随机效应模型合并估计", # 合并菱形的标签
cex = 0.9, # 字体大小
refline = 0 # 参考线（无效应）
)14.5.1 如何阅读森林图
- 每一行代表一个研究，方块是该研究的效应量估计值，横线是 95% 置信区间
- 方块大小反映该研究的权重——方块越大，权重越高（通常意味着样本量更大或方差更小）
- 底部菱形是所有研究合并后的总体效应量，菱形的宽度是 95% 置信区间
- 竖直虚线（x = 0）是无效应参考线。如果某个研究的置信区间与这条线交叉，说明该研究的效应不显著
14.6 异质性分析与亚组比较
当 I² 较高时，说明研究之间存在较大的异质性。这时我们可以通过亚组分析（subgroup analysis）来探索异质性的来源。
例如，我们可以检验不同生态系统类型之间的效应是否存在差异：
# 按生态系统类型进行亚组分析（混合效应模型）
res_sub <- rma(yi, vi, mods = ~ ecosystem, data = meta_es)
summary(res_sub)这个模型将 ecosystem 作为调节变量（moderator），检验不同生态系统类型之间的效应量是否存在显著差异。如果 QM（调节变量的 Q 检验）显著，说明生态系统类型能够解释部分异质性。
14.7 发表偏倚检验
荟萃分析的一个重要假设是我们纳入了所有相关研究。但现实中，显著结果更容易被发表（“文件抽屉问题”），这可能导致发表偏倚（publication bias），使合并效应被高估。
14.7.1 漏斗图
漏斗图（funnel plot）是检测发表偏倚的直观工具：
funnel(res, xlab = "对数响应比 (lnRR)")在没有发表偏倚的情况下，散点应该围绕总体效应量对称分布，形成一个倒漏斗形。如果图形明显不对称，可能存在发表偏倚。
14.7.2 Egger 回归检验
除了目视检查漏斗图，还可以用 Egger 检验进行统计检验：
regtest(res)如果 p 值显著（通常 < 0.05），则提示可能存在发表偏倚。但需要注意，当纳入研究数量较少（< 10）时，这些检验的统计效力较低。
14.8 完整分析流程总结
14.9 课后练习
14.9.1 练习 1：理解效应量
以下是两项研究的数据：
| 研究 | 处理组均值 | 处理组SD | 处理组n | 对照组均值 | 对照组SD | 对照组n |
|---|---|---|---|---|---|---|
| 研究 A | 25.3 | 4.2 | 20 | 30.1 | 5.0 | 20 |
| 研究 B | 18.7 | 3.1 | 8 | 20.5 | 3.8 | 8 |
（a）手动计算两项研究的对数响应比 lnRR（提示：\(lnRR = \ln(\bar{X}_T / \bar{X}_C)\)）
（b）用 escalc() 函数验证你的手动计算结果
（c）哪项研究的效应量方差更大？为什么？
14.9.2 练习 2：完整荟萃分析
下面的数据来自 6 项研究，考察了有机肥施用对土壤微生物生物量碳（MBC）的影响：
exercise_data <- data.frame(
study = c("赵 2018", "钱 2019", "孙 2020", "李 2021", "周 2017", "吴 2022"),
mean_t = c(285, 310, 245, 330, 270, 295),
sd_t = c(42, 55, 38, 60, 45, 50),
n_t = c(10, 15, 12, 8, 20, 14),
mean_c = c(220, 250, 210, 240, 230, 235),
sd_c = c(35, 48, 30, 45, 40, 42),
n_c = c(10, 15, 12, 8, 20, 14)
)请完成以下任务：
（a）使用 escalc() 计算每项研究的对数响应比
（b）使用 rma() 拟合随机效应模型，报告总体效应量及其 95% 置信区间
（c）将总体效应量转换为百分比变化，并用一句话描述结果
（d）绘制森林图和漏斗图
14.9.3 练习 3：思考题
（a）为什么荟萃分析通常使用随机效应模型而非固定效应模型？试从生态学研究的特点出发解释。
（b）一位同学在做荟萃分析时，只在 Google Scholar 上用英文关键词搜索了文献。这种做法可能带来什么问题？如何改进？
（c）如果漏斗图显示明显的不对称性，是否一定意味着存在发表偏倚？还有哪些因素可能导致漏斗图不对称？