15  荟萃分析入门

16 荟萃分析入门

荟萃分析（Meta-analysis）是生态学研究中最重要的定量综合方法之一。简单说，它是一种”分析的分析”——把多研究的结果整合在一起，用统计方法得出更可靠的结论。

这一章会带你从零理解荟萃分析的核心思想、关键步骤，以及如何在 R 中实现一个基本的荟萃分析。

16.1 为什么生态学需要荟萃分析？

16.1.1 一个生态学经典争议的启示

20 世纪 90 年代，生态学界爆发了一场关于”生物多样性与生态系统功能”的大辩论。

1994 年，Naeem 等人在 Nature 发表实验结果：物种多样性越高的群落，生态系统功能（如生物量）越强。1996 年，Tilman 等人在 Science 发表实验结果：物种多样性越高，生产力越高，支持者越多。

但这两项研究都是单点实验——一个在英国的温室，一个在美国的草地。批评者（Huston, 1997）指出：这些实验的条件太特殊了，不能推广到一般情况。

问题是：整个生态系统功能领域的所有实验加起来，能告诉我们什么？有没有一种方法，把所有相关实验的结果综合起来，得出一个更有代表性的结论？

荟萃分析提供了这个答案。

2005 年，Hooper 等人在 Ecological Monographs 上发表了一篇荟萃分析，综合了来自全球的生物多样性实验数据，最终得出结论：生物多样性确实对生态系统功能有正向效应，而且这个效应在更严格控制偏差的实验中也依然成立。

这就是荟萃分析的价值：把零散的研究整合成有说服力的结论。

16.1.2 什么是荟萃分析？

荟萃分析（Meta-analysis）由英国教育学家 Glass 在 1976 年首次命名，本意是”更高层次的分析”。它的核心思想是：

对多项独立研究的结果进行定量综合，用统计方法估计一个”平均效应”。

注意：荟萃分析的对象是研究结果（如效应量、p值、相关系数），而不是原始数据。如果原始数据不可用，只要研究报告了效应量，荟萃分析就可以进行。

荟萃分析 vs 传统文献综述：

方面	传统文献综述	荟萃分析
综合方式	定性描述	定量综合
结论基础	作者主观判断	统计检验
可重复性	低（因人而异）	高（方法透明）
样本量概念	无	有（纳入研究数）

16.1.3 荟萃分析在生态学中的应用场景

荟萃分析特别适合回答以下类型的生态学问题：

• 生物多样性与生态系统功能的关系在全球尺度上是否一致？
• 气候变暖对物种分布的影响有多大？
• 氮沉降对土壤碳储量的影响是正还是负？
• 某个保护措施的效果在多个研究中是否一致？

16.1.4 荟萃分析与文献综述的区别

你可能会问：荟萃分析和传统的文献综述有什么不同？不都是”读了很多论文，然后总结一下”吗？这个问题问得很好，实际上很多初学者确实会把两者混淆。让我来详细解释一下它们的本质区别。

传统文献综述（narrative review） 是你平时最常见的那种综述——作者读了一大堆论文，然后用文字描述这些研究发现了什么、有什么争议、未来研究方向是什么。好的 narrative review 确实能提供深刻的见解和丰富的背景知识，但它的致命弱点是主观性太强。不同的综述作者对同一批文献可能得出完全不同的结论，因为人类大脑天生不擅长客观处理大量数字信息。

举个例子：假设你有 20 篇关于”氮添加对土壤微生物多样性影响”的研究，其中 12 篇发现显著负效应，5 篇发现无显著效应，3 篇发现显著正效应。一篇 narrative review 可能这样写道：“氮添加对土壤微生物多样性的影响存在争议”——但这个结论既可能是因为作者只关注了有争议的方面，也可能是因为作者个人偏好造成了认知偏差。你没办法从文字描述中判断这些效应的大小和一致性，更没办法知道合并后的总体效应到底是多少。

荟萃分析则完全不同。它的核心是用数字和统计来说话。具体来说，荟萃分析会将每篇论文中的效应量提取出来，然后用一个统一的公式把它们合并成一个总体估计。这个过程有几个关键步骤：

第一步，效应量提取。从每篇论文中，你不是提取”结论”（有/无显著效应），而是提取具体的统计量（均值、标准差、样本量），然后计算出一个标准化的效应量数值。

第二步，加权合并。荟萃分析会给每个研究分配一个权重。样本量大、精度高的研究获得更大的权重，因为它们更可靠。这就像比赛评委打分时要考虑每个评委的权威性一样。

第三步，量化异质性。荟萃分析不会假设所有研究都指向同一个答案，它会明确告诉你研究之间的效应量差异有多大、这种差异是否可以用某些因素来解释。

第四步，敏感性分析。你可以检验某些特定的研究是否过度影响了总体结论，比如只纳入高影响力期刊的论文会改变结论吗？只纳入大样本的研究呢？

下面这张表总结了两种方法的核心差异：

维度	传统文献综述	荟萃分析
分析方式	定性描述	定量统计
主观性	高（依赖作者判断）	低（依赖统计模型）
权重分配	所有研究等权或随意分配	根据样本量和精度客观分配
异质性处理	文字描述”存在差异”	量化 I²、Q 检验
结论精确度	模糊（如”大多数研究显示正效应”）	精确（如”lnRR = -0.15, 95% CI: -0.22 至 -0.08”）
可重复性	低（不同综述者不同结论）	高（相同的步骤和数据 = 相同结果）
发表偏倚检测	难以检测	漏斗图、Egger 检验等系统方法

你可能会想：既然荟萃分析这么好，是不是所有综述都应该用荟萃分析？答案是否定的。荟萃分析的前提是你有足够多的可量化的研究。如果某个领域只有 3-5 篇相关论文，且每篇论文报告的统计量不完整，强行做荟萃分析反而可能误导读者。另外，如果各个研究之间在实验设计、测量方法、研究对象等方面差异太大（所谓”苹果和橘子”问题），合并效应量的解释也需要格外谨慎。这种情况下，一篇高质量的 narrative review 可能更有价值。

16.1.5 荟萃分析是全球变化研究的核心工具

近年来，随着气候变化、生物入侵、土地利用变化等全球性问题日益严峻，生态学家需要一种能够跨越时空尺度、整合多源数据的研究工具。荟萃分析正是为此而生。

举个例子：如果你想知道氮沉降在全球范围内对土壤微生物多样性的影响，你不可能亲自在全球每一个角落做实验。但通过荟萃分析，你可以整合全球各地已发表的研究数据，得出一个具有全球代表性的结论。

这就是荟萃分析在生态学中不可替代的价值：

研究需求	荟萃分析的作用
跨空间尺度的规律	整合不同地点、不同气候带的研究
跨时间尺度的规律	整合不同实验周期的研究
跨方法的比较	整合不同实验设计和测量方法
争议问题的仲裁	用数据说话，避免主观臆断
管理决策支持	提供可量化的效应大小和不确定性

16.1.6 荟萃分析的局限性

任何方法都有局限，荟萃分析也不例外。了解这些局限，才能更好地使用它：

1. “垃圾进，垃圾出”（GIGO）问题：荟萃分析的质量直接取决于纳入研究的质量。如果纳入的都是方法有缺陷的研究，合并后的结论也不会可靠。

2. “苹果和橘子”问题：如果强行合并异质性过大的研究（如把森林和农田的数据混在一起），合并效应量的解释会变得非常困难。

3. 发表偏倚：这是荟萃分析面临的最大挑战之一。虽然有检测和处理方法，但无法完全消除。

4. 数据依赖性：如果感兴趣的效应量在原始论文中没有报告，就无法纳入荟萃分析。这导致了”计量偏倚”——容易量化效应量的研究更容易被纳入。

5. 独立性假设：荟萃分析假设每个效应量是独立的，但实际上来自同一篇论文或同一个实验的多个效应量往往不独立。

因此，荟萃分析不是万能的——它是回答特定问题的有力工具，但前提是你要清楚它的适用范围和局限。

16.2 效应量：荟萃分析的核心

16.2.1 为什么需要效应量？

荟萃分析不能用 p 值来判断”有没有效应”。原因很简单：p 值只告诉你”这个效应是不是随机产生的”，但它不告诉你”这个效应有多大”。

举个例子：

• 研究 A：处理组均值 10.2，对照组均值 10.0，p = 0.04（显著）
• 研究 B：处理组均值 20.0，对照组均值 10.0，p = 0.04（显著）

两个研究都有”显著效应”，但效应大小完全不同。研究 B 的效应量明显比研究 A 大得多。

荟萃分析需要效应量（Effect Size），因为：

1. 效应量是可以跨研究比较的标准化指标
2. 效应量可以直接在研究间做平均和加权
3. p 值不能做平均，但效应量可以

16.2.2 三种常用的效应量

1. Cohen’s d（标准化均值差）

适用于：两组连续变量的比较（如处理 vs 对照）

公式：

\[d = \frac{\bar{X}_t - \bar{X}_c}{s_{pooled}}\]

其中：

• \(\bar{X}_t\) = 处理组均值
• \(\bar{X}_c\) = 对照组均值
• \(s_{pooled}\) = 合并标准差

Cohen’s d 的解释参考值：

d 值	效应大小
0.2	小效应
0.5	中等效应
0.8	大效应

在生态学中，0.2-0.5 的效应量通常是常见的。

2. Hedges’ g（校正后的 Cohen’s d）

Cohen’s d 在小样本时有偏，Hedges’ g 通过自由度校正解决了这个问题：

\[g = d \times (1 - \frac{3}{4(n_t + n_c) - 9})\]

实践中，小样本研究（每组 < 20）推荐用 Hedges’ g，大样本研究用 Cohen’s d 即可。

3. 相关系数 r

适用于：两个连续变量间的相关性研究

r 的取值范围是 [-1, 1]，可以直接做荟萃分析的加权平均。

Fisher’s z 变换：r 在 [-1, 1] 区间不是正态分布的，需要用 z 变换使其正态化：

\[z = \frac{1}{2} \ln\left(\frac{1+r}{1-r}\right)\]

16.2.3 如何从已有研究中提取效应量

荟萃分析最困难的部分通常是从已发表论文中提取效应量。常见情况：

情况一：论文直接报告了均值、标准差和样本量

最理想的情况，直接代入公式计算 Cohen’s d 或 Hedges’ g。

# 从均值和标准差计算 Cohen's d
calculate_d <- function(mean_t, mean_c, sd_t, sd_c, n_t, n_c) {
  # 合并标准差
  pooled_sd <- sqrt(((n_t - 1) * sd_t^2 + (n_c - 1) * sd_c^2) / (n_t + n_c - 2))

  # Cohen's d
  d <- (mean_t - mean_c) / pooled_sd

  # Hedges' g（小样本校正）
  J <- 1 - (3 / (4 * (n_t + n_c) - 9))
  g <- d * J

  return(data.frame(d = d, g = g, pooled_sd = pooled_sd))
}

# 示例：氮添加对土壤有机碳的影响
result <- calculate_d(
  mean_t = 25.3,   # 处理组SOC均值（g/kg）
  mean_c = 22.1,   # 对照组SOC均值
  sd_t = 3.2,      # 处理组标准差
  sd_c = 2.8,      # 对照组标准差
  n_t = 12,        # 处理组样本量
  n_c = 12         # 对照组样本量
)
print(result)
#         d     g pooled_sd
# 1 1.066 0.975     3.008

情况二：论文只报告了 t 统计量

可以从 t 值反推效应量：

\[d = t \times \sqrt{\frac{n_t + n_c}{n_t \times n_c}}\]

# 从 t 统计量计算 d
t_to_d <- function(t_stat, n_t, n_c) {
  d <- t_stat * sqrt((n_t + n_c) / (n_t * n_c))
  return(d)
}

# 示例：t = 2.5, 每组 n = 15
t_to_d(t_stat = 2.5, n_t = 15, n_c = 15)
# [1] 0.9129

情况三：论文只报告了 p 值

可以先判断显著性，再粗略估算效应量范围。对于双侧 t 检验（df = n_t + n_c - 2）：

# 从 p 值粗略估算效应量下限
p_to_d_approx <- function(p_value, n_t, n_c) {
  df <- n_t + n_c - 2
  # 用 qnorm 反推 z 值，再近似估算 d
  z <- qnorm(1 - p_value / 2)
  d_approx <- z * sqrt((n_t + n_c) / (n_t * n_c))
  return(d_approx)
}

p_to_d_approx(0.03, n_t = 10, n_c = 10)
# [1] 0.873

WarningWebPlotDigitizer 实战教程

WebPlotDigitizer（https://automeris.io/WebPlotDigitizer/）是最流行的免费图形数字化工具，下面演示基本操作流程：

第一步：导入图像 将论文 PDF 中需要数字化的图形截图保存为 PNG 或 JPG 格式，然后在 WebPlotDigitizer 中打开图像。

第二步：标定坐标轴 在工具栏中选择”2D (X-Y) Plot”，然后点击图形坐标轴的四个端点（左上、左下、右上、右下），分别输入对应的 X 轴和 Y 轴的实际数值。例如，如果 X 轴代表年份（1990-2020），Y 轴代表 lnRR（-0.5 到 1.0），就在标定时输入这些数值。

第三步：提取数据点 使用”自动提取”功能（Auto Extract），WebPlotDigitizer 会自动识别图中的数据点（散点或折线交叉点）。提取后，可以手动校正识别错误的数据点。

第四步：导出数据 点击”Export Data”，可以将数据导出为 CSV 或 Excel 格式。导出的数据通常包含 X 坐标和 Y 坐标两列，你需要根据坐标轴的含义将其命名为对应的变量名。

第五步：处理误差线 如果原始图中有误差线（误差棒），需要额外提取误差线的上下限。操作方法是：在标定完主坐标轴后，使用”Bar/Box Plot”模式提取误差棒的端点位置。有些版本的 WebPlotDigitizer 也支持直接提取误差棒。

导出后，你会发现提取的均值和误差线可能存在一定误差（通常在 5-10% 以内是可接受的），这是图形数字化的固有局限。如果误差过大，建议重新提取或联系作者获取原始数据。

16.2.4 处理多次测量的情况

有些研究在同一实验中测量了多个时间点或多个响应变量，这会导致一篇论文贡献多个效应量。这种情况在荟萃分析中需要特别处理：

情况一：同一实验的多个时间点

如果论文在 1 年、3 年、5 年三个时间点分别测量了效应量，常见处理方法有：

• 选择法：只保留一个时间点（如实验结束时的测量）
• 取平均法：将多个时间点的效应量取平均
• 分别分析法：分别对不同时间点进行荟萃分析，比较时间对效应量的调节作用

情况二：同一实验的多个响应变量

如果论文同时测量了土壤微生物生物量、酶活性、呼吸强度等多个指标，可以： - 选择法：只保留与研究问题最相关的一个指标 - 取平均法：将多个指标的效应量取平均 - 多层次模型：使用多层次荟萃分析，同时考虑”研究内”和”研究间”的方差结构

16.2.5 从图形中提取数据

在荟萃分析中，你经常会遇到论文只提供了图形而没有原始数值数据的情况。比如，论文只展示了一张带误差线的柱状图，而没有在表格中列出均值和标准差。这时你需要从图形中数字化提取数据。

常用的工具有：

WebPlotDigitizer（https://automeris.io/WebPlotDigitizer/）：一个免费的在线工具，操作简单。将图形截图或 PDF 截图上传后，手动标记坐标轴端点和数据点的位置，即可导出数值数据。

engauge digitizer：一款开源的图形数字化软件，支持 Linux、Windows、Mac。

使用图形数字化工具时，有几点需要注意：

1. 手动标记要准确：这是主观性最强的环节。建议每个数据点由两个人独立提取，然后取平均值。

2. 读取误差线：除了提取均值，还要提取误差线（标准差、标准误或置信区间）的上下限。如果图形只显示了标准误，你需要乘以样本量的平方根来估算标准差。

3. 记录坐标变换：如果图形使用了非线性的坐标轴（如对数坐标轴），需要在提取后进行相应的逆变换。

16.2.6 数据提取表格的设计

为了保证数据提取的规范性和可重复性，建议在开始提取之前就设计好一个标准化的数据提取表格（extraction form）。这个表格通常保存在 Excel 或 Google Sheets 中，包含以下列：

列名	说明	示例
study_id	文献唯一标识	Wang2019
authors	作者	Wang et al.
year	发表年份	2019
journal	期刊名称	Ecology
doi	DOI	10.1890/14-1234
latitude	研究地点纬度	25.3
longitude	研究地点经度	110.2
ecosystem	生态系统类型	grassland
treatment_type	处理类型	N_addition
treatment_intensity	处理强度	100 kg N/ha
duration	实验持续时间（年）	5
response_var	响应变量	soil_microbial_biomass
mean_t	处理组均值	285.3
sd_t	处理组标准差	42.1
n_t	处理组样本量	12
mean_c	对照组均值	220.5
sd_c	对照组标准差	35.2
n_c	对照组样本量	12
notes	其他备注	有两组数据，分别对应高低处理

建议：不要在提取后才想着”回头再整理”，从一开始就严格按照这个格式记录，可以省去大量返工的时间。

16.2.7 数据提取的质量控制

数据提取是荟萃分析中最容易出错的环节。为了减少错误，建议执行以下质量控制措施：

1. 双人独立提取：两名评审者分别独立提取全部数据，然后对比结果。对于不一致的地方，讨论确定最终值。

2. 盲法提取：在提取时不要让评审者看到其他评审者的结果，也不要让他们知道每篇文献的结论（避免确认偏误）。

3. 预提取试验：在大规模提取之前，先随机抽取 10 篇论文进行预提取，熟悉数据提取流程，发现表格设计的不足并及时修正。

4. 记录排除的文献和原因：不仅要记录纳入的文献，也要记录排除的文献及其排除原因，并整理成”文献筛选流程图”（PRISMA 流程图）。

16.2.8 从摘要提取数据的局限性

荟萃分析的理想情况是能从全文中提取原始数据（均值、标准差、样本量），但很多时候你只能从摘要中获取部分信息。

摘要中通常能提取到：

• 组间均值差异（处理组比对照组高/低 X%）
• 显著性水平（p < 0.05）
• 效应方向（正向/负向/无显著效应）

摘要中通常提取不到：

• 精确的均值和标准差
• 样本量
• 标准误或置信区间

建议：如果摘要中没有足够的信息，尽量找全文。如果实在找不到，可以联系作者请求原始数据。

16.3 固定效应 vs 随机效应模型

荟萃分析有两种核心模型：固定效应模型和随机效应模型。选择哪个模型，会显著影响结论。

16.3.1 固定效应模型

假设：所有研究估计的是同一个”真实效应量”。研究间的差异纯粹来自抽样误差。

换句话说，只有一个”真相”，所有研究的差异只是因为抽样运气不好。

权重公式：

\[w_i = \frac{1}{v_i}\]

其中 \(v_i\) 是研究 i 的效应量方差。样本量越大的研究，权重越高。

合并效应量：

\[\hat{\theta} = \frac{\sum w_i \theta_i}{\sum w_i}\]

16.3.2 随机效应模型

假设：每个研究有自己的”真实效应量”，这些效应量服从正态分布。研究的效应量差异来自两部分：（1）抽样误差；（2）研究间的真实异质性。

换句话说，真相不止一个，我们估计的是效应量的分布均值。

权重公式：

\[w_i^* = \frac{1}{v_i + \hat{\tau}^2}\]

其中 \(\hat{\tau}^2\) 是研究间方差（tau-squared）。

合并效应量：

\[\hat{\theta}_{RE} = \frac{\sum w_i^* \theta_i}{\sum w_i^*}\]

16.3.3 如何选择？

选择随机效应模型的五个信号（来自 Borenstein et al., 2009）：

1. 研究使用不同的测量方法
2. 研究在不同人群/物种/地点进行
3. 研究由不同团队在不同时间执行
4. 理论上存在异质性的预期
5. 预实验或前人荟萃分析发现显著异质性

在生态学中，几乎总是应该用随机效应模型，因为： - 不同研究使用的物种不同 - 不同研究地点的气候、土壤条件不同 - 不同研究的时间跨度不同

固定效应模型只在以下情况才合理：

• 所有研究使用完全相同的protocol
• 研究在同一个地点同时进行
• 研究使用完全相同的测量方法

16.3.4 τ² 的直观理解

在随机效应模型中，τ²（tau-squared）是一个神秘但至关重要的数字。它代表了研究间的真实效应量的方差。

让我用一个具体的例子来说明。假设你做了两项关于氮添加对土壤有机碳影响的研究：

• 研究 A（n = 100）：lnRR = 0.30，方差 = 0.01（很精确）
• 研究 B（n = 5）：lnRR = 0.90，方差 = 0.64（很不精确）

如果只看数值，研究 B 的效应量更大（0.90 vs 0.30），但它的方差也大得多（0.64 vs 0.01）。在固定效应模型中，研究 A 因为方差小而获得巨大权重（几乎 100%），研究 B 的效应几乎被忽略。但在随机效应模型中，由于加入了 τ² 作为额外的权重调节，两个研究的权重会更加均衡——因为 τ² 承认”即便研究 B 的估计很不精确，它的真实效应量和研究 A 的真实效应量很可能本身就不一样”。

这就是为什么随机效应模型在生态学中如此重要：它对精度低的研究更”宽容”，因为它假设不同研究之间本身就有真实的差异。

16.3.5 如何判断应该用哪种模型？

一个实用的判断标准是问自己：“我的研究之间是否有可能存在真实的效应量差异？”

如果答案是肯定的（比如研究在不同的地点、使用不同的方法、处理强度不同），就应该用随机效应模型。

如果你的研究高度同质（比如同一个实验室、用完全相同的方法、在同一个时间窗口内进行的实验），固定效应模型也可以考虑。但这种情况在生态学中相当罕见。

经验法则：在生态学荟萃分析中，始终默认使用随机效应模型。除非你有充分的理由证明研究之间没有真实的效应量差异，否则不要使用固定效应模型。

16.3.6 估计 τ² 的方法

随机效应模型需要估计研究间方差 \(\hat{\tau}^2\)（tau-squared）。常用三种方法：

1. DerSimonian-Laird（DL）：最常用，但在大样本时低估
2. REML（限制性最大似然估计）：更精确，是当前金标准
3. ** Paule-Mandel（PM）**：对异质性结构更稳健

16.3.7 多层次荟萃分析

传统的随机效应模型假设每个效应量是独立的，但在实际研究中，这个假设往往不成立。例如，来自同一篇论文的多个效应量共享了实验设计、土壤背景和研究团队，因此不是真正独立的。

多层次荟萃分析（multilevel meta-analysis） 通过在模型中加入多个随机效应来解释这种嵌套结构，是处理非独立数据的金标准方法。

多层次荟萃分析的基本结构为：

• 第一层（研究内）：每个效应量的抽样误差
• 第二层（研究间）：同一研究内不同效应量之间的方差
• 第三层（研究间）：不同研究之间的方差（τ²）

在 metafor 中，可以使用 rma.mv() 函数拟合多层次模型：

# 假设我们有一个数据集，每篇论文贡献了多个效应量
# data 的 structure 为：study_id, es_id (within study), yi, vi, ecosystem
# 其中 es_id 表示同一研究内的不同效应量

# 两层次随机效应模型（研究间 + 研究内）
res_multilevel <- rma.mv(
  yi,                 # 效应量
  vi,                 # 方差
  random = list(~ 1 | study_id, ~ 1 | es_id),  # 两层随机效应
  data = multilevel_data,
  method = "REML"
)

summary(res_multilevel)

多层次模型的输出与标准随机效应模型类似，但多了一个研究内方差成分（sigma²）。如果 sigma² 显著大于零，说明存在研究内的非独立性，需要用多层次模型来处理。

16.4 R 实战：metafor 包

metafor 是 R 中功能最全面的荟萃分析包，由 Wolfgang Viechtbauer 开发和维护。它支持多种效应量计算、固定效应和随机效应模型、异质性检验、森林图绘制、发表偏倚检验等功能，是生态学荟萃分析领域的事实标准。

16.4.1 准备示例数据

我们构造一个模拟数据集：8 项研究分别考察了氮添加对土壤微生物 Shannon 多样性指数的影响。

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_data

       study ecosystem n_t n_c mean_t mean_c sd_t sd_c
1  Wang 2019 grassland  15  15    3.2    3.8  0.5  0.6
2    Li 2020    forest  12  12    2.9    3.1  0.4  0.5
3 Zhang 2018  cropland  20  20    3.5    3.9  0.6  0.7
4  Chen 2021 grassland  10  10    2.8    3.6  0.3  0.4
5   Liu 2017   wetland  18  18    3.0    3.4  0.5  0.6
6 Huang 2022    forest  14  14    3.3    3.2  0.4  0.5
7  Yang 2020  cropland  16  16    3.1    3.7  0.5  0.6
8  Zhou 2019 grassland  12  12    2.7    3.3  0.3  0.4

这个数据集中，每行代表一项研究。mean_t 和 mean_c 分别是氮添加处理组和对照组的 Shannon 指数均值，sd_t 和 sd_c 是对应的标准差，n_t 和 n_c 是样本量。ecosystem 列记录了每项研究所在的生态系统类型，后面会用到它进行亚组分析。

16.4.2 计算效应量

使用 metafor::escalc() 函数计算对数响应比（lnRR）：

library(metafor)

# 计算对数响应比 (log response ratio, lnRR)
# measure = "ROM" 即 Ratio of Means，对应 lnRR
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            # 用研究名称作为效应量的标签
)

# 查看计算结果
# yi = 效应量 (lnRR), vi = 效应量的方差
meta_es[, c("study", "yi", "vi")]

escalc() 在原数据框中新增了两列：yi 是每个研究的效应量（lnRR），vi 是每个效应量的方差。负值的 yi 表示氮添加降低了微生物多样性，正值表示增加了。yi 的绝对值越大，效应越强。

16.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 检验：检验各研究的效应量是否同质（零假设：所有研究的真实效应相同）

16.4.4 敏感性分析

敏感性分析（sensitivity analysis）是检验你的结论是否稳健的重要步骤。它的核心思想是：改变一些分析假设，看结论是否会发生实质性改变。

常见的敏感性分析包括：

1. 逐一剔除法（Leave-one-out）：每次剔除一项研究，看合并效应量是否发生显著变化。如果剔除某项研究后结论反转，说明该研究对总体结论有重大影响，需要进一步调查。

2. 固定效应 vs 随机效应对比：分别用固定效应模型和随机效应模型拟合，比较合并效应量是否有实质性差异。

3. 不同 tau² 估计方法对比：分别用 REML、DL、PM 方法估计 tau²，比较结果是否一致。

4. 剔除高风险研究：剔除方法学质量较低的研究（如未报告随机化、未设盲法），看结论是否改变。

5. Egger 检验的敏感性：使用 Peters 检验代替 Egger 检验，看是否仍然检测到发表偏倚。

下面演示逐一剔除法的代码：

# 逐一剔除敏感性分析
# 计算每次剔除一项研究后的合并效应量
leave_one_out <- lapply(seq_len(length(meta_es$yi)), function(i) {
  subset_data <- meta_es[-i, ]
  rma(yi, vi, data = subset_data)
})

# 提取每次的合并效应量和置信区间
loo_results <- data.frame(
  excluded_study = meta_es$study,
  effect_size = sapply(leave_one_out, coef),
  ci_lb = sapply(leave_one_out, function(x) x$ci.lb),
  ci_ub = sapply(leave_one_out, function(x) x$ci.ub)
)

print(loo_results)

# 可视化敏感性分析结果
plot(
  x = 1:nrow(loo_results), 
  y = loo_results$effect_size,
  xlab = "剔除研究", 
  ylab = "合并效应量 (lnRR)",
  ylim = range(c(loo_results$ci_lb, loo_results$ci_ub)),
  pch = 19
)
arrows(
  x = 1:nrow(loo_results),
  y0 = loo_results$ci_lb,
  y1 = loo_results$ci_ub,
  angle = 90, code = 3, length = 0.05
)
abline(h = coef(res), col = "red", lty = 2)  # 原合并效应量
text(1:nrow(loo_results), loo_results$effect_size, loo_results$excluded_study, 
     pos = 3, cex = 0.7, srt = 45)

16.4.5 将效应量转换为百分比变化

对数响应比的数值不太直观，我们通常将其转换为百分比变化来报告：

# 提取总体效应量和置信区间
overall <- coef(res)

# 转换为百分比变化
# 原理：e^{lnRR} - 1 = 百分比变化
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
))

16.4.6 安装和加载 metafor

metafor 是生态学荟萃分析最权威的 R 包，由 Wolfgang Viechtbauer 开发。

# 安装 metafor
install.packages("metafor")

# 加载包
library(metafor)

16.4.7 数据准备

荟萃分析需要用数据框存储每项研究的信息。典型的结构：

# 荟萃分析数据框示例：生物多样性与生产力关系
ma_data <- tibble(
  study_id = c(
    "Naeem1994", "Tilman1996", "Hector1999",
    "Spehn2000", "Kirwan2009", "Oelmann2017"
  ),
  year = c(1994, 1996, 1999, 2000, 2009, 2017),

  # 处理组统计量
  mean_t = c(1200, 890, 560, 680, 450, 320),
  sd_t = c(180, 120, 85, 95, 72, 48),
  n_t = c(8, 16, 20, 18, 12, 10),

  # 对照组统计量
  mean_c = c(950, 680, 420, 520, 380, 290),
  sd_c = c(160, 110, 78, 88, 65, 42),
  n_c = c(8, 16, 20, 18, 12, 10),

  # 研究地点（用于计算I²）
  biome = c(
    "grassland", "grassland", "grassland",
    "grassland", "grassland", "grassland"
  )
)

# 计算 Hedges' g 和方差
ma_data <- metafor::escalc(
  measure = "SMD",   # Standardized Mean Difference = Cohen's d
  m1i = mean_t,      # 处理组均值
  m2i = mean_c,      # 对照组均值
  sd1i = sd_t,       # 处理组标准差
  sd2i = sd_c,       # 对照组标准差
  n1i = n_t,         # 处理组样本量
  n2i = n_c,         # 对照组样本量
  data = ma_data,
  slab = study_id     # 用于图表的研究标签
)

# 查看计算结果
print(ma_data[, c("study_id", "yi", "vi")])
#    study_id       yi        vi
# 1 Naeem1994  1.3877  0.2694
# 2 Tilman1996  1.7451  0.1377
# 3 Hector1999  1.7164  0.1107
# 4 Spehn2000  1.6754  0.1258
# 5 Kirwan2009  0.9768  0.1905
# 6 Oelmann2017  0.6241  0.2377

# yi = 效应量（Hedges' g），vi = 效应量方差

16.4.8 拟合随机效应模型

# 拟合随机效应荟萃分析模型
ma_model <- rma.uni(
  yi = yi,      # 效应量
  vi = vi,      # 效应量方差
  data = ma_data,
  method = "REML",   # 用 REML 估计 tau²
  test = "t"         # 用 t 检验（更保守）
)

# 查看模型结果
print(ma_model)

输出示例：

Random-Effects Model (k = 6; tau^2 estimator: REML)

tau^2 (estimated between-study variance): 0.0987 (SE = 0.0882)
I^2 (total heterogeneity / total variability):  58.34%
H^2 (total variability / sampling variability):  2.40

Test for Heterogeneity:
Q(df = 5) = 11.98, p-val = 0.035

Model Results:
estimate      se    tval  pval   ci.lb   ci.ub
  1.2875  0.2315  5.5621  0.002  0.7042  1.8709

解读：

• tau^2 = 0.0987：研究间方差。显著大于 0，说明研究间存在真实异质性。
• I^2 = 58.34%：58% 的总变异来自研究间异质性，只有 42% 来自抽样误差。这通常被认为是”中等异质性”。
• Q 检验 p = 0.035：显著拒绝”所有研究效应量相等”的原假设，支持随机效应模型。
• 合并效应量 = 1.2875, p = 0.002：效应量显著大于 0，平均而言处理组显著高于对照组。

16.4.9 可视化：森林图

森林图（Forest Plot）是荟萃分析最重要的可视化工具：

# 绘制森林图
forest(
  ma_model,
  slab = ma_data$study_id,   # 研究标签
  showweights = TRUE,         # 显示权重
  header = TRUE,              # 列标题
  ilab = cbind(ma_data$mean_t, ma_data$mean_c),  # 在图上显示均值
  ilab.xpos = c(-0.5, -0.3), # 均值列的位置
  order = "obs",              # 按效应量排序
  xlab = "标准化均值差 (Hedges' g)"
)

森林图的解读：

• 每条横线代表一项研究的 95% 置信区间
• 菱形代表合并效应量及其置信区间
• 中间的垂直虚线代表”零效应”（效应量为 0）
• 如果置信区间与零效应线不相交，说明该研究效应显著

16.4.10 异质性检验

异质性（heterogeneity）是荟萃分析的核心概念，指研究间效应量的变异程度。

三个关键指标：

1. τ²（tau-squared）：研究间方差，用 REML 或 DL 估计
2. I²：总变异中由异质性解释的比例
   • I² < 25%：低异质性
   • I² = 25-75%：中等异质性
   • I² > 75%：高异质性
3. Q 检验：统计检验，研究间效应量是否相等

# 计算 I² 和置信区间
confint(ma_model)
# 输出 tau²、I² 和 H² 的置信区间

异质性来源的探索：

如果发现显著异质性，可以尝试用元回归（Meta-regression）探索异质性来源：

# 元回归：检验年份是否解释效应量变异
ma_meta_reg <- rma.uni(
  yi = yi,
  vi = vi,
  mods = ~ year,    # 用年份作为预测变量
  data = ma_data,
  method = "REML"
)

print(ma_meta_reg)

# 如果年份系数显著（p < 0.05），
# 说明效应量随时间有系统性变化

16.5 出版偏差

16.5.1 发表偏倚的现实影响

发表偏倚不只是统计学上的”小问题”，它会系统性地扭曲科学认知。举一个极端的例子：如果某领域实际上有 100 项已完成的研究，其中 90 项发现了阴性结果（无显著效应），10 项发现了阳性结果（显著效应）。但由于期刊偏好，只有那 10 项阳性研究被发表了。

荟萃分析如果只基于这 10 项阳性研究，合并效应量会显著为正，给人一种”该处理肯定有效”的印象。但实际上，当你纳入那 90 项阴性研究后，总体结论可能是”该处理没有显著效应”。

这就是发表偏倚的危险之处：它让你看到的不是真实的科学图景，而是一个被扭曲的版本。

16.5.2 更多的发表偏倚检测方法

除了漏斗图和 Egger 检验之外，还有一些补充方法：

Peters 检验：Egger 检验的变体，使用更稳健的权重方案。当研究数量较少时，Peters 检验通常比 Egger 检验更可靠。

#| label: peters-test

# Peters 检验
# 注意：metafor 的 summary.peter.rma() 函数实现此检验
# 这里仅展示结构，实际运行时需要完整数据
rma(yi, vi, data = meta_es) |>
  summary(peters = TRUE)

Beggs 秩相关检验：不依赖效应量的分布假设，是一种更稳健的非参数方法。

#| label: begg-test

# Begg 检验
ranktest(res)

敏感性分析：逐一剔除法

另一种直观的方法是逐一剔除（leave-one-out）分析：每次剔除一项研究，看合并效应量是否发生显著变化。如果剔除某一项研究后，合并效应量发生大幅变化，说明该研究对总体结论有重大影响，可能需要进一步调查。

#| label: leave-one-out

# 逐一剔除敏感性分析
influence <- influence(res)
print(influence)

16.5.3 处理发表偏倚的策略

一旦检测到发表偏倚，可以采取以下策略：

1. 报告偏倚的存在：在论文中明确报告检测结果，并讨论可能的偏倚方向和影响。

2. 使用剪补法校正：剪补法可以在假设”缺失研究为阴性结果”的前提下，对合并效应量进行一定程度的校正。但要注意，剪补法的结果本身也有争议，不应过度依赖。

3. 扩大文献检索范围：主动检索灰色文献（学位论文、机构报告、会议摘要），这些通常比正式期刊更可能包含阴性结果。

4. 联系作者：尝试联系领域内的知名研究者，询问是否有未发表的相关研究。

5. 调整结论的确定性：如果存在中度至高度的发表偏倚风险，在结论中应适当降低确定性，用更保守的语言表述。

16.5.4 什么是出版偏差？

发表偏倚（Publication Bias）是荟萃分析最严重的威胁之一。

它的意思是：阳性结果（p < 0.05）的研究更容易被发表，而阴性结果（p > 0.05）的研究更容易被束之高阁。

后果是：荟萃分析纳入的研究是”有偏的样本”，合并出来的效应量会系统性地高估真实效应。

16.5.5 检测出版偏差：漏斗图

漏斗图（Funnel Plot）是最直观的出版偏差可视化工具：

# 绘制漏斗图
funnel(ma_model, backend = "metafor")

# 添加 Egger 回归线（检测对称性）
# 如果漏斗不对称，说明可能存在出版偏差

解读：

• 纵轴：效应量的精度（与样本量相关，精度 = 1/SE）
• 横轴：效应量
• 每个点：一项研究
• 理想情况：点呈倒 V 形（对称）分布，大样本研究在上方（精度高），小样本研究在下方（精度低）
• 如果不对称（缺少左下角或右下角的点），说明存在出版偏差

16.5.6 Egger 检验

Egger 检验是漏斗图不对称性的统计检验：

# Egger 检验
regtest(ma_model)

# 输出示例：
# Regression Test for Funnel Plot Asymmetry

# Model:     weighted regression
# zval =     2.345
# p-val =    0.019
# SE =       0.412
# intercept = 1.23 (95% CI: 0.21 to 2.25)

如果 p < 0.05，拒绝”漏斗对称”的原假设，说明可能存在出版偏差。

16.5.7 剪补法（Trim-and-Fill）

剪补法（Trim-and-Fill）是针对出版偏差的补救方法：

# 剪补法
taf_result <- trimfill(ma_model)
print(taf_result)

# 重新绘制漏斗图（包含估计的缺失研究）
funnel(taf_result, backend = "metafor")

剪补法的思想是：

1. 剪掉不对称部分的”多余”研究（剪）
2. 计算漏斗中心（填补）
3. 估计缺失的阴性结果研究数量（填充）
4. 重新计算合并效应量（校正后）

警告：剪补法只能给出校正效应量的估计，不能真正解决出版偏差问题。如果剪补后效应量大幅下降，说明出版偏差是严重问题，结果应谨慎解读。

16.5.8 针对出版偏差的预防措施

在荟萃分析设计阶段就应采取措施：

1. 预先注册：在 PROSPERO 或 OSF 上预先注册荟萃分析方案
2. 检索灰色文献：不要只检索期刊文章，也要检索学位论文、会议论文、技术报告
3. 联系作者：尝试联系作者获取未发表数据
4. 敏感性分析：分别分析已发表和未发表数据，比较结果差异

16.6 生态学案例：生物多样性与生产力荟萃分析

要说生态学领域最著名的荟萃分析，就不得不提 “生物多样性与生态系统功能”（Biodiversity and Ecosystem Functioning, BEF）这个研究领域。这个领域的研究试图回答一个核心问题：生物多样性越丰富，生态系统的功能（如生产力）是否就越强？

这个问题在 1990 年代引发了激烈的学术争论，而荟萃分析在其中发挥了关键的”裁判”作用。

16.6.1 经典案例：Hooper et al. (2005) 的全球荟萃分析

2005 年，Hooper 等人在 Ecological Monographs 上发表了一篇里程碑式的荟萃分析，题为 “Effects of biodiversity on ecosystem functioning: a consensus of current knowledge”。这篇论文综合了从 1990 年代到 2005 年间全球范围内开展的生物多样性实验数据，试图回答：

• 生物多样性对生态系统功能的正向效应在全球尺度上是否一致？
• 这种效应的强度有多大？
• 哪些因素调节了这个效应的大小？

这项荟萃分析纳入了来自草原、森林、海洋等多种生态系统的数百项实验，使用标准化均值差（Hedges’ d）作为效应量，合并分析后发现：生物多样性对生态系统功能（尤其是生产力）有显著的正向效应，且这个效应在控制了实验设计偏差之后依然稳健。

这篇论文的重要意义在于：它平息了当时关于 BEF 关系的激烈争论——荟萃分析的结果表明，早期单点实验发现的正向效应并不是偶然的，而是一个可以在全球尺度上重复验证的规律。

16.6.2 我们来做一个模拟演示

虽然无法复制 Hooper 等人的完整数据集，但我们可以模拟一个类似的分析流程，来感受生物多样性荟萃分析的核心思想。

假设我们整理了 12 项研究的数据，这些研究考察了植物物种多样性（物种数）对地上生物量的影响：

library(metafor)

# 模拟 12 项生物多样性实验的数据
bef_data <- data.frame(
  study = paste0("Exp", 1:12),
  ecosystem = rep(c("grassland", "forest", "marine"), c=as.numeric(table(sample(1:3, 12, replace=TRUE)))),
  n_t = sample(8:30, 12, replace=TRUE),    # 多样性高的处理组样本量
  n_c = sample(8:30, 12, replace=TRUE),    # 多样性低的对照组样本量
  mean_t = round(rnorm(12, 450, 80), 1),   # 高多样性组的生物量均值 (g/m2)
  mean_c = round(rnorm(12, 320, 70), 1),   # 低多样性组的生物量均值
  sd_t = round(rnorm(12, 80, 20), 1),
  sd_c = round(rnorm(12, 70, 15), 1)
)

bef_data

现在计算每项研究的标准化均值差（Hedges’ g）并合并：

# 计算 Hedges' g（标准化均值差）
bef_es <- escalc(
  measure = "SMD",    # Standardized Mean Difference
  m1i = mean_t,
  m2i = mean_c,
  sd1i = sd_t,
  sd2i = sd_c,
  n1i = n_t,
  n2i = n_c,
  data = bef_data,
  slab = study
)

# 拟合随机效应模型
bef_res <- rma(yi, vi, data = bef_es)
summary(bef_res)

#| fig-width: 9
#| fig-height: 7
#| label: fig-bef-forest
#| fig-cap: "植物多样性对生物量影响的森林图"

forest(
  bef_res,
  header = TRUE,
  xlab = "标准化均值差 (Hedges' g)",
  mlab = "随机效应模型合并估计",
  cex = 0.85,
  refline = 0,
  order = "size"    # 按权重（研究精度）排序
)

森林图中，菱形落在零效应线右侧，说明高多样性组的生物量显著高于低多样性组。这就是荟萃分析的魅力——它用数据说话，而不是凭感觉。

16.6.3 探索异质性来源

发现显著的正向总体效应之后，科学家们自然会问：是什么因素调节了这个效应的大小？ 比如，在草原和森林中，生物多样性效应的强度一样吗？

这就需要用到亚组分析或元回归了：

# 按生态系统类型进行亚组分析
bef_subgroup <- rma(yi, vi, mods = ~ ecosystem - 1, data = bef_es)
summary(bef_subgroup)

如果 ecosystemgrassland 和 ecosystemforest 的效应量显著不同（置信区间不重叠），就说明生态系统类型是效应量异质性的重要来源。

16.6.4 荟萃分析给你的启示

通过这个模拟案例，你应该能感受到荟萃分析在生态学中的核心价值：

1. 从局部到整体的跨越：单点实验告诉你”在我的样地里发生了什么”，荟萃分析告诉你”在全世界范围内，这个规律是否普遍成立”。

2. 量化而非猜测：荟萃分析用具体的数字（效应量大小、置信区间、I²）来描述规律，而不是模糊的”有影响”或”没有影响”。

3. 管理不确定性：通过异质性分析，你可以知道效应的大小在哪些条件下强、在哪些条件下弱，这对生态系统管理非常重要。

4. 累积证据：荟萃分析是”站在前人的肩膀上”——它把所有相关研究的结果整合起来，让后来的研究者能够系统地了解一个领域的知识现状。

16.7 生态学荟萃分析案例

16.7.1 案例：生物多样性与生产力荟萃分析

以下用一个简化的示例数据集演示完整的荟萃分析流程：

# 加载必要的包
library(metafor)
library(tidyverse)

# 示例数据：物种多样性-生产力关系的8项实验
# 数据来自虚拟场景，仅供教学演示
biodiv_data <- tibble(
  study = paste0("Study_", 1:8),
  biodiversity = c("species_richness", "species_richness", "functional_diversity",
                   "species_richness", "functional_diversity", "phylogenetic_diversity",
                   "species_richness", "functional_diversity"),
  biome = c("grassland", "grassland", "grassland", "forest",
            "forest", "grassland", "wetland", "grassland"),
  n_t = c(8, 12, 10, 15, 8, 12, 10, 14),
  mean_t = c(580, 920, 1250, 890, 760, 1100, 640, 850),
  sd_t = c(85, 130, 180, 145, 110, 155, 95, 125),
  n_c = c(8, 12, 10, 15, 8, 12, 10, 14),
  mean_c = c(420, 680, 820, 720, 590, 790, 510, 650),
  sd_c = c(78, 120, 165, 138, 102, 148, 88, 118)
)

# 计算效应量
biodiv_ma <- escalc(
  measure = "SMD",
  m1i = mean_t, m2i = mean_c,
  sd1i = sd_t, sd2i = sd_c,
  n1i = n_t, n2i = n_c,
  data = biodiv_data,
  slab = study
)

# 拟合随机效应模型
biodiv_model <- rma.uni(
  yi = yi, vi = vi,
  data = biodiv_ma,
  method = "REML",
  test = "t"
)

# 查看结果
summary(biodiv_model)

# 森林图
forest(
  biodiv_model,
  slab = biodiv_ma$study,
  showweights = TRUE,
  header = TRUE,
  order = "obs",
  xlab = "标准化效应量 (Hedges' g)"
)

# 异质性检验
confunity(biodiv_model)

# 元回归：检验不同生物多样性指标是否有不同效应
biodiv_meta_reg <- rma.uni(
  yi = yi, vi = vi,
  mods = ~ factor(biodiversity),
  data = biodiv_ma,
  method = "REML"
)

# 不同生物多样性指标的平均效应量
biodiv_meta_reg

# Egger 检验
regtest(biodiv_model)

# 结论整理
cat("\n========== 荟萃分析结论 ==========\n")
cat("合并效应量 (Hedges' g):", round(biodiv_model$b[1], 3), "\n")
cat("95% CI:", round(biodiv_model$ci.lb, 3), "~", round(biodiv_model$ci.ub, 3), "\n")
cat("p 值:", format.pval(biodiv_model$pval, digits = 3), "\n")
cat("tau²:", round(biodiv_model$tau2, 4), "\n")
cat("I²:", round(biodiv_model$I2, 1), "%\n")
cat("==================================\n")

16.7.2 报告荟萃分析的标准格式

荟萃分析的结果应该包含以下内容（参考 PRISMA 声明）：

【方法】
文献检索策略：在[数据库]检索"[检索词]"，时间截止[日期]。
纳入标准：[具体描述]
数据提取：[提取的变量和效应量类型]
统计分析：[模型类型、异质性检验方法]

【结果】
文献筛选：共检索到[X]篇，经筛选后纳入[Y]项研究
合并效应量：Hedges' g = [值]，95% CI = [下限] ~ [上限]，p = [值]
异质性检验：Q = [值]，df = [值]，p = [值]；I² = [值]%
出版偏差：[Egger 检验结果]

【结论】
[用通俗语言描述发现，避免仅重复数字]

16.8 常见误区与注意事项

16.8.1 误区一：把荟萃分析当成”高级文献综述”

荟萃分析不是”把所有研究读一遍然后写总结”，而是一种严格的统计方法。它需要：

• 预先制定明确的纳入和排除标准
• 系统检索（而不是随意选几篇文献）
• 完整提取效应量和方差
• 用统计方法合并效应量

随手选十几篇文献算个平均，不是荟萃分析。

16.8.2 误区二：忽视异质性

许多新手荟萃分析报告一个合并效应量就结束了，完全忽略了异质性。

但没有异质性的荟萃分析才是不正常的。生态学研究在不同地点、不同物种、不同时间进行，效应量有差异是必然的。

正确做法：报告 I² 和 τ²，探索异质性来源，讨论异质性对结论的影响。

16.8.3 误区三：忽视出版偏差

如果你的荟萃分析没有检验出版偏差，审稿人会直接质疑。

标准做法：至少画一个漏斗图，做一个 Egger 检验。如果发现出版偏差，使用剪补法进行敏感性分析。

16.8.4 误区四：用固定效应模型处理生态学数据

在生态学和医学领域，随机效应模型几乎总是更合理的选择。固定效应模型假设所有研究估计同一个效应，这个假设在生态学中很少成立。

规则：如果你的研究使用了不同的物种、地点、方法，用随机效应模型。

16.8.5 误区五：过度依赖统计显著性

在荟萃分析中，很多初学者会把”p < 0.05”作为判断”有没有效应”的唯一标准。但这是一个严重的误解。统计显著性受样本量影响极大——一个样本量极小的研究，即便效应量几乎为零，也可能得到 p < 0.05；一个样本量极大的研究，即便效应量有实质意义，也可能得到 p > 0.05。

荟萃分析应该以效应量及其置信区间为核心论据。p 值只是告诉你”这个效应是否显著不为零”，但它不告诉你”这个效应有多大”。你应该这样思考：合并效应量是 0.15（lnRR），95% CI 是 0.08 到 0.22——即便 p = 0.001，这个 0.15 的效应在生态学上是否有意义？0.15 意味着处理组比对照组高约 16%，这个变化在生态系统管理中是否值得关注？

置信区间给了你更丰富的信息：如果 CI 是 [0.08, 0.22]，你可以说”我们很有把握真实效应在 8% 到 22% 之间”；如果 CI 是 [-0.05, 0.35]，虽然合并效应量是 0.15，但真实效应可能低至 -5%（微小负效应）也可能高达 35%（强正效应），不确定性很大。

16.8.6 误区六：将相关性当作因果性

荟萃分析本质上是对已有研究的综合，而已有研究绝大多数都是观测性研究（而非随机对照实验）。观测性研究的固有局限是：无法控制混杂变量，因此不能直接推断因果关系。

举例来说，一项荟萃分析可能发现”土壤微生物多样性越高的样地，植物生产力也越高”，并得出”提高微生物多样性可以增加生产力”的结论。但这个结论是危险的——也许存在第三个变量（如土壤肥力）同时驱动了微生物多样性和生产力的提高，形成了一种”虚假相关”。

荟萃分析的结论应始终谨慎声明因果关系。如果没有来自随机对照实验的证据，应该使用”正相关”而非”提高 X 导致 Y 增加”这样的语言。

16.8.7 误区七：忽略非独立性

荟萃分析假设每个效应量是独立的——但这个假设在实践中经常被违反。

常见的非独立情况包括：

• 一篇论文报告了多个处理水平：如”低氮处理”和”高氮处理”来自同一篇论文，它们共享了实验设计、土壤背景和部分作者。

• 同一个实验发表在多篇论文中：研究者将同一批数据分拆成多篇论文发表，导致同一个效应被计算了两次。

• 嵌套设计：如多个样地内的多个子样方，样地之间有共享的环境条件。

忽略非独立性会导致标准误被低估，进而使置信区间过窄、显著性检验过于敏感。处理方法包括：

1. 在每个论文中只保留一个效应量：例如只保留效应最强的处理水平，或将多个处理水平取平均值。

2. 使用多层次模型：将效应量的嵌套结构（论文 > 效应量）纳入模型，同时考虑论文内和论文间的方差成分。

3. 在敏感性分析中排除非独立效应量：通过逐一剔除或子集分析来评估非独立性对结果的影响。

16.8.8 注意事项：荟萃分析的报告规范

为了保证荟萃分析报告的透明度和可重复性，国际学术界制定了 PRISMA（Preferred Reporting Items for Systematic Reviews and Meta-Analyses）声明，这是荟萃分析报告的”金标准”。

PRISMA 清单要求荟萃分析论文必须包含以下内容：

1. 标题：明确说明这是一项系统综述和/或荟萃分析

2. 摘要：结构化摘要，包括背景、目的、数据来源、纳入研究、结果、局限性、结论

3. 引言：研究问题的理论基础，综述的目的或假设

4. 方法：检索策略（数据库、检索词、时间范围）、纳入排除标准、研究筛选流程、数据提取方法、效应量选择、统计分析方法

5. 结果：文献筛选流程图（PRISMA flow diagram）、纳入研究的基本特征、合并效应量及其置信区间、异质性检验结果、发表偏倚检测结果、敏感性分析结果

6. 讨论：主要发现的总结、与已有研究的比较、局限性、结论

7. 资金来源：声明资金资助和利益冲突

建议：在你开始做荟萃分析之前，就把 PRISMA 清单打印出来，逐项对照检查，确保没有遗漏。

16.8.9 注意事项：效应量的可解释性

Cohen’s d 或 Hedges’ g 是标准化效应量，它们没有单位。有时审稿人会质疑”这个 d = 0.5 到底意味着什么”。

解释方法：

1. 与经验标准对比：Cohen 的经验标准（0.2 = 小，0.5 = 中，0.8 = 大）
2. 与领域经验对比：查阅该领域已发表的荟萃分析，比较效应量范围
3. 换算为原始单位：在讨论部分，将标准化效应量换算为原始变量的变化幅度

16.9 课后练习

16.9.1 基础题

练习 1：效应量计算

某荟萃分析纳入两项关于氮添加对土壤有机碳影响的研究。研究数据如下：

• 研究 A：处理组 SOC = 28.5 g/kg (SD = 4.2, n = 15)，对照组 SOC = 24.1 g/kg (SD = 3.8, n = 15)
• 研究 B：处理组 SOC = 32.1 g/kg (SD = 5.8, n = 8)，对照组 SOC = 27.3 g/kg (SD = 4.9, n = 8)

请计算：

1. 每项研究的 Cohen’s d
2. 每项研究的 Hedges’ g
3. 哪个效应量更保守？为什么？

练习 2：模型选择判断

以下场景应该用固定效应模型还是随机效应模型？说明理由：

1. 6 项研究都在同一个实验站进行，使用完全相同的protocol，只有土壤类型略有不同
2. 12 项研究来自全球不同地点（欧洲、北美、亚洲），使用了不同的采样方法和测量技术
3. 3 项研究由同一个实验室在不同年份完成，使用相同的设备和操作流程

练习 3：解读荟萃分析结果

某荟萃分析报告了以下结果，请判断每个陈述是否正确：

随机效应模型 (k = 15)
tau² = 0.45, I² = 72.3%
合并效应量 = 0.68, 95% CI [0.41, 0.95], p < 0.001
Egger 检验：z = 1.23, p = 0.22

判断题：

1. 72.3% 的总变异来自研究间异质性
2. 合并效应量在统计上显著大于 0
3. Egger 检验说明存在严重的出版偏差
4. tau² = 0.45 说明研究间方差较大

16.9.2 应用题

练习 4：完成一个小型荟萃分析

从以下场景出发，完成一个完整的荟萃分析练习（使用 R 和 metafor 包）：

场景：你对”土壤动物多样性对有机质分解速率的影响”感兴趣，在文献检索中找到 5 项相关研究，请完成以下步骤：

1. 创建包含以下信息的数据框（可用虚拟数据）

study	mean_t	sd_t	n_t	mean_c	sd_c	n_c
Study1	0.85	0.12	12	0.72	0.10	12
Study2	0.91	0.15	8	0.75	0.11	8
Study3	0.78	0.09	15	0.68	0.08	15
Study4	0.88	0.14	10	0.70	0.12	10
Study5	0.82	0.11	11	0.73	0.09	11

2. 使用 escalc() 计算 Hedges’ g 和方差
3. 拟合随机效应模型
4. 解读 tau²、I² 和合并效应量
5. 绘制森林图
6. 进行 Egger 检验并解释结果

练习 5：异质性来源分析

使用上题的数据，探索异质性来源：

1. 假设 Study 1 和 Study 2 来自”热带”生态系统，Study 3-5 来自”亚热带”生态系统
2. 用元回归检验生态系统类型是否解释效应量的变异
3. 分别计算热带和亚热带研究的平均效应量

16.9.3 思考题

练习 6：荟萃分析的局限性

荟萃分析是一种强大的方法，但它也有局限性。请讨论以下问题：

1. 如果你想研究”马尾松混交林对土壤碳储量的影响”，荟萃分析是否适合？需要考虑哪些问题？

2. 如果某领域的已发表文献只有 3-4 项研究，还适合做荟萃分析吗？应该怎么办？

3. 荟萃分析的结果能用来指导具体的管理决策吗？为什么或为什么不？

练习 7：搜索一项真实的生态学荟萃分析

1. 在 Google Scholar 中检索一项生态学领域的荟萃分析（关键词如 “meta-analysis biodiversity ecosystem functioning”）
2. 阅读摘要，判断它的问题是什么、用了哪种效应量、纳入了多少研究
3. 查看全文（或 methods 部分），判断它是否进行了出版偏差检验
4. 评价这项荟萃分析的质量（优缺点）

16.9.4 进阶练习

练习 8：解读荟萃分析的真实论文

请在 Google Scholar 中搜索以下论文并完成分析：

1. “Global meta-analysis of soil microbial biomass responses to organic amendments” (Sime et al., 2020, Soil Biology and Biochemistry)

2. 回答以下问题：

   • 这项荟萃分析的研究问题是什么？
   • 纳入了多少项研究？覆盖了多少个国家？
   • 使用了哪种效应量？为什么？
   • 合并后的总体效应是什么？（方向和大小）
   • 异质性如何？哪些因素解释了效应量的变异？
   • 是否进行了发表偏倚检验？结果如何？
   • 你认为这项荟萃分析的主要局限性是什么？

练习 9：设计自己的荟萃分析

假设你想做一项关于”凋落物去除对土壤碳通量影响”的荟萃分析，请回答：

1. 你会如何制定检索策略？列出至少 5 个相关数据库和 10 个检索词组合。

2. 你的纳入和排除标准是什么？（至少列出 5 条纳入标准和 5 条排除标准）

3. 你会提取哪些变量？设计一个数据提取表格。

4. 你会选择哪种效应量？理由是什么？

5. 你计划如何检测和处理发表偏倚？

6. 你的亚组分析计划包括哪些调节变量？

16.9.5 实践练习

练习 10：完成一个真实的荟萃分析小练习

本练习使用 R 的 metafor 包和内置的教学数据集，带你走一遍完整的荟萃分析流程。

数据说明：我们将使用 metafor 包内置的 dat.bornmann2007 数据集，该数据集来自 Bornmann et al. (2007) 关于同行评审有效性的荟萃分析，共包含 33 项研究。

1. 加载数据并查看结构

library(metafor)

# 加载内置数据集
data("dat.bornmann2007")

# 查看数据维度
dim(dat.bornmann2007)

# 查看前几行
head(dat.bornmann2007[, 1:8])

2. 计算效应量（如果数据集中没有的话）

数据集已包含效应量（yi）和方差（vi）列。检查这些列的含义：

# 查看效应量的分布
hist(dat.bornmann2007$yi, breaks = 15,
     main = "效应量分布", xlab = "Hedges' g")

3. 拟合随机效应模型并解释结果

# 拟合随机效应模型
res <- rma(yi, vi, data = dat.bornmann2007)
summary(res)

回答： - 合并效应量是多少？95% CI 是多少？ - p 值是多少？效应是否显著？ - tau² 和 I² 是多少？异质性如何？ - Egger 检验结果如何？

4. 绘制森林图

#| fig-width: 8
#| fig-height: 12

forest(
  res,
  header = TRUE,
  slab = paste(dat.bornmann2007$author, dat.bornmann2007$year),
  xlab = "标准化均值差 (Hedges' g)",
  cex = 0.7
)

5. 绘制漏斗图并进行 Egger 检验

#| fig-width: 6
#| fig-height: 5

funnel(res, main = "漏斗图")

# Egger 检验
regtest(res)